统计计算 课后习题答案 第五章习题答案

第五章EM算法习题5.1141

1

1，， ， 2 44 44 44(0,1125，18，20，34197次实EM算法估计参数。解：设y1=125，y2=18，y3=20，y4=34，n=197，给出似然函数211y213y4

211y213y4)

4

4

44 4n ， 先使用极大似然估计法估计参数值。LnL()y1Ln(2)y2Ln(1)y3Ln(1)y4Ln(y1y2y3y4)Ln4，dLnL() y2

y3y40，此时手算求方程的解不容易。现用EMd 2

1

1 算法求θ的估计值。，将第一个结果分为两部分，这两部分的概率分别为11，出现的次数分，4 z，yz

1，出现的次1 1 1数分别为z2，y2z2，此时似然函数111111y2z213z2y4)

，44yz11234(yyz11234z2y4121 44 4 44 4 4先使用极大似然估计法估计参数值。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)输出结果：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结果分析：可以看出第19次和第20次的θ估计值为几乎相同，可以认为估计值为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)输出结果为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，即求得的参数θ值。2、有两个袋子，每个袋子中都有红白两种颜色的球，现从中选取一个袋子，从10，然后放回，再任意抽取10为w(0)0.5,p(0)q(0)0.5。ZwXX1表示红色，0X1X2,Xnw,p,q。i i iX的密度函数为wpi1pi1)qi1qiX和X1分别代入密度函数，得到(wp1)q)i(1p)11qi。i i innL(,p,q)(wpi1pi1)qi1qi)i1n(wp1)q)i(1p)11qii1nn记twp1)qL(,p,q)ti1ti，取对数得：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算法求出参数w,p,q。给定初值w(0)0.5,p(0)0.3,q(0)0.5，结果为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这就是第一轮结果。输出结果：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]结果分析：从第一个袋子中取球的概率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]结果分析：硬币A正面向上的概率为0.4064，硬币B正面向上的概率为0.5368，硬币C正面向上的概率为0.6432。0.5，0.5，0.5，则输出结果为：0.5[0.6,0.6]EM算法对初值十分敏感。习题5.21、由三个一维正态分布的线性组合组成的混合正态分布1X12X23X3，为权重，XN(,2)，i1,2,,3，使用观测的数据估计参数,,。.i i i i i i i取初值为0.1N(2,1)0.6N(6,5)0.3N(8,10)，从三个正态分布N(6(8,10)1：6：3的数据，由这些数据估计参数。算法：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步。求出参数(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按照算法，当给定初值(0),(0),(0)，可以使用EM算法求出参数,,。k k k输出最后五个结果：μ=[1.914881915496394,5.851337772747595,8.080449143374715]σ=[0.8578868776616748,4.897071949343551,9.701635337069677]权重=[0.08914767443684753,0.5910554280516397,0.31979689751151275]μ=[1.9148819102044665,5.851337748629688,8.08044903272749]σ=[0.8578868792303709,4.897071857491882,9.701635172058804]权重=[0.08914767460447753,0.5910554051075384,0.3197969202879842]μ=[1.914881905014679,5.851337724977286,8.080448924215876]σ=[0.8578868807687922,4.897071767413045,9.701635010232819]权重=[0.0891476747688724,0.5910553826062808,0.3197969426248468]μ=[1.9148818999250605,5.8513377017814046,8.080448817798652]σ=[0.8578868822775229,4.89707167907282,9.70163485153025]权重=[0.08914767493009461,0.5910553605393197,0.31979696453058565]μ=[1.914881894933677,5.851337679033228,8.080448713435395]σ=[0.857886883757136,4.8970715924376504,9.701634695890807]权重=[0.08914767508820534,0.5910553388982729,0.31979698601352174]结果分析：1=1.91σ1=0.62=5.85σ2=4.90=8.08，σ3=9.70，权重分别为0.089，0.59，0.32，与真实参数值很接近。2、由三个二维正态分布的线性组合组成的混合正态分布1X12X23X3，为权重XN(,22,3，使用观测的数据估计参数,。.i i i i i i i取初值为1N()1N()1N()，从三个正态分布3 3 3N(,0N(,N(,0)111据估计参数。协方差矩阵均为1 0，均值矩阵分别为060。0 11 0 9 0 11 0 9 输出最后四个结果：

μ=[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])]权重=[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])]权重=[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])]权重=[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])]权重=[0.3339991174702885,0.33266754919525543,0.3333333333344561]结果分析：由最后五轮的结果可以看出，估计的参数μ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，假设真实参数值为0.6，0.25，

0.15，iGa(,i22i

1 2 3Ga(,从这三个伽玛分布 1 1Ga(,

1 1Ga(,

1 1)，Ga(,

)中各取等比例的数据，221

222

223由这些数据估计参数i23EM算法估计参数。输出结果：运行次数为500次，则最后6次输出结果为：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]结果分析：λ值分别为0.5670，0.2893，0.1436。475，试估计两个分量的高斯混合模型的5个参数。分析：通过数据可以看出，数据可以分为两组，一组均值在-50附近，另一30mu=[30,-50]0.5100。输出结果：迭代次数 21μ=[32.98488732-57.51107686]σ=[20.723376749.49999356]权重[0.866827720.13317228][32.98488732-57.51107686]迭代次数 22μ=[32.98488732-57.51107686]σ=[20.723376729.49999356]权重[0.866827720.13317228][32.98488732-57.51107686]迭代次数 23μ=[32.98488732-57.51107686]σ=[20.723376729.49999356]权重[0.866827720.13317228][32.98488732-57.51107686]迭代次数 24μ=[32.98488732-57.51107686]σ=[20.723376729.49999356]权重[0.866827720.13317228][32.98488732-57.51107686]迭代次数 25μ=[32.98488732-57.51107686]σ=[20.723376729.49999356]权重[0.866827720.13317228]输出最后五次得结果，可以看出两个正态分布分别为N（32.98,20722N(57.51,9.52)0.8668，0.13325、已知服从混合正态分布的数据如下：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，试求出两个正态分布的参数以及权重。2.518mu=[2.5,18]0.5，标准差可以取的两者相同，但是要大一些。输出结果：迭代次数 16 μ= [3.370035946757747, 18.19709035808891] σ=[0.6537780076019535, 1.505948027842325] 权重= [0.6649055223372459,0.33509447766275424]迭代次数 17 μ= [3.371166666666667, 18.27316666666667] σ=[0.6542533281660868, 0.17629637058763784] 权重= [0.6666666666666666,0.3333333333333333]迭代次数 18 μ= [3.371166666666667, 18.27316666666667] σ=[0.6542523510763174, 0.15903712005552567] 权重= [0.6666666666666666,0.3333333333333333]迭代次数 19 μ= [3.371166666666667, 18.27316666666667] σ=[0.6542523510763174, 0.15903712005552567] 权重= [0.6666666666666666,0.3333333333333333]迭代次数 20 μ= [3.371166666666667, 18.27316666666667] σ=[0.6542523510763174, 0.15903712005552567] 权重= [0.6666666666666666,0.3333333333333333]通过最后五次输出结果，可知：（3.37,0.652，(18.27,0.162)0.6667，0.3333。6（EM估计XN(,1nXXXmXmXmXmn丢失，也可以认为是截尾数据。目前有18个数据：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个数据丢失。184.525个数据的均值3.8237EM算法估计总体的均值。（2）EME步求期望的过程使用蒙特卡洛方法得到的平均值近似EMEM算法估计总体的均值。注：XN(,2)，如果限制X的取值在区间[a,b]上，此时X的密度函数为f(x;a,b,,2)函数。

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

,x[a,b]，其中f(x)，F(x)为X的密度函数和分布(x)截尾正态分布的密度函数为f(x;a,b,,2) ,x[a,b](b)(a) (b)(a)期望为E(X) (b)(a) 解：缺失数据的密度函数为p(zi)

(zi)21 1 e 21 1

,zia，因此似然函1(a) i12m1 (yi)2nm 1 1 (zi)2i12数为L()2i2

e 1(a) e

2,zia，对数似然函数为：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输出结果：迭代次数1μ=3.9906502194746533.82372迭代次数2μ=4.0025208966110573.82372迭代