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

Question

我想用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 函数的结果进行比较，但我的代码给出了错误的结果。你能看看我的代码并给我一些意见吗？

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)

Answer 1

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

获得 num_generated 泊松结果的有效实现应该类似于以下伪代码：

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
}

Answer 2

您不知道 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 个，您对这个变量执行的任何操作都会容易得多。

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

How to generate random number that satisfying poisson distribution

algorithm

statistics

matlab

probability