ロジットモデル

事後シミュレーション比較のプログラム

#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);
}