離散型隨機變數的模擬－逆變換法

01-24

離散型隨機變數一般用分布律來描述其分布，即用隨機變數的取值及其取值的概率來描述。比如貝努力分布、二項分布、泊松分布、幾何分布等。服從這些分布的隨機變數都可以用R軟體進行隨機模擬。

比如二項分布的隨機變數的模擬函數為rbinom(n,size,prob），其中為生成隨機數的個數，（size，prob)為而二項分布的參數。假設二項分布為 $B(n,p)$ ,那麼n＝size，prob＝p。泊松分布隨機變數模擬函數為rpois(n,lambda），幾何分布隨機變數模擬函數為rgeom(n, prob)。

對於這些已經熟知的分布，我們不需要設計演算法進行實現，直接拿去用就可以了。但是對於那些在R軟體中找不到模擬函數的分布，我們必須找到相應的演算法給予實現。在隨機變數的模擬中，我們給出了逆變換演算法，即只要給出分布函數 $F$ ，我們通過計算分布函數的逆函數 $F^{-1}$ ,就可以得到服從 $F$ 的隨機變數 $F^{-1}(U)$ ,其中U是隨機數。該方法可以用於離散型隨機變數的模擬。下面我們給出具體的離散型隨機變數的模擬演算法，證明我這裡就不再給出，讀者可以自行證明。

假設一個離散型分布為

那麼模擬的隨機變數為

演算法為

輸入：F

所用的隨機數： $Usim U[0,1]$

輸出： $Xsim F$

步驟：

（1）生成隨機數 $U$ ;

（2）若 $F(x_{i-1})leq U < F(x_i)$ ,則 $X=x_i$

例1 假設一離散性分布取值為1，2，3，4.其取值概率為

$P(X=1)=0.2,P(X=2)=0.15,P(X=3)=0.25,P(X=4)=0.4$ ,請模擬一隨機變數使得其分布為上述分布。

演算法步驟：

(1) 生成隨機數 $U$

(2) 如果 $U<0.2,$ 則令 $X=1$ ;否則

如果 $U<0.35,$ 則令 $X=2;$ 否則

如果 $U<0.6,$ 則令 $X=3$ ；否則

令 $X=4$ 。

我們按照上述演算法步驟用R語言生成10個隨機數

x<-rep(1,10)for (i in 1:10){ u=runif(1) if (u<0.2) x[i]=1 else if(u<0.35) x[i]=2 else if(u<0.6) x[i]=3 else x[i]=4 }x

運行結果為

[1] 4 4 3 4 4 2 1 3 2 1

如何檢驗所模擬的隨機變數服從給定的分布？如果給定的分布為離散型分布，可以考慮用 $chi^2$ 檢驗。下面我們生成100個隨機數，檢驗它們是否服從給定的分布。

x<-rep(1,100)p<-c(0.2,0.15,0.25,0.4)for (i in 1:100){ u=runif(1) if (u<0.2) x[i]=1 else if(u<0.35) x[i]=2 else if(u<0.6) x[i]=3 else x[i]=4 }chisq.test(table(x),p)

檢驗的結果如下：

Pearson"s Chi-squared test

data: table(x) and p

X-squared = 12, df = 9, p-value = 0.2133

由於p值為0.2133>0.05，因此接受所生成的隨機數為給定的分布。

當然我們也可以利用R中的sample函數來實現。

x<-1:4p<-c(0.2,0.15,0.25,0.4)sample(x,10,replace=T,p)

假設我們打算模擬一個有限狀態空間上服從均勻分布的隨機變數，利用上面的方法是可以的實現的。事實上對逆變換法稍作改進，就可以得到一個非常簡單的模擬方法。具體推導和做法如下：

假設一個隨機變數 $X$ 的狀態空間為 ${1,2,cdots,n}$ ,且 $P{X=j}=frac{1}{n}, i=1,2,cdots,n.$ 易知 $X$ 的分布函數為 $F(j)=frac{j}{n}, j=1,2,...,n$ .那麼由前面給出演算法可知

該式等價於 $X=j$ if $j-1leq nU<j$ ;換句話說為

其中Int(x)表示x的整數部分(即小於或者等於x的最大整數)。

例子2 假設 $X$ 服從 ${1,2,cdots,10}$ 上的均勻分布，請用R語言模擬100個 $X$ 的隨機數。

R語言實現如下：

u<-runif(100)x<-floor(10*u)+1x

運行結果如下：

> x

[1] 7 10 4 9 2 2 7 8 10 2 6 2 9 6 3 7 3 2 5 9 1 3 10 10 9 4 1 3 1 6 10 1 4 9 5 6 8 8

[39] 8 2 10 10 6 3 9 2 4 9 2 3 1 9 8 10 1 4 6 8 10 9 6 7 1 7 9 2 9 5 1 7 6 2 4 4 10 8

[77] 7 3 3 3 7 8 7 6 9 6 3 3 7 6 8 2 9 4 10 4 10 10 3 3

作業1 用逆變換法模擬下面隨機變數 $X$ ，其分布為

$P{X=1}=0.3,P{X=2}=0.2$ , $P{X=3}=0.35,P{X=4}=0.15$ ，用 $chi^2$ 檢驗模擬的效果。

參考答案：

invmo<-function(n){#產生n個隨機數； x<-c() for (i in 1:n){ u<-runif(1) if (u<=0.3) x<-c(x,1) else if(u<=0.5) x<-c(x,2) else if(u<=0.85) x<-c(x,3) else x<-c(x,4) } x }

> invmo(100)#輸出100個服從上述分布的隨機數

[1] 3 1 1 3 3 2 1 2 1 3 1 1 3 3 3 1 3 3 2 2 2 1 4 3 3 1 1 3 2 1 2 4 2 2

[35] 2 4 1 3 3 3 2 4 4 2 4 1 2 3 1 3 1 1 1 1 2 1 2 1 1 3 3 3 2 3 3 3 3 4

[69] 1 1 1 3 2 3 1 1 1 1 2 1 3 3 2 2 4 1 1 3 2 1 4 1 1 4 1 4 4 1 3 2

下面給出 $chi^2$ 檢驗，檢驗這些生成的隨機數是否服從給定的分布。

　result<-invmo(100) p<-c(0.3,0.2,0.35,0.15) Y<-table(result) chisq.test(Y,p=p)

運行的結果：

Chi-squared test for given probabilities

data: Y

X-squared = 0.70714, df = 3, p-value = 0.8715

P值＝0.8715,因此認為模擬的數據服從給定的分布。

作業2 請給出生成1，2，...，n的隨機排列，請給出演算法與R程序。

提示：先隨機地從1，2，...，n位置中任選一個位置r，然後把第r個位置上的數和第n個位置上的數進行互換；然後在從1，2，...，n－1位置中任選一個位置r，然後把第r個位置上的數和第n－1個位置上的數進行互換，依此類推。

參考答案：

perm<-function(n){ x<-1:n #x的第i個位置放i for (i in 1:n){ r<-ceiling((n-i+1)*runif(1))#從1，2，...,n-i+1中任選一個數 y<-x[n-i+1] x[n-i+1]<-x[r] x[r]<-y } x }

運行結果，一個1:100的隨機排列：

> perm(100)

[1] 22 56 17 48 57 84 80 11 61 89 73 88 1 81 45 87 29

[18] 23 27 50 68 37 60 98 40 63 18 19 39 78 66 9 35 16

[35] 59 10 51 41 21 25 13 5 96 3 34 93 6 75 76 72 91

[52] 58 79 62 44 86 47 71 77 53 28 67 99 8 82 43 90 20

[69] 4 52 30 74 33 2 24 64 69 36 31 14 70 15 55 92 65

[86] 54 94 12 95 49 85 97 46 32 38 42 7 26 100 83

作業3 用逆變換法模擬下面隨機變數 $X$ ,其分布為

$P{X=4}=0.2,P{X=6}=0.3 , P{X=8}=0.1,P{X=123}=0.4$ 。