2015年6月27日星期六

EM alogrithm for mixed Normal distribtuion

Data are from Robert Hogg Exercise 6.6.8

library(car)
x<-c(119.0,96.0,146.2,138.6,143.4,98.2,124.5,114.1,136.2,136.4,
184.8,79.8,151.9,114.2,145.7,95.9,97.3,136.4,109.2,103.2)
densityPlot(x)
p_0<-0.6
u1_0<-105
u2_0<-130
sigma1_0<-15
sigma2_0<-25
p<-c()
pi<-c()
u1<-c()
u2<-c()
sigma1_2<-c()
sigma2_2<-c()
s<-c()
s1<-c()
diff<-c()
j<-1   #this is imporant
s<-c(p_0,u1_0,u2_0,sigma1_0^2,sigma2_0^2) #initial guess of p,u1,u2,sigma1, sigma2
p_i_0<-p_0*dnorm(x,u2_0,sigma2_0)/((1-p_0)*dnorm(x,u1_0,sigma1_0)+p_0*dnorm(x,u2_0,sigma2_0))
p[j]<-mean(p_i_0)
u1[j]<-sum((1-p_i_0)*x)/sum(1-p_i_0)
u2[j]<-sum(p_i_0*x)/sum(p_i_0)
sigma1_2[j]<-sum((1-p_i_0)*(x-u1_0)^2)/sum(1-p_i_0)
sigma2_2[j]<-sum(p_i_0*(x-u2_0)^2)/sum(p_i_0)
s1<-c(p[j],u1[j],u2[j],sigma1_2[j],sigma2_2[j])
for(i in 1:5){
diff[i]<-abs(s1[i]-s[i])
}
for(i in 1:5){     #in fact we don't need this, we can mainly focus on p
while(diff[i]>0.000001){
s<-c(p[j],u1[j],u2[j],sigma1_2[j],sigma2_2[j])
p_i<-p[j]*dnorm(x,u2[j],(sigma2_2[j])^0.5)/((1-p[j])*dnorm(x,u1[j],sigma1_2[j]^0.5)+p[j]*dnorm(x,u2[j],sigma2_2[j]^0.5))
#p_i will be updated by p,u1,u2,sigma1 and sigma2
u1[j+1]<-sum((1-p_i)*x)/sum(1-p_i)
u2[j+1]<-sum(p_i*x)/sum(p_i)
sigma1_2[j+1]<-sum((1-p_i)*(x-u1[j+1])^2)/sum(1-p_i)
sigma2_2[j+1]<-sum(p_i*(x-u2[j+1])^2)/sum(p_i)
p[j+1]<-mean(p_i)
s1<-c(p[j+1],u1[j+1],u2[j+1],sigma1_2[j+1],sigma2_2[j+1])
diff[i]<-abs(s1[i]-s[i])
j=j+1
}
}

#results:
After 301 iterations all parameters are within 0.00000 accuracy

For p:

 [1] 0.6446007 0.6508217 0.6549901 0.6582507 0.6608042 0.6629207 0.6649003
  [8] 0.6668472 0.6687443 0.6705665 0.6723035 0.6739539 0.6755194 0.6770028
 [15] 0.6784075 0.6797370 0.6809951 0.6821853 0.6833114 0.6843770 0.6853854
 [22] 0.6863398 0.6872434 0.6880990 0.6889093 0.6896769 0.6904043 0.6910936
 [29] 0.6917470 0.6923664 0.6929539 0.6935110 0.6940395 0.6945410 0.6950169
 [36] 0.6954685 0.6958972 0.6963042 0.6966907 0.6970577 0.6974063 0.6977374
 [43] 0.6980519 0.6983508 0.6986348 0.6989047 0.6991612 0.6994050 0.6996367
 [50] 0.6998571 0.7000666 0.7002657 0.7004552 0.7006353 0.7008066 0.7009696
 [57] 0.7011245 0.7012720 0.7014122 0.7015457 0.7016726 0.7017934 0.7019083
 [64] 0.7020177 0.7021217 0.7022208 0.7023150 0.7024047 0.7024901 0.7025713
 [71] 0.7026487 0.7027223 0.7027924 0.7028591 0.7029226 0.7029830 0.7030405
 [78] 0.7030953 0.7031475 0.7031971 0.7032444 0.7032894 0.7033322 0.7033730
 [85] 0.7034119 0.7034489 0.7034841 0.7035176 0.7035496 0.7035800 0.7036089
 [92] 0.7036365 0.7036627 0.7036877 0.7037115 0.7037342 0.7037558 0.7037764
 [99] 0.7037959 0.7038146 0.7038323 0.7038493 0.7038654 0.7038807 0.7038953
[106] 0.7039092 0.7039225 0.7039351 0.7039471 0.7039585 0.7039694 0.7039798
[113] 0.7039897 0.7039991 0.7040081 0.7040166 0.7040247 0.7040325 0.7040399
[120] 0.7040469 0.7040536 0.7040600 0.7040660 0.7040718 0.7040773 0.7040826
[127] 0.7040876 0.7040923 0.7040968 0.7041012 0.7041053 0.7041092 0.7041129
[134] 0.7041165 0.7041199 0.7041231 0.7041261 0.7041291 0.7041319 0.7041345
[141] 0.7041370 0.7041394 0.7041417 0.7041439 0.7041460 0.7041480 0.7041499
[148] 0.7041516 0.7041534 0.7041550 0.7041565 0.7041580 0.7041594 0.7041608
[155] 0.7041620 0.7041633 0.7041644 0.7041655 0.7041666 0.7041676 0.7041685
[162] 0.7041694 0.7041703 0.7041711 0.7041719 0.7041727 0.7041734 0.7041740
[169] 0.7041747 0.7041753 0.7041759 0.7041765 0.7041770 0.7041775 0.7041780
[176] 0.7041784 0.7041789 0.7041793 0.7041797 0.7041801 0.7041804 0.7041808
[183] 0.7041811 0.7041814 0.7041817 0.7041820 0.7041822 0.7041825 0.7041827
[190] 0.7041830 0.7041832 0.7041834 0.7041836 0.7041838 0.7041840 0.7041842
[197] 0.7041843 0.7041845 0.7041846 0.7041848 0.7041849 0.7041850 0.7041852
[204] 0.7041853 0.7041854 0.7041855 0.7041856 0.7041857 0.7041858 0.7041859
[211] 0.7041860 0.7041860 0.7041861 0.7041862 0.7041863 0.7041863 0.7041864
[218] 0.7041864 0.7041865 0.7041866 0.7041866 0.7041867 0.7041867 0.7041867
[225] 0.7041868 0.7041868 0.7041869 0.7041869 0.7041869 0.7041870 0.7041870
[232] 0.7041870 0.7041871 0.7041871 0.7041871 0.7041871 0.7041872 0.7041872
[239] 0.7041872 0.7041872 0.7041872 0.7041873 0.7041873 0.7041873 0.7041873
[246] 0.7041873 0.7041873 0.7041874 0.7041874 0.7041874 0.7041874 0.7041874
[253] 0.7041874 0.7041874 0.7041874 0.7041875 0.7041875 0.7041875 0.7041875
[260] 0.7041875 0.7041875 0.7041875 0.7041875 0.7041875 0.7041875 0.7041875
[267] 0.7041875 0.7041875 0.7041875 0.7041875 0.7041875 0.7041876 0.7041876
[274] 0.7041876 0.7041876 0.7041876 0.7041876 0.7041876 0.7041876 0.7041876
[281] 0.7041876 0.7041876 0.7041876 0.7041876 0.7041876 0.7041876 0.7041876
[288] 0.7041876 0.7041876 0.7041876 0.7041876 0.7041876 0.7041876 0.7041876
[295] 0.7041876 0.7041876 0.7041876 0.7041876 0.7041876 0.7041876 0.7041876
For u1
u1
  [1] 105.37122 104.28497 103.19437 102.19782 101.34808 100.69007 100.21667  99.88967  99.66755
 [10]  99.51631  99.41138  99.33623  99.28022  99.23661  99.20121  99.17139  99.14549  99.12247
 [19]  99.10163  99.08255  99.06490  99.04849  99.03315  99.01877  99.00527  98.99257  98.98060
 [28]  98.96932  98.95867  98.94862  98.93912  98.93015  98.92166  98.91364  98.90605  98.89886
 [37]  98.89206  98.88562  98.87952  98.87374  98.86826  98.86307  98.85815  98.85348  98.84905
 [46]  98.84485  98.84087  98.83709  98.83350  98.83009  98.82686  98.82379  98.82087  98.81810
 [55]  98.81547  98.81297  98.81059  98.80833  98.80619  98.80415  98.80221  98.80037  98.79861
 [64]  98.79695  98.79536  98.79386  98.79243  98.79106  98.78977  98.78853  98.78736  98.78624
 [73]  98.78518  98.78417  98.78321  98.78230  98.78143  98.78060  98.77981  98.77906  98.77835
 [82]  98.77767  98.77702  98.77640  98.77582  98.77526  98.77473  98.77422  98.77374  98.77329
 [91]  98.77285  98.77243  98.77204  98.77166  98.77130  98.77096  98.77064  98.77033  98.77003
[100]  98.76975  98.76949  98.76923  98.76899  98.76876  98.76854  98.76833  98.76813  98.76794
[109]  98.76776  98.76759  98.76743  98.76727  98.76712  98.76698  98.76685  98.76672  98.76660
[118]  98.76648  98.76637  98.76626  98.76616  98.76607  98.76598  98.76589  98.76581  98.76573
[127]  98.76565  98.76558  98.76551  98.76545  98.76539  98.76533  98.76527  98.76522  98.76517
[136]  98.76512  98.76507  98.76503  98.76499  98.76495  98.76491  98.76487  98.76484  98.76481
[145]  98.76478  98.76475  98.76472  98.76469  98.76467  98.76464  98.76462  98.76460  98.76457
[154]  98.76455  98.76454  98.76452  98.76450  98.76448  98.76447  98.76445  98.76444  98.76442
[163]  98.76441  98.76440  98.76439  98.76438  98.76437  98.76436  98.76435  98.76434  98.76433
[172]  98.76432  98.76431  98.76430  98.76430  98.76429  98.76428  98.76428  98.76427  98.76427
[181]  98.76426  98.76425  98.76425  98.76425  98.76424  98.76424  98.76423  98.76423  98.76423
[190]  98.76422  98.76422  98.76422  98.76421  98.76421  98.76421  98.76420  98.76420  98.76420
[199]  98.76420  98.76419  98.76419  98.76419  98.76419  98.76419  98.76419  98.76418  98.76418
[208]  98.76418  98.76418  98.76418  98.76418  98.76418  98.76417  98.76417  98.76417  98.76417
[217]  98.76417  98.76417  98.76417  98.76417  98.76417  98.76417  98.76417  98.76417  98.76416
[226]  98.76416  98.76416  98.76416  98.76416  98.76416  98.76416  98.76416  98.76416  98.76416
[235]  98.76416  98.76416  98.76416  98.76416  98.76416  98.76416  98.76416  98.76416  98.76416
[244]  98.76416  98.76416  98.76416  98.76416  98.76416  98.76416  98.76416  98.76416  98.76416
[253]  98.76415  98.76415  98.76415  98.76415  98.76415  98.76415  98.76415  98.76415  98.76415
[262]  98.76415  98.76415  98.76415  98.76415  98.76415  98.76415  98.76415  98.76415  98.76415
[271]  98.76415  98.76415  98.76415  98.76415  98.76415  98.76415  98.76415  98.76415  98.76415
[280]  98.76415  98.76415  98.76415  98.76415  98.76415  98.76415  98.76415  98.76415  98.76415
[289]  98.76415  98.76415  98.76415  98.76415  98.76415  98.76415  98.76415  98.76415  98.76415
[298]  98.76415  98.76415  98.76415  98.76415
for u2

 [1] 133.5728 133.8861 134.2721 134.6356 134.9464 135.1737 135.3096 135.3706 135.3799 135.3572
 [11] 135.3157 135.2642 135.2078 135.1499 135.0923 135.0362 134.9820 134.9302 134.8809 134.8340
 [21] 134.7895 134.7473 134.7073 134.6695 134.6337 134.5997 134.5676 134.5371 134.5083 134.4809
 [31] 134.4550 134.4305 134.4072 134.3851 134.3641 134.3442 134.3254 134.3075 134.2905 134.2743
 [41] 134.2590 134.2445 134.2307 134.2175 134.2051 134.1932 134.1820 134.1713 134.1611 134.1515
 [51] 134.1423 134.1335 134.1252 134.1173 134.1098 134.1027 134.0959 134.0895 134.0833 134.0775
 [61] 134.0719 134.0666 134.0616 134.0568 134.0523 134.0479 134.0438 134.0399 134.0361 134.0326
 [71] 134.0292 134.0260 134.0229 134.0200 134.0172 134.0146 134.0121 134.0097 134.0074 134.0052
 [81] 134.0032 134.0012 133.9993 133.9975 133.9958 133.9942 133.9927 133.9912 133.9898 133.9885
 [91] 133.9872 133.9860 133.9849 133.9838 133.9827 133.9818 133.9808 133.9799 133.9791 133.9782
[101] 133.9775 133.9767 133.9760 133.9754 133.9747 133.9741 133.9735 133.9730 133.9725 133.9720
[111] 133.9715 133.9710 133.9706 133.9702 133.9698 133.9694 133.9691 133.9687 133.9684 133.9681
[121] 133.9678 133.9675 133.9673 133.9670 133.9668 133.9665 133.9663 133.9661 133.9659 133.9657
[131] 133.9655 133.9654 133.9652 133.9651 133.9649 133.9648 133.9646 133.9645 133.9644 133.9643
[141] 133.9642 133.9640 133.9639 133.9639 133.9638 133.9637 133.9636 133.9635 133.9634 133.9634
[151] 133.9633 133.9632 133.9632 133.9631 133.9631 133.9630 133.9630 133.9629 133.9629 133.9628
[161] 133.9628 133.9627 133.9627 133.9627 133.9626 133.9626 133.9626 133.9625 133.9625 133.9625
[171] 133.9625 133.9624 133.9624 133.9624 133.9624 133.9623 133.9623 133.9623 133.9623 133.9623
[181] 133.9623 133.9622 133.9622 133.9622 133.9622 133.9622 133.9622 133.9622 133.9622 133.9621
[191] 133.9621 133.9621 133.9621 133.9621 133.9621 133.9621 133.9621 133.9621 133.9621 133.9621
[201] 133.9621 133.9621 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620
[211] 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620
[221] 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620
[231] 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620
[241] 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620 133.9620
[251] 133.9620 133.9620 133.9620 133.9620 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619
[261] 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619
[271] 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619
[281] 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619
[291] 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619 133.9619
[301] 133.9619
For sigma1 square:
 [1] 220.12489 205.10096 183.53226 160.65836 141.18629 127.45956 118.93847 113.86190 110.77411
 [10] 108.80379 107.47407 106.52342 105.80421 105.23069 104.75194 104.33723 103.96770 103.63162
 [19] 103.32155 103.03264 102.76167 102.50639 102.26517 102.03679 101.82027 101.61481 101.41971
 [28] 101.23438 101.05826 100.89085 100.73171 100.58039 100.43649 100.29964 100.16948 100.04567
 [37]  99.92790  99.81586  99.70927  99.60785  99.51136  99.41954  99.33217  99.24902  99.16990
 [46]  99.09459  99.02292  98.95471  98.88979  98.82799  98.76917  98.71318  98.65988  98.60914
 [55]  98.56084  98.51486  98.47108  98.42941  98.38973  98.35196  98.31599  98.28174  98.24914
 [64]  98.21810  98.18854  98.16039  98.13359  98.10807  98.08377  98.06063  98.03860  98.01762
 [73]  97.99764  97.97861  97.96049  97.94324  97.92681  97.91116  97.89626  97.88207  97.86855
 [82]  97.85568  97.84342  97.83175  97.82063  97.81005  97.79996  97.79036  97.78122  97.77251
 [91]  97.76421  97.75631  97.74879  97.74162  97.73479  97.72829  97.72210  97.71621  97.71059
[100]  97.70524  97.70015  97.69530  97.69067  97.68627  97.68208  97.67809  97.67429  97.67067
[109]  97.66722  97.66393  97.66080  97.65782  97.65498  97.65228  97.64970  97.64725  97.64491
[118]  97.64269  97.64057  97.63855  97.63663  97.63480  97.63305  97.63139  97.62981  97.62830
[127]  97.62687  97.62550  97.62420  97.62296  97.62178  97.62065  97.61958  97.61856  97.61759
[136]  97.61666  97.61578  97.61494  97.61414  97.61338  97.61265  97.61196  97.61130  97.61067
[145]  97.61008  97.60951  97.60897  97.60845  97.60796  97.60749  97.60704  97.60662  97.60621
[154]  97.60583  97.60546  97.60511  97.60478  97.60446  97.60416  97.60387  97.60360  97.60334
[163]  97.60309  97.60285  97.60263  97.60241  97.60221  97.60201  97.60183  97.60165  97.60148
[172]  97.60132  97.60117  97.60102  97.60088  97.60075  97.60062  97.60051  97.60039  97.60028
[181]  97.60018  97.60008  97.59999  97.59990  97.59981  97.59973  97.59965  97.59958  97.59951
[190]  97.59944  97.59938  97.59932  97.59926  97.59921  97.59915  97.59910  97.59906  97.59901
[199]  97.59897  97.59893  97.59889  97.59885  97.59881  97.59878  97.59875  97.59872  97.59869
[208]  97.59866  97.59863  97.59861  97.59859  97.59856  97.59854  97.59852  97.59850  97.59848
[217]  97.59846  97.59845  97.59843  97.59841  97.59840  97.59839  97.59837  97.59836  97.59835
[226]  97.59834  97.59832  97.59831  97.59830  97.59829  97.59829  97.59828  97.59827  97.59826
[235]  97.59825  97.59825  97.59824  97.59823  97.59823  97.59822  97.59822  97.59821  97.59821
[244]  97.59820  97.59820  97.59819  97.59819  97.59818  97.59818  97.59818  97.59817  97.59817
[253]  97.59817  97.59816  97.59816  97.59816  97.59816  97.59815  97.59815  97.59815  97.59815
[262]  97.59814  97.59814  97.59814  97.59814  97.59814  97.59814  97.59813  97.59813  97.59813
[271]  97.59813  97.59813  97.59813  97.59813  97.59813  97.59812  97.59812  97.59812  97.59812
[280]  97.59812  97.59812  97.59812  97.59812  97.59812  97.59812  97.59812  97.59812  97.59812
[289]  97.59811  97.59811  97.59811  97.59811  97.59811  97.59811  97.59811  97.59811  97.59811
[298]  97.59811  97.59811  97.59811  97.59811
For sigma2 square

 [1] 559.5878 525.6849 505.7973 488.0545 472.0863 458.8586 449.1209 442.7654 439.1009 437.3237
 [11] 436.7687 436.9631 437.5961 438.4720 439.4702 440.5184 441.5732 442.6096 443.6137 444.5780
 [21] 445.4991 446.3760 447.2092 447.9998 448.7495 449.4602 450.1336 450.7719 451.3767 451.9501
 [31] 452.4935 453.0088 453.4975 453.9609 454.4006 454.8177 455.2135 455.5891 455.9457 456.2843
 [41] 456.6057 456.9110 457.2010 457.4764 457.7380 457.9866 458.2229 458.4474 458.6608 458.8636
 [51] 459.0564 459.2397 459.4140 459.5797 459.7373 459.8872 460.0297 460.1653 460.2943 460.4170
 [61] 460.5337 460.6447 460.7503 460.8509 460.9465 461.0375 461.1241 461.2065 461.2850 461.3596
 [71] 461.4307 461.4983 461.5627 461.6240 461.6823 461.7378 461.7907 461.8410 461.8889 461.9345
 [81] 461.9779 462.0192 462.0586 462.0960 462.1317 462.1657 462.1980 462.2288 462.2581 462.2861
 [91] 462.3126 462.3380 462.3621 462.3850 462.4069 462.4277 462.4475 462.4664 462.4844 462.5015
[101] 462.5178 462.5333 462.5481 462.5622 462.5756 462.5883 462.6005 462.6121 462.6231 462.6336
[111] 462.6436 462.6532 462.6622 462.6709 462.6791 462.6869 462.6944 462.7015 462.7083 462.7147
[121] 462.7209 462.7267 462.7323 462.7376 462.7427 462.7475 462.7521 462.7564 462.7606 462.7646
[131] 462.7683 462.7719 462.7754 462.7786 462.7817 462.7847 462.7875 462.7902 462.7927 462.7952
[141] 462.7975 462.7997 462.8018 462.8038 462.8057 462.8075 462.8093 462.8109 462.8125 462.8140
[151] 462.8154 462.8167 462.8180 462.8193 462.8204 462.8216 462.8226 462.8236 462.8246 462.8255
[161] 462.8264 462.8272 462.8280 462.8288 462.8295 462.8302 462.8308 462.8315 462.8321 462.8326
[171] 462.8332 462.8337 462.8342 462.8346 462.8351 462.8355 462.8359 462.8363 462.8366 462.8370
[181] 462.8373 462.8376 462.8379 462.8382 462.8385 462.8387 462.8390 462.8392 462.8394 462.8397
[191] 462.8399 462.8401 462.8402 462.8404 462.8406 462.8407 462.8409 462.8410 462.8412 462.8413
[201] 462.8414 462.8416 462.8417 462.8418 462.8419 462.8420 462.8421 462.8422 462.8422 462.8423
[211] 462.8424 462.8425 462.8425 462.8426 462.8427 462.8427 462.8428 462.8428 462.8429 462.8429
[221] 462.8430 462.8430 462.8431 462.8431 462.8432 462.8432 462.8432 462.8433 462.8433 462.8433
[231] 462.8434 462.8434 462.8434 462.8434 462.8435 462.8435 462.8435 462.8435 462.8435 462.8436
[241] 462.8436 462.8436 462.8436 462.8436 462.8436 462.8437 462.8437 462.8437 462.8437 462.8437
[251] 462.8437 462.8437 462.8437 462.8437 462.8438 462.8438 462.8438 462.8438 462.8438 462.8438
[261] 462.8438 462.8438 462.8438 462.8438 462.8438 462.8438 462.8438 462.8438 462.8438 462.8438
[271] 462.8439 462.8439 462.8439 462.8439 462.8439 462.8439 462.8439 462.8439 462.8439 462.8439
[281] 462.8439 462.8439 462.8439 462.8439 462.8439 462.8439 462.8439 462.8439 462.8439 462.8439
[291] 462.8439 462.8439 462.8439 462.8439 462.8439 462.8439 462.8439 462.8439 462.8439 462.8439
[301] 462.8439