《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 names){
  phat[i+1]<-sum(x==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)

　　圖形如下：

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