統(tǒng)計計算 課后習(xí)題答案 第五章習(xí)題答案

第五章EM算法習(xí)題5.1141

1

1，， ， 2 44 44 44(0,1125，18，20，34197次實EM算法估計參數(shù)。解：設(shè)y1=125，y2=18，y3=20，y4=34，n=197，給出似然函數(shù)211y213y4

211y213y4)

4

4

44 4n ， 先使用極大似然估計法估計參數(shù)值。LnL()y1Ln(2)y2Ln(1)y3Ln(1)y4Ln(y1y2y3y4)Ln4，dLnL() y2

y3y40，此時手算求方程的解不容易?，F(xiàn)用EMd 2

1

1 算法求θ的估計值。，將第一個結(jié)果分為兩部分，這兩部分的概率分別為11，出現(xiàn)的次數(shù)分，4 z，yz

1，出現(xiàn)的次1 1 1數(shù)分別為z2，y2z2，此時似然函數(shù)111111y2z213z2y4)

，44yz11234(yyz11234z2y4121 44 4 44 4 4先使用極大似然估計法估計參數(shù)值。LnL()(z2y4)Ln(y2z1)Ln(1)(y1y2y3y4)Ln4zb(y,1)，E(z)1

y，z

b(y,)，E(z

)y

，這里的θ是1 12

1 21

2 31

2 1 31(i) (i)指的本輪的輸入，因此E(z1)2(i)y1，E(z2)1(i)y3。E步：E(LnL())(E(z2)y4)Ln(y2E(z1))Ln(1)(y1y2y3y4)Ln4(i)( (i)y3y4(y2

1(i)，(i)y1)Ln(1)(y1y2y3y4)Ln4，1 2M步： dE(LnL())

(i)

1 1 1(i) ，d 1(i)y3y41(y22(i)y1)0

1(i)

(i)1(i)y3y4 。(i)2(i)y1y21(i)y3y4EMθ的估計值。1EM算法：importnumpyasnp#使用EM算法求解x=[125,18,20,34]theta0=0.5n=np.sum(x)foriinrange(20):s=x[3]+(theta0/(theta0+1))*x[2]t=((1-theta0)/(2-theta0))*x[0]+x[1]theta0=s/(s+t)print('使用EM算法得到θ估計值為',theta0)輸出結(jié)果：EM算法得到θ估計值為0.4053156146179402EM算法得到θ估計值為0.3809850048870182EM算法得到θ估計值為0.37524953892338536EM算法得到θ估計值為0.37392273631605066EM算法得到θ估計值為0.3736171055046297EM算法得到θ估計值為0.37354677154393484EM算法得到θ估計值為0.3735305894083948EM算法得到θ估計值為0.373526866483774EM算法得到θ估計值為0.3735260099834931EM算法得到θ估計值為0.3735258129365993EM算法得到θ估計值為0.37352576760391815EM算法得到θ估計值為0.3735257571746663EM算法得到θ估計值為0.3735257547753092EM算法得到θ估計值為0.37352575422331236EM算法得到θ估計值為0.3735257540963198EM算法得到θ估計值為0.37352575406710387EM算法得到θ估計值為0.3735257540603824EM算法得到θ估計值為0.37352575405883603EM算法得到θ估計值為0.3735257540584803EM算法得到θ估計值為0.37352575405839844結(jié)果分析：可以看出第19次和第20次的θ估計值為幾乎相同，可以認(rèn)為估計值為0.3735257代碼２使用極大似然估計法：fromscipy.optimizeimportrootdata=[125,18,20,34]deff(x):y=-data[0]/(2-x)-data[1]/(1-x)+data[2]/(1+x)+data[3]/xreturnya1=root(f,[0.3])print(a1)輸出結(jié)果為fjac:array([[-1.]])fun:array([0.])message:'Thesolutionconverged.'nfev:8qtf:array([-8.6407681e-10])r:array([347.40621023])status:1success:Truex:array([0.37352575])從而看出方程的解為0.37352575，即求得的參數(shù)θ值。2、有兩個袋子，每個袋子中都有紅白兩種顏色的球，現(xiàn)從中選取一個袋子，從10，然后放回，再任意抽取10為w(0)0.5,p(0)q(0)0.5。ZwXX1表示紅色，0X1X2,Xnw,p,q。i i iX的密度函數(shù)為wpi1pi1)qi1qiX和X1分別代入密度函數(shù)，得到(wp1)q)i(1p)11qi。i i innL(,p,q)(wpi1pi1)qi1qi)i1n(wp1)q)i(1p)11qii1nn記twp1)qL(,p,q)ti1ti，取對數(shù)得：i1n n ni iLnL(,p,q)Lnti1tixLt 1xLn1t)，i1

i1n n

i1dLnL(t)

(1) 0，

i1 i1 0，得到X，wp(1w)qX。

dt t

1t使用EM算法，由貝葉斯公式求取出的球的顏色已知的條件下該球從第一個wpi1pi

) (i) (i)袋子中取出的概率uwpx(1p)1x

(1w)qx(1q)1x

，已知w,p

q ，則u(i)

i i i i，此時w(i)(p(i))i1p(i)i，此時w(i)(p(i))i1p(i)i1w(i)(q(i))i1q(i)in n nu(i)

u(i)x

(1u(i))xi innw(i1)i1 ,p(i1)i1 ，q(i1)i1 w(0p(0)q(0時，nnn i1

u(i)

i1

(1u(i))可以使用EM算法求出參數(shù)w,p,q。給定初值w(0)0.5,p(0)0.3,q(0)0.5，結(jié)果為1，1，0，1，0，0，1，0，1 wpi1piXi=1uiwpx(1p)1x

(1w)qx(1q)1x

，得u10.375，同理得i i i iu20.375，u30.5833，……，0.375，（1）[0.37,0.37,0.583,0.37,0.583,0.583,0.37,0.583,0.37,0.37]，n n nu(i)

u(i)x

(1u(i))xi innw(1)i1 0.45832,p(1)i1 =0.4909，q(1)i1 0.6923，nnn這就是第一輪結(jié)果。輸出結(jié)果：0.374999999999999940.374999999999999940.58333333333333340.374999999999999940.58333333333333340.58333333333333340.374999999999999940.58333333333333340.37499999999999994

i1

u(i)

i1

(1u(i))0.374999999999999940.4583333333333333[0.49090909090909085,0.6923076923076923]0.4583333333333332[0.49090909090909096,0.6923076923076923]0.4583333333333333[0.49090909090909085,0.6923076923076923]0.4583333333333332[0.49090909090909096,0.6923076923076923]0.4583333333333333[0.49090909090909085,0.6923076923076923]結(jié)果分析：從第一個袋子中取球的概率w為0.4583，從第一個袋子中取到的紅色的球的概率p為0.4909，從第二個袋子中取到的紅色的球的概率q為0.6923。3w,p,q0.4，0.6，ABC0.40641711229946526[0.5368421052631579,0.6432432432432432]0.40641711229946526[0.5368421052631579,0.6432432432432431]0.40641711229946526[0.536842105263158,0.6432432432432431]0.40641711229946526[0.5368421052631581,0.6432432432432431]0.40641711229946526[0.5368421052631581,0.6432432432432431]結(jié)果分析：硬幣A正面向上的概率為0.4064，硬幣B正面向上的概率為0.5368，硬幣C正面向上的概率為0.6432。0.5，0.5，0.5，則輸出結(jié)果為：0.5[0.6,0.6]EM算法對初值十分敏感。習(xí)題5.21、由三個一維正態(tài)分布的線性組合組成的混合正態(tài)分布1X12X23X3，為權(quán)重，XN(,2)，i1,2,,3，使用觀測的數(shù)據(jù)估計參數(shù),,。.i i i i i i i取初值為0.1N(2,1)0.6N(6,5)0.3N(8,10)，從三個正態(tài)分布N(6(8,10)1：6：3的數(shù)據(jù)，由這些數(shù)據(jù)估計參數(shù)。算法：k k k已知(j),(j),(j)，則k k kp(x,z,) (j)N(x,,)（1）E步。概率wik

(z)

ip(x,z,)

k i k k 。m(j)mizk1

kN(xi,k,k)k k k（2）M步。求出參數(shù)(j1),(j1),(j1)k k knw nw (x)(j2iki ki1 nikw(ji1w x(jni1 ，

(j1)ik i(ji1 ，(j) 。wnk k wn(j1)iki1k k k按照算法，當(dāng)給定初值(0),(0),(0)，可以使用EM算法求出參數(shù),,。k k k輸出最后五個結(jié)果：μ=[1.914881915496394,5.851337772747595,8.080449143374715]σ=[0.8578868776616748,4.897071949343551,9.701635337069677]權(quán)重=[0.08914767443684753,0.5910554280516397,0.31979689751151275]μ=[1.9148819102044665,5.851337748629688,8.08044903272749]σ=[0.8578868792303709,4.897071857491882,9.701635172058804]權(quán)重=[0.08914767460447753,0.5910554051075384,0.3197969202879842]μ=[1.914881905014679,5.851337724977286,8.080448924215876]σ=[0.8578868807687922,4.897071767413045,9.701635010232819]權(quán)重=[0.0891476747688724,0.5910553826062808,0.3197969426248468]μ=[1.9148818999250605,5.8513377017814046,8.080448817798652]σ=[0.8578868822775229,4.89707167907282,9.70163485153025]權(quán)重=[0.08914767493009461,0.5910553605393197,0.31979696453058565]μ=[1.914881894933677,5.851337679033228,8.080448713435395]σ=[0.857886883757136,4.8970715924376504,9.701634695890807]權(quán)重=[0.08914767508820534,0.5910553388982729,0.31979698601352174]結(jié)果分析：1=1.91σ1=0.62=5.85σ2=4.90=8.08，σ3=9.70，權(quán)重分別為0.089，0.59，0.32，與真實參數(shù)值很接近。2、由三個二維正態(tài)分布的線性組合組成的混合正態(tài)分布1X12X23X3，為權(quán)重XN(,22,3，使用觀測的數(shù)據(jù)估計參數(shù),。.i i i i i i i取初值為1N()1N()1N()，從三個正態(tài)分布3 3 3N(,0N(,N(,0)111據(jù)估計參數(shù)。協(xié)方差矩陣均為1 0，均值矩陣分別為060。0 11 0 9 0 11 0 9 輸出最后四個結(jié)果：

μ=[array([0.08674084,-1.01081046]),array([5.97164231,0.02066913]),array([-1.72512975e-03,9.03690716e+00])]σ=[array([0.975921,1.01688392]),array([0.98415799,0.99543767]),array([1.02560005,0.99113734])]權(quán)重=[0.3339991174706056,0.3326675491949382,0.3333333333344561]μ=[array([0.08674084,-1.01081046]),array([5.97164231,0.02066913]),array([-1.72512975e-03,9.03690716e+00])]σ=[array([0.975921,1.01688392]),array([0.98415799,0.99543767]),array([1.02560005,0.99113734])]權(quán)重=[0.3339991174703147,0.33266754919522923,0.3333333333344561]μ=[array([0.08674084,-1.01081046]),array([5.97164231,0.02066913]),array([-1.72512975e-03,9.03690716e+00])]σ=[array([0.975921,1.01688392]),array([0.98415799,0.99543767]),array([1.02560005,0.99113734])]權(quán)重=[0.3339991174702905,0.3326675491952534,0.3333333333344561]μ=[array([0.08674084,-1.01081046]),array([5.97164231,0.02066913]),array([-1.72512975e-03,9.03690716e+00])]σ=[array([0.975921,1.01688392]),array([0.98415799,0.99543767]),array([1.02560005,0.99113734])]權(quán)重=[0.3339991174702885,0.33266754919525543,0.3333333333344561]結(jié)果分析：由最后五輪的結(jié)果可以看出，估計的參數(shù)μ1=[-0.0255,-0.9627]，σ1=[1.0111,1.0089]，μ2=[5.901,0.0099]，σ2=[0.9749,1.0064]，μ3=[-0.00379,9.029]，σ3=[1.0319,1.0031]0.3329，0.3338，0.3333。3、由三個獨立的gamma分布的線性組合組成的分布X1X2X3，其中X 1 1

)，i1,2,,3，假設(shè)真實參數(shù)值為0.6，0.25，

0.15，iGa(,i22i

1 2 3Ga(,從這三個伽瑪分布 1 1Ga(,

1 1Ga(,

1 1)，Ga(,

)中各取等比例的數(shù)據(jù)，221

222

223由這些數(shù)據(jù)估計參數(shù)i23EM算法估計參數(shù)。輸出結(jié)果：運行次數(shù)為500次，則最后6次輸出結(jié)果為：lambda=[0.14364879643841305,0.2893515589824163,0.5669996445791706]lambda=[0.5669996445791635,0.28935155898242576,0.14364879643841072]lambda=[0.14364879643841288,0.28935155898241693,0.5669996445791702]lambda=[0.5669996445791639,0.2893515589824252,0.14364879643841086]lambda=[0.14364879643841277,0.28935155898241743,0.5669996445791698]lambda=[0.5669996445791643,0.2893515589824247,0.14364879643841094]結(jié)果分析：λ值分別為0.5670，0.2893，0.1436。475，試估計兩個分量的高斯混合模型的5個參數(shù)。分析：通過數(shù)據(jù)可以看出，數(shù)據(jù)可以分為兩組，一組均值在-50附近，另一30mu=[30,-50]0.5100。輸出結(jié)果：迭代次數(shù) 21μ=[32.98488732-57.51107686]σ=[20.723376749.49999356]權(quán)重[0.866827720.13317228][32.98488732-57.51107686]迭代次數(shù) 22μ=[32.98488732-57.51107686]σ=[20.723376729.49999356]權(quán)重[0.866827720.13317228][32.98488732-57.51107686]迭代次數(shù) 23μ=[32.98488732-57.51107686]σ=[20.723376729.49999356]權(quán)重[0.866827720.13317228][32.98488732-57.51107686]迭代次數(shù) 24μ=[32.98488732-57.51107686]σ=[20.723376729.49999356]權(quán)重[0.866827720.13317228][32.98488732-57.51107686]迭代次數(shù) 25μ=[32.98488732-57.51107686]σ=[20.723376729.49999356]權(quán)重[0.866827720.13317228]輸出最后五次得結(jié)果，可以看出兩個正態(tài)分布分別為N（32.98,20722N(57.51,9.52)0.8668，0.13325、已知服從混合正態(tài)分布的數(shù)據(jù)如下：2.809,3.938,18.265,3.765,3.003,2.602,3.774,4.290,4.454,3.616,18.197,3.039,2.477,18.167,2.687,18.052,18.501,18.457，試求出兩個正態(tài)分布的參數(shù)以及權(quán)重。2.518mu=[2.5,18]0.5，標(biāo)準(zhǔn)差可以取的兩者相同，但是要大一些。輸出結(jié)果：迭代次數(shù) 16 μ= [3.370035946757747, 18.19709035808891] σ=[0.6537780076019535, 1.505948027842325] 權(quán)重= [0.6649055223372459,0.33509447766275424]迭代次數(shù) 17 μ= [3.371166666666667, 18.27316666666667] σ=[0.6542533281660868, 0.17629637058763784] 權(quán)重= [0.6666666666666666,0.3333333333333333]迭代次數(shù) 18 μ= [3.371166666666667, 18.27316666666667] σ=[0.6542523510763174, 0.15903712005552567] 權(quán)重= [0.6666666666666666,0.3333333333333333]迭代次數(shù) 19 μ= [3.371166666666667, 18.27316666666667] σ=[0.6542523510763174, 0.15903712005552567] 權(quán)重= [0.6666666666666666,0.3333333333333333]迭代次數(shù) 20 μ= [3.371166666666667, 18.27316666666667] σ=[0.6542523510763174, 0.15903712005552567] 權(quán)重= [0.6666666666666666,0.3333333333333333]通過最后五次輸出結(jié)果，可知：（3.37,0.652，(18.27,0.162)0.6667，0.3333。6（EM估計XN(,1nXXXmXmXmXmn丟失，也可以認(rèn)為是截尾數(shù)據(jù)。目前有18個數(shù)據(jù)：2.809,3.938,4.265,3.765,3.003,2.602,3.774,4.290,4.454,3.616,4.197,3.039,2.477,4.167,2.687,4.052,3.501,3.457，截斷點為4.5，7個數(shù)據(jù)丟失。184.525個數(shù)據(jù)的均值3.8237EM算法估計總體的均值。（2）EME步求期望的過程使用蒙特卡洛方法得到的平均值近似EMEM算法估計總體的均值。注：XN(,2)，如果限制X的取值在區(qū)間[a,b]上，此時X的密度函數(shù)為f(x;a,b,,2)函數(shù)。

f(x)F(b)F(a)

,x[a,b]，其中f(x)，F(xiàn)(x)為X的密度函數(shù)和分布(x)截尾正態(tài)分布的密度函數(shù)為f(x;a,b,,2) ,x[a,b](b)(a) (b)(a)期望為E(X) (b)(a) 解：缺失數(shù)據(jù)的密度函數(shù)為p(zi)

(zi)21 1 e 21 1

,zia，因此似然函1(a) i12m1 (yi)2nm 1 1 (zi)2i12數(shù)為L()2i2

e 1(a) e

2,zia，對數(shù)似然函數(shù)為：mLnL() (

1 (yi)

nm)

(zi)Ln

2 Ln a2i12

i1

1( ) 22 (zi)2E(LnL())

m 1(Ln

(yi)

nm )

2 2 ELn a2i1

i1

1(

)zimm

z(i)ii(Ln 1

(y)2 2i ) 2

( )

( )) adz) ai1

2 a

1(

a

(i) i1( )zi

z(i)iidE(LnL()m

(y)mdLnii

( )

( )) adz) add id

ai1

1(

a

(i) i1( )z(i)m nm

(i )(y)

(z) dzii1

a ii1

a(i) i1( ) m nm yim (Ezi)i1m

i1yim(nm)(Ez)0i1myi(nm)(Ez1)i1 ，E(Z)

(a)n 1 1(a)將a=4.5代入，給出算法如下：（1）初值為(0)=3.8237，計算E(Z(i))(i1)m

(4.5(i1))1(a(i1))yi

(nm)E(z(i))（2）再計算?(i)i1 ，迭代直至收斂即可。n輸出結(jié)果：迭代次數(shù)1μ=3.9906502194746533.82372迭代次數(shù)2μ=4.0025208966110573.82372迭代