C++ Boost用泊松-伽马混合生成实参负二项随机数的错误排查

你的判断完全准确，错误就出在泊松分布的初始化语句上：

• boost::random::poisson_distribution<int> my_poi(my_gamma(rng)); 这行代码仅在分布初始化阶段调用了1次gamma采样，把这一次生成的随机值当成了泊松分布固定不变的lambda参数，后续100万次泊松采样全程使用这个固定值，完全没有实现“每个样本都先抽gamma再抽泊松”的混合逻辑。
• 按你传入的参数计算，初始化时抽到的gamma值恰好落在0.5附近，因此最终模拟均值稳定在0.5，和你观测到的结果完全吻合。

用泊松-伽马混合生成实r值负二项分布的核心要求是：每生成1个最终的负二项样本，都要先从伽马分布抽取1个随机lambda，再用这个lambda作为参数抽取泊松样本，不能提前把lambda固定死。
你之前的参数推导是完全正确的：对应负二项分布（均值r*(1-p)/p、方差r*(1-p)/p²）的伽马分布参数为shape=r、scale=(1-p)/p，和你代码里的计算一致。
另外原代码用int类型存储采样总和存在溢出风险，大样本量下建议改用long long类型。

int nr_sim = 1000000;
double mean = 2.0;
double variance = 15.0;
double r = mean * mean / (variance - mean);
double p = mean / variance;
double beta = (1 - p) / p;

typedef boost::mt19937 RNGType;
RNGType rng(5);

boost::random::gamma_distribution<double> my_gamma(r, beta);
boost::random::poisson_distribution<int> my_poi;
long long simulated_sum = 0;
for (int i = 0; i < nr_sim; i++) {
    double current_lambda = my_gamma(rng);
    simulated_sum += my_poi(rng, boost::random::poisson_distribution<int>::param_type(current_lambda));
}
double my_result = static_cast<double>(simulated_sum) / nr_sim;

运行上述代码，得到的模拟均值会稳定在2.0左右，和理论值一致。如果不想手动传param_type，也可以在每次循环内临时构造泊松分布对象，只是运行效率会比复用分布对象、动态传参稍低。

内容的提问来源于stack exchange，提问作者dan