ロジットモデル
事後シミュレーション比較のプログラム
#include <oxstd.h> #import <maximize> // mx: independent variable, vy: dependent variable // We need to use them in the following subroutine. So set them static. static decl mx,vy; fLK(const dX, const adFunc, const avScore, const amHess){ // fLK will calculate the value of log posterior density as adFunc[0] decl vb;vb=dX; adFunc[0]=-vb'*vb/(2*1.0)+sumc((vy.*mx))*vb-sumc(log(1+exp(mx*vb))); // 1st derivative vector if(avScore){ avScore[0]=-vb/1.0+sumc(vy.*mx)'-sumc((exp(mx*vb)./(1+exp(mx*vb))).*mx)';} // 2nd derivative matrix, Hessian matrix if(amHess){ amHess[0]=-(1/1.0)*unit(rows(vb))-mx'*((exp(mx*vb)./(1+exp(mx*vb)).^2).*mx);} return 1; } // main(){ decl cburn,cm,cn,cnobs,crepeat; decl car_accept,car_count,cmh_accept,cmh_count; decl dh_n,dh_o,dlhat,dl_n,dl_o,dfrac; decl vb,vb1,vb_o,vb_n,veig,vbhat,vl1hat; decl file,mdata,mbeta_s,mB1,mB1root,ml2hat; // Read Data File // PUB12 PUB34 PUB5 PRIV YEARS SCHOOL LOGINC PTCON YESVM // We use PRIV SCHOOL LOGINC PTCON to explain YESVM //ranseed(111); cnobs=95; file = fopen("choice.txt");fscan(file,"%#m",cnobs,9,&mdata);fclose(file); // dependent varible: yesvm vy=mdata[][8]; // yesvm //independent variables: Constant PRIV SCHOOL LOGINC PTCON // // In posterior simulation comparison, convergence seems to be slow // due to the LOGINC & PTCON (their sample means are much larger than // those of other independent variables). So we subtract their mean // just for our convenience. // mx=ones(cnobs,1)~mdata[][3]~mdata[][5]~(mdata[][6]-meanc(mdata[][6])) ~(mdata[][7]-meanc(mdata[][7])); cm=columns(mx); // Suppose that the prior for beta is set as follows: vb ~ N(0,I). vbhat=zeros(cm,1);MaxBFGS(fLK,&vbhat,&dlhat,0,1); fLK(vbhat,&dlhat,&vl1hat,&ml2hat); // vbhat : posterior mode // dlhat : log posterior evaluated at mode (constant term ignored) // vl1hat: 1st derivative vector of log posterior evaluated at mode // ml2hat: 2nd derivative matrix of log posterior evaluated at mode // Hessian matrix mB1=-invert(ml2hat); // check whether mB1 is positive definite (just in case) // veig: eigenvalues of mB1. should be > 0 for positive definite matrix eigen(mB1,&veig); if(min(veig)<=0) {println("mB1 is not positive definite! Something may be wrong with the program.");} mB1root=choleski(mB1);vb1=mB1*(vl1hat+invert(mB1)*vbhat); // // Burn-in period: 10,000, Number of samples: 100,000 crepeat=100000;cburn=10000; mbeta_s=zeros(crepeat,cm); car_accept=car_count=cmh_accept=cmh_count=0; vb=vbhat; for(cn=-cburn;cn<crepeat;++cn){ // Calculate mode for different vy vbhat=zeros(cm,1); MaxNewton(fLK,&vbhat,&dlhat,0,0); fLK(vbhat,&dlhat,&vl1hat,&ml2hat); mB1=-invert(ml2hat); eigen(mB1,&veig); if(min(veig)<=0) {println("mB1 is not positive definite! Something may be wrong with the program.");} mB1root=choleski(mB1);vb1=mB1*(vl1hat+invert(mB1)*vbhat); if(fmod(cn,100)==0){println("cn=",cn);} // // A-R Step vb_o=vb; fLK(vb_o,&dl_o,0,0); dh_o=dlhat+vl1hat'*(vb_o-vbhat)+0.5*(vb_o-vbhat)'*ml2hat*(vb_o-vbhat); // generate new beta do{ vb_n=vb1+mB1root*rann(cm,1); fLK(vb_n,&dl_n,0,0); dh_n=dlhat+vl1hat'*(vb_n-vbhat)+0.5*(vb_n-vbhat)'*ml2hat*(vb_n-vbhat); car_count+=1; }while(ranu(1,1)>=exp(dl_n-dh_n)); car_accept+=1; // M-H Step dfrac=exp(dl_n+min(dl_o|dh_o)-dl_o-min(dl_n|dh_n)); if(ranu(1,1)<=dfrac){vb=vb_n;cmh_accept+=1;} cmh_count+=1; // Add the following lines for posterior sumlation comparison // Generate y|beta ~ Bernulli (p) where p=exp(x*beta)/(1+exp(x*beta)) decl dprob; dprob=exp(mx*vb)./(1+exp(mx*vb)); decl i; vy=zeros(cnobs,1);for(i=0;i<cnobs;++i){if(ranu(1,1)<=dprob[i]){vy[i]=1;}} // or you can use following lines // vy=floor(ranu(cnobs,1)-dprob+1.0); // if(cn>=0){mbeta_s[cn][]=vb';} } format(100); file = fopen("beta_psc.txt","w");fprint(file,"%15.10f",mbeta_s);fclose(file); // check whether the argorithm was efficient enough. println("AR acceptance=",car_accept/car_count); println("MH acceptance=",cmh_accept/cmh_count); }