右截尾数据线性回归EM算法

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由上面计算 */ sita1=(ex2+x5)/(ex2+x3+x4+x5); /* M- 步得到的迭代公式 */ if abs(sita1-sita)17 可以用其它的标识： da

2、ta 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, interc

3、ept); /* 创建一个值来自 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.1415926)*(-0

4、.5)*exp(-0.5*(c-beta01-beta11*v)/sigma 1)*2)/(1-probnorm(c-beta01-beta11*v)/sigma1);/*EZi =0(j) 1( j)vi ( j)H(j)ci0(j)viezv=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;

5、 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*(sum

6、vsq)-(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;si

7、gma=&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;