右截尾数据线性回归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)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 n17 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 n17 then do c=t;ez=beta01+beta11*v+sigma1*(2*3.)*(-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.)*(-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;