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

Question

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

1. 生成统一的 r.v. z1 , z2 , ...
2. 在 z1.z2..zm<=exp(-lamda*T) 时停止
3. 分配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

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