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

#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]=sumc((vy.*mx))*vb-sumc(log(1+exp(mx*vb)));
// 1st derivative vector
if(avScore){
 avScore[0]=sumc(vy.*mx)'-sumc((exp(mx*vb)./(1+exp(mx*vb))).*mx)';}
// 2nd derivative matrix, Hessian matrix
if(amHess){
 amHess[0]=-mx'*((exp(mx*vb)./(1+exp(mx*vb)).^2).*mx);}
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);
fLK(vbhat,&dlhat,&vl1hat,&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)));      
}