神奇的卜瓦松分配(Poisson distribution)

過去在寫亂數的程式時，有一種情形一直讓我很困擾，就是希望可以隨機產一串數字後，這些數字的平均值要接近一個預期的平均值。可是，經常會發生結果與預期的平均值相差太多的情形。過去使用的方法是均勻分配(uniform distribution)的作法，例如預期平均每λ時間單位會發生一次事件，就以1/λ當作每個時間發生事件的機率，接著就在每個時間產生一個介於0到1之間的亂數，當產生的亂數小於機率值時就視為事件發生，否則就是事件未發生。當模擬時間太短時，此種作法的執行結果和預期平均值常會有很大的差距。後來看到很多文獻都有提到卜瓦松分配，所以就試著以卜瓦松分配的概念來產生亂數，這樣的方法比較不會有結果和預期平均值的相差太多的情形。詳細的方法如下：
1. 定義一個卜瓦松分配的有效區間。（本人是設定為λ*3，這是為了寫程式方便增加此區間，原卜瓦松分配沒有定義這個區間）
2. 依據卜瓦松分配，計算有效區間內的每一個機率質量函數的值。
3. 計算機率質量函數的以下累積值。
4. 將時間設為0。
5. 產生一個0到1之間的亂數，尋找和此亂數最接近的以下累積值，此累積值的對應變數就是到下一個事件所需的時間。
6. 調整時間到下一個事件發生的時間。
7. 重覆5和6一直到時間超過模擬時間。
以下是使用java撰寫的完整程式，包含了均勻分配及卜瓦松分配的比較。

import java.util.Random;
public class compare {
  public static void main(String[] args) {
    int lambda=20;
    Random r = new Random();
    System.out.println("lambda=="+lambda);
    int simulation_time;
    int count;
    double p=(double)1/(double)lambda;
    double[] d;
    d=new double[lambda*3];
    for(int i=0;i < lambda*3; i++) {
      d[i]=Math.exp(-lambda)*Math.pow(lambda, i);
      for(int j=1;j<=i;j++) {
        d[i]=d[i]/j;
      }
    }
    for(int i=1;i < lambda*3;i++) {
      d[i]=d[i-1]+d[i];
    }
    for(simulation_time=20;simulation_time<=200;simulation_time+=20) {
      System.out.println("simulation time=="+Integer.toString(simulation_time));
      count=0;
      for(int i=0;i<=simulation_time;i++) {
        double x=r.nextDouble();
        if(x < p)
          count++;
      }
      float ave=(float)simulation_time/(float)count;
      System.out.println("simple random average=="+Float.toString(ave));
      count=0;
      int t=0;
      while(t        int next=0;
        double rn=r.nextDouble();
        while(rn>=d[next] && next          next++;
        }
        if(d[next]-rn>rn-d[next-1])
          t=t+next-1;
        else
          t=t+next;
        if(t          count++;    
      }
      double difference=Math.abs(lambda-ave);
      ave=(float)simulation_time/(float)count;
      System.out.println("Poisson random average=="+Float.toString(ave));
      difference-=Math.abs(lambda-ave);
      if(difference>=0)
        System.out.println("Poisson is better");     
      else
        System.out.println("uniform is better");     
    }
  }
}