从低均值泊松分布中绘制数字的性能

2024-03-17

为了在C++中从泊松分布中抽取随机数,通常建议使用

RNG_type rng;
std::poisson_distribution<size_t> d(1e-6);
auto r = d(rng);

每次呼叫时std::poisson_distribution对象,整个随机位序列被消耗(例如 32 位标准::mt19937 http://www.cplusplus.com/reference/random/mt19937/, 64 位为std::mt19937_64 http://www.cplusplus.com/reference/random/mt19937_64/)。让我惊讶的是,如此低的平均值(mean = 1e-6),绝大多数时候,只需几个位就足以确定要返回的值为 0。然后可以缓存其他位以供以后使用。

假设设置为 true 的位序列与泊松分布的高返回值相关联,当使用以下平均值时1e-6,任何不以 19 个 true 开头的序列都必然返回零!的确,

1 - 1/2^19 < P(0, 1e-6) < 1 - 1/2^20

, where P(n, r)表示抽签的概率n来自均值的泊松分布r。不浪费位的算法将在一半的时间使用一位,四分之一的时间使用两位,八分之一的时间使用三位,......

是否有一种算法可以通过在绘制泊松数时消耗尽可能少的位来提高性能?与其他方法相比,是否有其他方法可以提高性能std::poisson_distribution当我们考虑低均值时?


回应@Jarod42 的评论,他说

想知道使用更少的位数是否不会破坏等概率......

我不认为这会破坏等概率。在一次模糊的测试中,我用简单的伯努利分布考虑了同样的问题。我以概率抽样真实1/2^4并以概率采样错误1 - 1/2^4。功能drawWithoutWastingBits一旦在缓存和函数中看到 true 就停止drawWastingBits无论这些位是什么,都会消耗 4 位。

#include <iostream>
#include <vector>
#include <string>
#include <algorithm>
#include <random>

bool drawWithoutWastingBits(std::vector<bool>& cache, size_t& cache_index)
{
    /* 
        Get a true with probability 1/2^4 (=1/16=0.0625) and a false otherwise
    */

    size_t nbTrues = 0;
    while (cache[cache_index])
    {
        ++nbTrues;
        ++cache_index;
        if (nbTrues == 4)
        {
            return true;
        }
    }
    ++cache_index;
    return false;
}


bool drawWastingBits(std::vector<bool>& cache, size_t& cache_index)
{
    /* 
        Get a true with probability 1/2^4 (=1/16=0.0625) and a false otherwise
    */

    bool isAnyTrue = false;
    for (size_t i = 0 ; i < 4; ++i)
    {
        if (cache[cache_index])
        {
            isAnyTrue = true;
        }
        ++cache_index;
    }
    return !isAnyTrue;
}

int main()
{
    /*
        Just cache a lot of bits in advance in `cache`. The same sequence of bits will be used by both function.
        I am just caching way enough bits to make sure they don't run out of bits below
        I made sure to have the same number of zeros and ones so that any deviation is caused by the methodology and not by the RNG
    */

    // Produce cache
    std::vector<bool> cache;
    size_t nbBitsToCache = 1e7;
    cache.reserve(nbBitsToCache);
    for (size_t i = 0 ; i < nbBitsToCache/2 ; ++i)
    {
        cache.push_back(false);
        cache.push_back(true);
    }
    // Shuffle cache
    {
        std::mt19937 mt(std::random_device{}());
        std::shuffle(cache.begin(), cache.end(), mt);
    }


    // Draw without wasting bits
    {
        size_t nbDraws = 1e6;
        size_t cache_index = 0;
        std::pair<size_t, size_t> outcomes = {0,0};
        for (size_t r = 0 ; r < nbDraws ; ++r)
        {
            drawWithoutWastingBits(cache, cache_index) ? ++outcomes.first : ++outcomes.second;
            assert(cache_index <= cache.size());
        }   

        assert(outcomes.first + outcomes.second == nbDraws);
        std::cout << "Draw Without Wasting Bits: prob true = " << (double)outcomes.first / nbDraws << "\n";
    }


    // Draw wasting bits
    {
        size_t nbDraws = 1e6;
        size_t cache_index = 0;
        std::pair<size_t, size_t> outcomes = {0,0};
        for (size_t r = 0 ; r < nbDraws ; ++r)
        {
            drawWastingBits(cache, cache_index) ? ++outcomes.first : ++outcomes.second;
            assert(cache_index <= cache.size());
        }   

        assert(outcomes.first + outcomes.second == nbDraws);
        std::cout << "Draw Wit Wasting Bits: prob true = " << (double)outcomes.first / nbDraws << "\n";
    }
}

可能的输出

Draw Without Wasting Bits: prob true = 0.062832
Draw Wit Wasting Bits: prob true = 0.062363

德夫罗耶的非均匀随机变量生成 http://luc.devroye.org/rnbookindex.html,第 505 页和第 86 页,提到了顺序搜索算法的反演。

根据该算法,如果您知道mean远小于 1,那么如果生成均匀随机变量u在 [0, 1] 中,泊松变量将为 0,如果u <= exp(-mean),否则大于 0。

如果平均值较低并且您可以容忍近似分布,那么您可以使用以下方法(参见“附录 A”)差分隐私的离散高斯 https://arxiv.org/pdf/2004.00010.pdf"):

  1. Express mean以有理数的形式,以有理数的形式numer/denom。例如,如果mean是一个固定值,那么numer and denom可以相应地预先计算,例如在编译时。
  2. 随机生成伯努利(numer / denom) 数字(生成 1 的概率numer / denom或 0 否则)。如果 1 是这样生成的,请使用伯努利重复此步骤(numer / (denom * 2)), 伯努利(numer / (denom * 3)),依此类推,直到这样生成0。使用一种可以最大限度地减少比特浪费的算法生成这些数字,例如 Lumbroso 的 Fast Dice Roller 论文(2013 年)附录 B 中提到的算法,或者从那里修改并在我的章节中给出的“ZeroToOne”方法布尔条件 https://github.com/peteroupc/peteroupc.github.io/blob/master/randomfunc.md#boolean-truefalse-conditions。也可以看看这个问题 https://stackoverflow.com/questions/60777414/uniformly-distributed-bit-sequence.
  3. 如果步骤 2 产生偶数个 1,则泊松变量恰好为 0。
  4. 如果步骤 2 产生奇数个 1,则泊松变量大于 0,并且需要“较慢”的算法,仅对大于 0 的泊松变量进行采样。

例如,假设平均值为 1e-6 (1/1000000),则生成伯努利 (1/1000000) 数,然后生成伯努利 (1/2000000) 等。直到以这种方式生成 0。如果生成偶数个,则泊松变量恰好为 0。否则,泊松变量为 1 或更大,并且需要“较慢”的算法。

下面的算法就是一个例子,它基于第 505 页和第 86 页中的算法,但仅对泊松变量 1 或更大的样本进行采样:

METHOD Poisson1OrGreater(mean)
 sum=Math.exp(-mean)
 prod=sum
 u=RNDRANGE(sum, 1)
 i=0
 while i==0 or u>sum
   prod*=mean/(i+1)
   sum+=prod
   i=i+1
 end
 return i
END METHOD

不过,这种方法不是很稳健,特别是因为它使用接近 1 的数字(浮点空间更稀疏)而不是接近 0 的数字。


请注意,总和为n独立泊松(mean) 随机变量是泊松(mean*n)分发(第 501 页)。因此,本答案中的上述讨论适用于n泊松随机变量只要n有时它们的平均值仍然很小。例如,要生成均值为 1e-6 的 1000 个泊松随机变量的总和,只需生成均值为 0.001 的单个泊松随机变量。这将大大节省对伪随机数生成器的调用。


还有另一种方法可以生成低均值(1 或更少)的泊松变量。 Duchon 和 Duvignau 在“在增长的均匀排列中保留长度为 k 的循环数”中对此进行了描述,组合学电子期刊 23(4),2016。

首先,生成一个 Poisson(1) 随机变量x = Poisson1()使用下面给出的算法,该算法仅使用整数算术(其中RNDINT(a)生成 [0,a]):

METHOD Poisson1()
  ret=1; a=1; b=0
  while true // until this method returns
    j=RNDINT(a)
    if j<a and j<b: return ret
    if j==a: ret=ret+1
    else
      ret=ret-1; b=a+1
    end
    a=a+1
  end
END METHOD

Now let mean是所需的平均值。抛硬币x次,其中硬币正面朝上的概率等于mean。 (换句话说,生成一个二项式(x, mean)随机变量。)正面的数量是泊松(mean)随机变量。

本文内容由网友自发贡献,版权归原作者所有,本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容,请联系:hwhale#tublm.com(使用前将#替换为@)

从低均值泊松分布中绘制数字的性能 的相关文章

随机推荐