プロビットモデルの最尤推定　（プログラム）

#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 vprob,vb;vb=dX;
vprob=probn(mx*vb);
adFunc[0]=sumc(vy.*log(vprob)+(1-vy).*log(1-vprob));
return 1;
}
//
main(){
decl cm,cnobs,dlhat,dzvalue;
decl vbhat,vl1hat,vse;
decl file,mdata,ml2hat;
// 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];cm=columns(mx);
vbhat=zeros(cm,1);MaxBFGS(fLK,&vbhat,&dlhat,0,1);
// calculate Hessian matrix numerically.
Num2Derivative(fLK,vbhat,&ml2hat);
vse=sqrt(diagonal(invert(-ml2hat)))';
dzvalue=vbhat./vse;
println("log likelihood=",dlhat);
println("%r",{"Const","PRIV","SCHOOL","LOGINC","PTCON"},
        "%c",{"Estimate","Std.Err.","z-value","p-value"},
        vbhat~vse~dzvalue~2*tailn(fabs(dzvalue)));      
}