精選優(yōu)質(zhì)文檔-----傾情為你奉上精選優(yōu)質(zhì)文檔-----傾情為你奉上專心---專注---專業(yè)專心---專注---專業(yè)精選優(yōu)質(zhì)文檔-----傾情為你奉上專心---專注---專業(yè)例17.1的相關(guān)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;得到的迭代結(jié)果：obssita110.6082520.6243230.6264940.6267850.62682將sita=0.62682代入得到y(tǒng)1-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。估計結(jié)果和實際值很按近。不用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算法和直接估計得到的結(jié)果：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;結(jié)果顯示：y1=129y2=18y3=18y4=31,EM算法的結(jié)果。y1=125y2=23y3=23y4=27y1=130y2=18y3=18y4=31y1=128y2=20y3=20y4=29y1=132y2=15y3=15y4=3417.2.2右截尾數(shù)據(jù)簡單線性回歸計算程序創(chuàng)建SAS數(shù)據(jù)集：dataa1;inputv1t1;cards;1701764170277217034441703542170378017048601705196190408190408190134419013441901440220408220408220504220504220504150806415080641508064150806415080641508064150806415080641508064150806417054481705448170544819016801901680190168019016801901680220528220528220528220528220528;run;按要求作數(shù)據(jù)變換，注意這里的條件n>17可以用其它的標(biāo)識：dataa;seta1;v=1000/(v1+273.2);t=log10(t1);n=_n_;/*用于和后有參數(shù)估計的數(shù)據(jù)集合并*/vsq=v**2;/*用于求參數(shù)beta0,beta1和sigma估計*/by_v=1;/*為了以后和sw合并*/ifn>17thenc=t;dropv1t1;/*直接回歸求得參數(shù)的初值，并將這些初值賦予宏變量beta01,beta11,sigma1*/procregdata=aoutest=estnoprint;modelt=v;dataest;setest;callsymput('beta01',intercept);/*創(chuàng)建一個值來自DATA步的宏變量beta01*/callsymput('beta11',v);/*創(chuàng)建一個值來自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;/*將每次迭代的結(jié)果放在一個數(shù)據(jù)集result中，先放入直接回歸得到參數(shù)估計的初值*/dataresult(keep=beta0beta1sigma);beta0=&beta01;/*第一個觀測為初值*/beta1=&beta11;sigma=&sigma1;optionsnodatenonotesnosource;/*宏B是迭代程序*/%macrob;%doi=1%to30;%A;/*調(diào)用宏A*/dataw;mergeasw;byby_v;renamebeta0=beta01beta1=beta11sigma=sigma1;dataresult;/*將每次迭代的結(jié)果放在一個數(shù)據(jù)集*/setresultsw;%end;%mendb;%b;run;optionsnocenter;procprintdata=result;迭代結(jié)果作圖：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;直接回歸結(jié)果：dataa;seta1;v=1000/(v1+273.