1、例17.1的相關SAS計算程序。EM算法計算得出丁二0.6268 :data a;sita=0.5; /*為要估計的參考sita賦初值0.5 */x3=18; /*已知條件 */x4=20;x5=34;do time=1 to 10;p=sita/(2+sita);/* p按上面公式計算*/ex2=125*p;/* x2的條件期望。x2的條件分布為二項分布，n=125, p由上面計算 */sita仁(ex2+x5)/(ex2+x3+x4+x5);/* M-步得到的迭代公式*/if abs(sita1-sita)<=0.00001 then stop;output;sita=sita1;e

2、n d;run;得到的迭代結果：obssita11 0.608252 0.624323 0.626494 0.626785 0.62682將sita=0.62682代入得到y1-y5的估計值：data a;sita=0.62682;y1= (1/2+sita/4)*197;y2=1/4*(1-sita)*197;y3=1/4*(1-sita)*197;y4=sita/4*197;format _nu meric_ 5.;put y1= y2= y3= y4=;run;y1=129 y2=18 y3=18 y4=31。估計結果和實際值很按近。不用EM算法，直接估計時會分別得到4個sita的估計值

3、：data;sita=4 *(125/197-1/2);put sita=;sita =1-4*18/197 ;put sita =;sita =1-4*20/197 ;put sita =;sita =4*34/197 ;put sita =;run;得到 sita 估計值：sita=0.538071066sita=0.6345177665sita=0.5939086294sita=0.6903553299計算用 EM 算法和直接估計得到的結果：data a;do sita=0.6268, 0.538071066, 0.6345177665, 0.5939086294, 0.69035532

4、99; y1=(1/2+sita/4)*197;y2=1/4*(1-sita)*197;y3=1/4*(1-sita)*197; y4=sita/4*197;format _numeric_ 5.;put y1= y2= y3= y4=;output;end;run;結果顯示：y1=129 y2=18 y3=18 y4=31, EM 算法的結果。y1=125 y2=23 y3=23 y4=27y1=130 y2=18 y3=18 y4=31y1=128 y2=20 y3=20 y4=29y1=132 y2=15 y3=15 y4=3417.2.2 右截尾數據簡單線性回歸計算程序創建 SAS 數

5、據集：data a1;input v1 t1;cards;170 1764170 2772170 3444170 3542170 3780170 4860170 5196190 408190 408190 1344190 1344190 1440220 408220 408220 504220 504220 504150 8064150 8064150 8064150 8064150 8064150 8064150 8064150 8064150 8064150 8064170 5448170 5448170 5448190 1680190 1680190 1680190 1680190 16

6、80220 528220 528220 528220 528220 528Jrun;按要求作數據變換，注意這里的條件 n>17 可以用其它的標識： data a;set a1;v=1000/(v1+273.2);t=log10(t1);n=_n_; /* 用于和后有參數估計的數據集合并 */vsq=v*2; /*用于求參數 beta0, beta1 和 sigma 估計 */ by_v=1; /*為了以后和 sw 合并*/if n>17 then c=t;drop v1 t1;/*直接回歸求得參數的初值，并將這些初值賦予宏變量betaO1,beta11,sigma1*/proc r

7、eg data=a outest=est n opri nt;model t=v;data est;set est;call symput('beta01', intercept); /* 創建一個值來自DATA 步的宏變量 betaOl*/call symput('beta11', v);/*創建一個值來自 DATA 步的宏變量 beta11*/call symput('sigma1', _rmse_);data w;set a ;beta01=&beta01;beta11=&beta11;sigma仁 &sigma1;

8、/*宏A求出迭代公式中的各項和，并得到迭代公式值，為下一步迭代提供值*/%macro A;data w;set w;if n >17 the n do c=t;ez=beta01+beta11*v+sigma1*(2*3.1415926)*(-0.5)*exp(-0.5*(c-beta01-beta11*v)/sigma1)*2)/(1-prob norm(c-beta01-beta11*v)/sigma1);/* EZi =罟) c(j)*/ezv=v*ez; t1=0;vt=0;hq=(2*3.1415926)*(-0.5)*exp(-0.5*(c-beta01-beta11*v)/

9、sigma1)*2)/(1-prob norm(c-beta01-beta11*v)/sigma1)*(c-beta01-beta11*v)/sigma1);*/tmu=0;en d;else do t1=t; vt=v*t; ezv=0; ez=0; hq=0;tmu=(t-beta01-beta11*v)*2;en d;proc means data=w n opri nt;var v ez ezv vt t1 hq vsq tmu sigma1;output out=sw sum=sumv sumez sumezv sumvt sumt1 sumhq sumvsq sumtmu sums

10、igma1 ;data sw;set sw;sigma1= sumsigma1/_freq_; beta0=(sumvsq)*(sumt1+sumez)-(sumv)*(sumvt+sumezv)/(40*(sumvsq)-(sumv)*2);beta1=-(sumv)*(sumt1+sumez)-40*(sumvt+sumezv)/( (40*(sumvsq)-(sumv)*2);sigma=(sumtmu/40+sigma1*2*(23+sumhq)/40)*0.5;by_v=1;keep beta0 beta1 sigma by_v;%mend A;/*將每次迭代的結果放在一個數據集re

11、sult 中，先放入直接回歸得到參數估計的初值 */data result(keep=beta0 beta1 sigma);beta0=&beta01; /* 第一個觀測為初值 */beta1=&beta11;sigma=&sigma1;options nodate nonotes nosource;/* 宏 B 是迭代程序 */%macro b;%do i=1 %to 30;%A; /* 調用宏 A*/data w;merge a sw;by by_v;rename beta0=beta01 beta1=beta11 sigma=sigma1;data result; /* 將每次迭代的結果放在一個數據集 */set result sw;%end;%mend b;%b;run;options nocenter;proc print data=result;迭代結果作圖：data result;set result;n=_n_;proc gplot data=result;symbol1 v=star i=join c=blue;symbol2 v=star i=join c=black;symbol3 v=star i=join c=red;plot beta0*n=1 beta1*n=2 sigma*n=3; run;直接回歸結