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

精选优质文档-----倾情为你奉上精选优质文档-----倾情为你奉上专心---专注---专业专心---专注---专业精选优质文档-----倾情为你奉上专心---专注---专业例17.1的相关SAS计算程序。EM算法计算得出：dataa;sita=0.5;/*为要估计的参考sita赋初值0.5*/x3=18;/*已知条件*/x4=20;x5=34;dotime=1to10;p=sita/(2+sita);/*p按上面公式计算*/ex2=125*p;/*x2的条件期望。x2的条件分布为二项分布，n=125,p由上面计算*/sita1=(ex2+x5)/(ex2+x3+x4+x5);/*M-步得到的迭代公式*/ifabs(sita1-sita)<=0.00001thenstop;output;sita=sita1;end;run;得到的迭代结果：obssita110.6082520.6243230.6264940.6267850.62682将sita=0.62682代入得到y1-y5的估计值：dataa;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.;puty1=y2=y3=y4=;run;y1=129y2=18y3=18y4=31。估计结果和实际值很按近。不用EM算法，直接估计时会分别得到4个sita的估计值：data;sita=4*(125/197-1/2);putsita=;sita=1-4*18/197;putsita=;sita=1-4*20/197;putsita=;sita=4*34/197;putsita=;run;得到sita估计值：sita=0.sita=0.sita=0.sita=0.计算用EM算法和直接估计得到的结果：dataa;dosita=0.6268,0.,0.,0.,0.;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.;puty1=y2=y3=y4=;output;end;run;结果显示：y1=129y2=18y3=18y4=31,EM算法的结果。y1=125y2=23y3=23y4=27y1=130y2=18y3=18y4=31y1=128y2=20y3=20y4=29y1=132y2=15y3=15y4=3417.2.2右截尾数据简单线性回归计算程序创建SAS数据集：dataa1;inputv1t1;cards;1701764170277217034441703542170378017048601705196190408190408190134419013441901440220408220408220504220504220504150806415080641508064150806415080641508064150806415080641508064150806417054481705448170544819016801901680190168019016801901680220528220528220528220528220528;run;按要求作数据变换，注意这里的条件n>17可以用其它的标识：dataa;seta1;v=1000/(v1+273.2);t=log10(t1);n=_n_;/*用于和后有参数估计的数据集合并*/vsq=v**2;/*用于求参数beta0,beta1和sigma估计*/by_v=1;/*为了以后和sw合并*/ifn>17thenc=t;dropv1t1;/*直接回归求得参数的初值，并将这些初值赋予宏变量beta01,beta11,sigma1*/procregdata=aoutest=estnoprint;modelt=v;dataest;setest;callsymput('beta01',intercept);/*创建一个值来自DATA步的宏变量beta01*/callsymput('beta11',v);/*创建一个值来自DATA步的宏变量beta11*/callsymput('sigma1',_rmse_);dataw;seta;beta01=&beta01;beta11=&beta11;sigma1=&sigma1;/*宏A求出迭代公式中的各项和，并得到迭代公式值，为下一步迭代提供值*/%macroA;dataw;setw;ifn>17thendoc=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;elsedot1=t;vt=v*t;ezv=0;ez=0;hq=0;tmu=(t-beta01-beta11*v)**2;end;procmeansdata=wnoprint;varvezezvvtt1hqvsqtmusigma1;outputout=swsum=sumvsumezsumezvsumvtsumt1sumhqsumvsqsumtmusumsigma1;datasw;setsw;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;keepbeta0beta1sigmaby_v;%mendA;/*将每次迭代的结果放在一个数据集result中，先放入直接回归得到参数估计的初值*/dataresult(keep=beta0beta1sigma);beta0=&beta01;/*第一个观测为初值*/beta1=&beta11;sigma=&sigma1;optionsnodatenonotesnosource;/*宏B是迭代程序*/%macrob;%doi=1%to30;%A;/*调用宏A*/dataw;mergeasw;byby_v;renamebeta0=beta01beta1=beta11sigma=sigma1;dataresult;/*将每次迭代的结果放在一个数据集*/setresultsw;%end;%mendb;%b;run;optionsnocenter;procprintdata=result;迭代结果作图：dataresult;setresult;n=_n_;procgplotdata=result;symbol1v=stari=joinc=blue;symbol2v=stari=joinc=black;symbol3v=stari=joinc=red;plotbeta0*n=1beta1*n=2sigma*n=3;run;直接回归结果：dataa;seta1;v=1000/(v1+273.