如何生成满足泊松分布的随机数答案

【问题标题】：How to generate random number that satisfying poisson distribution如何生成满足泊松分布的随机数
【发布时间】：2016-01-31 20:27:21
【问题描述】：

我想用lambda = 1和T=6生成500000个泊松分布的随机数，使用组合方法可以描述如下：

生成统一的 r.v. z1, z2, ...
z1.z2..zm<=exp(-lamda*T)时停止
分配k = m – 1

然后计算 10 个区间（[0,1],[2,3],...,[16,17],[18,∞)] 中每个区间有多少个数。

我知道 MATLAB 有一个内置函数 poissrnd 用于上述任务。但是，我想用上面的算法自己做。我尝试这样做并将其与poissrnd 函数的结果进行比较，但我的代码给出了错误的结果。你能看看我的代码，给我一些cmets吗？

num_generated = 500000;
lambda=1;T=6;
k_vec=[]; %% Store k
for i=1:number_generated
    multiple=1;
    for j=1:number_generated
        %% Step 1: Generate uniform in the interval [0,1]: z1,z2...
        z=rand(); 
        %% Step 2: Stop when z1z2...zm<=exp(-lambda*T)
        multiple=multiple*z;
        if(multiple<=exp(-lambda*T))
            k=j-1;
            k_vec=[k_vec k]; % Record k in vec
            break;
        end
    end
end
range_1 = sum( k_vec(:)==0 )+sum(k_vec(:)==1) % # number with in range [0,1]
range_2 = sum( k_vec(:)==2 )+sum( k_vec(:)==3) % # number with in range [2,3]
range_3 = sum( k_vec(:)==4 )+sum( k_vec(:)==5) % # number with in range [4,5]
range_4 = sum( k_vec(:)==6 )+sum( k_vec(:)==7) % # number with in range [6,7]
range_5 = sum( k_vec(:)==8 )+sum( k_vec(:)==9) % # number with in range [8,9]
range_6 = sum( k_vec(:)==10 )+sum( k_vec(:)==11) % # number with in range [10,11]
range_7 = sum( k_vec(:)==12 )+sum( k_vec(:)==13) % # number with in range [12,13]
range_8 = sum( k_vec(:)==14 )+sum( k_vec(:)==15) % # number with in range [14,15]
range_9 = sum( k_vec(:)==16 )+sum( k_vec(:)==17) % # number with in range [16,17]
range_10 = sum(k_vec(:)>=18)         % # number with in range [18,+infty)

【问题讨论】：

标签： algorithm matlab statistics probability

【解决方案1】：

您不知道multiple 需要多少随机值才能收敛，因此您需要将for 循环在j 上更改为一个while 循环，该循环将持续到multiple > exp(-lambda*T)。

通过将其更改为 while 循环，您现在需要 k 作为计数器并在循环的每次迭代中递增：

（警告：未经测试的代码）

for i = 1:number_generated
    multiple = 1;
    k = 0;   %// Initialize counter for each number generated
    while multiple > exp(-lambda*T)   %// replace `for` loop
        k = k + 1;    %// Increment counter
        %% Step 1: Generate uniform in the interval [0,1]: z1,z2...
        z = rand(); 
        %% Step 2: Stop when z1z2...zm<=exp(-lambda*T)
        multiple = multiple*z;
    end
    %// If we exit the loop, we know multiple <= exp(-lambda*T)
    k = k - 1;
    k_vec = [k_vec k]; % Record k in vec
end

您还应该不惜一切代价避免使用诸如range_1、range_2、...等顺序变量名称... Matlab 旨在处理数组和矩阵，因此您应该使用它们。在您的情况下，最简单的方法是：

range(1) = sum(...
range(2) = sum(...
...
range(10) = sum(...

现在您的工作区中有一个变量，而不是 10 个，您对这个变量执行的任何操作都会变得更加容易。

【讨论】：

【解决方案2】：

我不使用 Matlab，因此无法为您提供修复的确切语法。至少，您似乎忘记为每个新泊松重置 multiple 和 k。此外，您只生成一个z。

获得num_generated Poisson 结果的工作实现应该类似于以下伪代码：

threshold = Math.exp(-lambda * T)
loop num_generated times {
    %% Each time through this loop produces a single Poisson outcome 
    count = 0
    product = 1.0
    while (product = product * rand()) >= threshold {
      count += 1
    }
    %% count now has a valid Poisson value, do what you want with it
}

【讨论】：