右截尾数据线性回归EM算法
例17.1的相关SAS计算程序。EM算法计算得出：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由上面计算 */sita1=(ex2+x5)/(ex2+x3+x4+x5); /* M-步得到的迭代公式 */if abs(sita1-sita)<=0.00001 then stop; output;sita=sita1;end;run;得到的迭代结果：obs sita11 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 _numeric_ 5.;put y1= y2= y3= y4= ;run;y1=129 y2=18 y3=18 y4=31。估计结果和实际值很按近。不用EM算法，直接估计时会分别得到4个sita的估计值：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.6903553299;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数据集：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 1680220 528220 528220 528220 528220 528;run;按要求作数据变换，注意这里的条件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;/*直接回归求得参数的初值，并将这些初值赋予宏变量beta01,beta11,sigma1*/proc reg data=a outest=est noprint;model t=v;data est; set est; call symput('beta01', intercept); /*创建一个值来自DATA步的宏变量beta01*/call symput('beta11', v); /*创建一个值来自DATA步的宏变量beta11*/call symput('sigma1', _rmse_);data w;set a ;beta01=&beta01;beta11=&beta11;sigma1=&sigma1;/*宏A求出迭代公式中的各项和，并得到迭代公式值，为下一步迭代提供值*/%macro A;data w;set w;if n>17 then do c=t;ez=beta01+beta11*v+sigma1*(2*3.1415926)*(-0.5)*exp(-0.5*(c-beta01-beta11*v)/sigma1)*2)/(1-probnorm(c-beta01-beta11*v)/sigma1);/*=*/ezv=v*ez; t1=0; vt=0;hq=(2*3.1415926)*(-0.5)*exp(-0.5*(c-beta01-beta11*v)/sigma1)*2)/(1-probnorm(c-beta01-beta11*v)/sigma1)*(c-beta01-beta11*v)/sigma1);/*hq=*/tmu=0;end;else do t1=t; vt=v*t; ezv=0; ez=0; hq=0;tmu=(t-beta01-beta11*v)*2;end;proc means data=w noprint;var v ez ezv vt t1 hq vsq tmu sigma1;output out=sw sum=sumv sumez sumezv sumvt sumt1 sumhq sumvsq sumtmu sumsigma1 ;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;/*将每次迭代的结果放在一个数据集result中，先放入直接回归得到参数估计的初值*/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;直接回归结果：data a;set a1;v=1000/(v1+273.2);t=log10(t1);proc reg data=a;model t=v/dw;run;直接回归的参数估计值： -4.93051 3.74704 两种估计方法得到的误差：data a;set a;r1=t-4.93051+3.74704*v;r2=t-6.019+4.311*v;run;
