プロビットモデル

プログラム

#include <oxstd.h>
main()
{
decl cburn,cm,cn,cnobs,crepeat;
decl file,mdata,mbeta_s,mx,vy,vystar;
decl mB0,mB1;
decl i,vb0,vb1,vb;
// Read Data File
// PUB12        PUB34   PUB5    PRIV    YEARS   SCHOOL  LOGINC  PTCON   YESVM
// We use PRIV  SCHOOL  LOGINC  PTCON   to explain YESVM
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
mx=ones(cnobs,1)~mdata[][3]~mdata[][5:7];
//
// prior for beta vb ~ N(vb0,mB0). vb0=0,mB0=100*I
cm=columns(mx);vb0=zeros(cm,1);mB0=100*unit(cm);
// initialize latent variable: vystar
vystar=zeros(cnobs,1);
//
mB1=invert(mx'*mx+invert(mB0));
//       Burn-in period: 2000, Number of repetition: 10000
cburn=3000;crepeat=10000;mbeta_s=zeros(crepeat,cm);
for(cn=-cburn;cn<crepeat;++cn){
// vb | ystar  ~ N(vb1,mB1)
vb1=mB1*(mx'*vystar+invert(mB0)*vb0);
vb=vb1+choleski(mB1)*rann(cm,1);
// ystar | vb ~TN
 for(i=0;i<cnobs;++i){
  if(vy[i]==0){
     do{vystar[i]=mx[i][]*vb+rann(1,1);
           }while(vystar[i] >=0);
   }
  else{
     do{vystar[i]=mx[i][]*vb+rann(1,1);
           }while(vystar[i] <0);
   }
 }
// save samples for vb
//println(vb');
println(cn);
if(cn>=0){mbeta_s[cn][]=vb';}
}
 format(100);
 file = fopen("beta.txt","w");fprint(file,"%15.10f",mbeta_s);fclose(file);
}