产生超几何分布的随机变量答案

【问题标题】：Producing random variates distributed hypergeometrically产生超几何分布的随机变量
【发布时间】：2019-04-11 06:45:35
【问题描述】：

我们如何在 C++ 中为超几何分布生成随机变量？不幸的是，标准库不提供std::hypergeometric_distribution。有一个boost::math::hypergeometric_distribution；但是，它似乎不能用于生成随机变量。

【问题讨论】：

标签： c++ random boost distribution

【解决方案1】：

好吧，这里是简单的代码（可能有错误，尽管它编译、运行并产生了与事实相似的东西）——采样Hypergeometric distribution。符号与 Wiki 文章中的相同。在内部，它用值 [0...N-1] 填充数组，首先打乱它们中的n，如果打乱的值小于K（成功），则计数。

Shuffle 是经典的 Fisher-Yates-Knuth 之一。

#include <cstdint>
#include <iostream>
#include <random>
#include <vector>
#include <numeric>
#include <cmath>

#define func auto


func nCk(uint64_t n, uint64_t k) -> double { // computes binomial coefficient
    if (k > n)
        return 0.0;

    if (k * 2ULL > n)
        k = n - k;

    if (k == 0.0)
        return 1.0;

    return exp(lgamma(n+1) - (lgamma(k+1) + lgamma(n-k+1)));
}


func PMF(uint64_t N, uint64_t K, uint64_t n, uint64_t k) -> double { // Hypergeometric distribution PMF
    return nCk(K, k) * nCk(N - K, n - k) / nCk(N, n);
}


func sample_hyper(int N, int K, int n,            // sampling from Hypergeometric distribution
                  std::mt19937_64&  rng,
                  std::vector<int>& nums) -> int {

    int rc{ 0 };
    for (int k = 0; k != n; ++k) {
        std::uniform_int_distribution<int> uni{ k, N-1 };
        auto s = uni(rng);
        std::swap(nums[k], nums[s]);

        rc += (nums[k] < K);
    }
    return rc;
}


func main() -> int {
    auto rng = std::mt19937_64{1234567891ULL};

    int N = 500; // compare with Wiki
    int K = 50;
    int n = 100;

    auto nums = std::vector<int>( N );
    std::iota(nums.begin(), nums.end(), 0); // array to shuffle, filled with 0, 1, 2, 3 ... sequence

    auto h = std::vector<int>( K, 0 ); // histogram

    int NT = 100000; // number of trials
    for(int k = 0; k != NT; ++k) {
        int q = sample_hyper(N, K, n, rng, nums);
        h[q] += 1; // count it
    }

    for (int k = 0; k != 20; ++k) { // and print it
        double e = double(h[k]) / double(NT);
        std::cout << k << "   " << e << "     " << PMF(N, K, n, k) << '\n';
    }

    return 0;
}

更新

我已经为超几何分布实现了 PMF，我的采样似乎与它一致。请检查更新的代码。

【讨论】：