R語言與點(diǎn)估計(jì)學(xué)習(xí)筆記（EM算法與Bootstrap法）

一、EM算法

EM算法是一種在觀測(cè)到數(shù)據(jù)后，用迭代法估計(jì)未知參數(shù)的方法?？梢宰C明EM算法得到的序列是穩(wěn)定單調(diào)遞增的。這種算法對(duì)于截尾數(shù)據(jù)或參數(shù)中有一些我們不感興趣的參數(shù)時(shí)特別有效。

EM算法的步驟為：

E-step（求期望）：在給定y及theta=theta（i）的條件下，求關(guān)于完全數(shù)據(jù)對(duì)數(shù)似然關(guān)于潛在變量z的期望

M-step（求極值）：求上述期望關(guān)于theta的最大值theta（i+1）

重復(fù)以上兩步，直至收斂即可得到theta的MLE。

從上面的算法我們可以看到對(duì)于一個(gè)參數(shù)的情況，EM僅僅只是求解MLE的一個(gè)迭代算法。M-step做得就是optimize函數(shù)做得事情。對(duì)于EM算法，我們也沒有現(xiàn)成的求解函數(shù)（這個(gè)是自然的），我們一樣可以通過人機(jī)交互的辦法處理。

先舉一個(gè)一元的例子：

設(shè)一次實(shí)驗(yàn)可能有4個(gè)結(jié)果，發(fā)生概率分別為0.5+theta/4, 0.25-theta/4 ,0.25-theta/4 ,theta/4.其中theta在0,1之間?，F(xiàn)進(jìn)行了197次實(shí)驗(yàn)，結(jié)果發(fā)生的次數(shù)分別為：125,18,20,34，求theta的MLE。

計(jì)算出theta（i+1）=（195theta（i）+68）/（197theta（i）+144）

為什么是這個(gè)結(jié)果，請(qǐng)翻閱王兆軍《數(shù)理統(tǒng)計(jì)講義》p43-p44

我們用簡(jiǎn)單的循環(huán)就可以解決這個(gè)問題，程序及結(jié)果如下：

>fun<-function(error=1e-7){

+theta<-0.5

+k<-1

+while(T){

+k<-k+1

+theta[k]<-(159*theta[k-1]+68)/(197*theta[k-1]+144)

+if(abs(theta[k]-theta[k-1])<error) break

+}

+list(theta<-theta[k],iter<-k)

+}

>fun()

[[1]]

[1]0.6268215

[[2]]

[1]9

我們?cè)倏匆粋€(gè)二元的簡(jiǎn)單例子：

幼兒園里老師給a，b，c，d四個(gè)小朋友發(fā)糖吃，但老師有點(diǎn)偏心，不同小朋友得到糖的概率不同，p(a)=0.5,p(b)=miu, p(c)=2*miu, p(d)=0.5-3*miu 如果確定了參數(shù)miu，概率分布就知道了。我們可以通過觀察樣本數(shù)據(jù)來推測(cè)參數(shù)知道c和d二人得到的糖果數(shù)，也知道a與b二人的糖果數(shù)之和為h，如何來估計(jì)出參數(shù)miu呢？前面我們知道了，如果觀察到a，b，c，d就可以用ML估計(jì)出miu。反之，如果miu已知，根據(jù)概率期望 a/b=0.5/miu，又有a+b=h。由兩個(gè)式子可得到 a=0.5*h/(0.5+miu)和b=miu*h/(0.5+miu)。

># 已知條件

>

>h = 20

>c = 10

>d = 10

>

># 隨機(jī)初始兩個(gè)未知量

>miu = runif(1,0,1/6)

>b = round(runif(1,1,20))

>

>iter = 1

>nonstop=TRUE

>while (nonstop) {

+     # E步驟，根據(jù)假設(shè)的miu來算b

+     b = c(b,miu[iter]*h/(0.5+miu[iter]))

+     print(b)

+     # M步驟，根據(jù)上面算出的b再來計(jì)算miu

+     miu = c(miu,(b[iter+1] +c)/(6*(b[iter+1]+c+d))) 

+     print(miu)

+     # 記錄循環(huán)次數(shù)

+     iter = iter + 1

+     # 如果前后兩次的計(jì)算結(jié)果差距很小則退出

+     nonstop =((miu[iter]-miu[iter-1])>10^(-10))

+}

[1]3.000000 4.450531

[1]0.14310878 0.09850182

>print(cbind(miu,b))

miu        b

[1,]0.14310878 3.000000

[2,]0.09850182 4.450531

關(guān)于EM算法，及后續(xù)的發(fā)展GME的理論你可以在多數(shù)數(shù)理統(tǒng)計(jì)書上找到相關(guān)結(jié)論，也可以用類似辦法編寫函數(shù)處理它。

二、       自助法（bootstrap）

Bootstrap法是以原始數(shù)據(jù)為基礎(chǔ)的模擬抽樣統(tǒng)計(jì)推斷法,可用于研究一組數(shù)據(jù)的某統(tǒng)計(jì)量的分布特征,特別適用于那些難以用常規(guī)方法導(dǎo)出對(duì)參數(shù)的區(qū)間估計(jì)、假設(shè)檢驗(yàn)等問題?！癇ootstrap”的基本思想是:在原始數(shù)據(jù)的圍內(nèi)作有放回的再抽樣,樣本含量仍為n,原始數(shù)據(jù)中每個(gè)觀察單位每次被抽到的概率相等,為1，…,n,所得樣本稱為bootstrap樣本。于是可得到參數(shù)Η的一個(gè)估計(jì)值Η(b),這樣重復(fù)若干次,記為B。設(shè)B=1000,就得到該參數(shù)的1000個(gè)估計(jì)值,則參數(shù)Η的標(biāo)準(zhǔn)誤的bootstrap估計(jì)。簡(jiǎn)而言之就是：既然樣本是抽出來的，那我何不從樣本中再抽樣。

我們知道，如果分布函數(shù)F是已知的。在理論上就能夠計(jì)算出參數(shù)的估計(jì)量的均方誤差．若分布函數(shù)f未知，由格里文科-康特利定理知，當(dāng)M充分大時(shí)，經(jīng)驗(yàn)分布函數(shù)以概率1一致收斂到F。

我們舉一例：利用bootstrap法估計(jì)標(biāo)準(zhǔn)正態(tài)分布隨機(jī)變量的期望theta=E（X）

>gauss<-rnorm(100,2,6)

>boot<-0

>for(i in 1:1000){

+boot[i]=mean(sample(gauss,replace=T))

+}

>summary(boot)

Min. 1st Qu. Median    Mean 3rd Qu.    Max.

0.3345 1.9540  2.3350  2.3230 2.7020  4.2330

>summary(gauss)

Min. 1st Qu. Median    Mean 3rd Qu.    Max.

-13.380  -2.238  2.570   2.296   6.861 16.230

>sd(boot)

[1]0.599087

>sd(gauss)/sqrt(100)

[1]0.5906275

結(jié)果分析：

需要指出的是bootstrap法不是為了提高估計(jì)量的精度．而是一般用來對(duì)估計(jì)量的方差進(jìn)行估計(jì)。