右截尾数据线性回归EM算法_第1页
右截尾数据线性回归EM算法_第2页
右截尾数据线性回归EM算法_第3页
右截尾数据线性回归EM算法_第4页
右截尾数据线性回归EM算法_第5页
已阅读5页,还剩1页未读 继续免费阅读

下载本文档

版权说明:本文档由用户提供并上传,收益归属内容提供方,若内容存在侵权,请进行举报或认领

文档简介

精选优质文档-----倾情为你奉上精选优质文档-----倾情为你奉上专心---专注---专业专心---专注---专业精选优质文档-----倾情为你奉上专心---专注---专业例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.

温馨提示

  • 1. 本站所有资源如无特殊说明,都需要本地电脑安装OFFICE2007和PDF阅读器。图纸软件为CAD,CAXA,PROE,UG,SolidWorks等.压缩文件请下载最新的WinRAR软件解压。
  • 2. 本站的文档不包含任何第三方提供的附件图纸等,如果需要附件,请联系上传者。文件的所有权益归上传用户所有。
  • 3. 本站RAR压缩包中若带图纸,网页内容里面会有图纸预览,若没有图纸预览就没有图纸。
  • 4. 未经权益所有人同意不得将文件中的内容挪作商业或盈利用途。
  • 5. 人人文库网仅提供信息存储空间,仅对用户上传内容的表现方式做保护处理,对用户上传分享的文档内容本身不做任何修改或编辑,并不能对任何下载内容负责。
  • 6. 下载文件中如有侵权或不适当内容,请与我们联系,我们立即纠正。
  • 7. 本站不保证下载资源的准确性、安全性和完整性, 同时也不承担用户因使用这些下载资源对自己和他人造成任何形式的伤害或损失。

评论

0/150

提交评论