用均勻分布隨機變量生成泊松分布隨機變量

　　《R語言的科學編程與仿真》的第18章提到，所有的隨機變量可以通過處理U(0,1)隨機變量生成。該書在18.2里給出了一個模擬算法，具體內容摘抄如下：

　　假設X是在集合{0,1，...}取值的離散隨機變量，累積分布函數是F，概率質量函數是p。下面一段代碼給定一個均勻隨機變量U並返回一個累積分布函數是F的離散隨機變量X。

# given U~U(0,1)
X<-0
while(F(X)<U){
  X<-X+1
}

 　　當這個算法結束時，我們有F(X)≥U和F(X-1)<U，因此

　　P{X=x}=P{F(x-1)<U≤F(x)}=F(x)-F(x-1)=p(x)

正是所需要的。

　　基於上述想法，這里用均勻分布來生成泊松分布，並畫圖比較。首先是寫一個函數來產生由均勻分布生成的服從泊松分布的隨機變量的取值，其次是進行重復試驗，運用頻率來估計隨機變量取該值的概率，最后是畫圖進行比較。

#第一步，生成模擬函數
poiss_sim<-function(lambda){
  x<-0 # 初始值定為0
  px<-exp(-lambda)
  Fx<-px
  U<-runif(1)
  while(Fx<U){
    x<-x+1
    px<-px*lambda/x # 利用泊松分布的概率密度函數p(x)和p(x-1)存在的遞歸關系，簡化運算
    Fx<-Fx+px # Fx即為分布函數F(x) 
  }
  return(x)
}

#第二步，重復試驗，這里的實驗次數N定為10000次
N<-1e4
lambda<-2

x<-replicate(N,poiss_sim(lambda=lambda))
y<-table(x)
names<-as.numeric(names(y))
freq<-as.numeric(y)
phat<-vector(length=length(names))
for(i in 1:length(names)){
  phat[i]<-sum(x==names[i])/N
} # 用頻率估計概率

#第三步，和R自帶的函數dpois()進行圖形比較
par(las=1)
plot(names,dpois(names,lambda=lambda),xlab='x',ylab='p(x)',type='h')
points(names,dpois(names,lambda=lambda),pch=19)
points(names,phat,pch=3)
legend(6,0.25,legend=c('reality','simulation'),pch=c(19,3),bty='n',cex=2)

 　　圖形如下：

　　可以看出，真實值和模擬值還是很接近的！