ポアソン回帰モデルの最尤推定　（プログラム）

#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,vm,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(vm.*exp(mx*vb));
 adFunc[0]+=sumc(vy.*log(vm)-loggamma(vy+1));
// 1st derivative vector
if(avScore){
 avScore[0]=sumc(vy.*mx)'-sumc(vm.*exp(mx*vb).*mx)';}
// 2nd derivative matrix, Hessian matrix
if(amHess){
 amHess[0]=-mx'*(vm.*exp(mx*vb).*mx);}
//
return 1;
}
//
main(){
decl cm,cnobs,dlhat,dzvalue;
decl vbhat,vl1hat,vse;
decl file,mdata,ml2hat;
// Read Data File
// X1-X8, m_i, y
cnobs=34;file = fopen("ship.txt");fscan(file,"%#m",cnobs,10,&mdata);fclose(file);
// dependent varible: yesvm
vy=mdata[][9]; // y
//independent variables: Constant X1-X8
mx=ones(cnobs,1)~mdata[][0:7];cm=columns(mx);
vm=mdata[][8]; // m
// Suppose that the prior for beta is set as follows: vb ~ N(0,1000I).
// To use improper prior, edit the subroutine fLK as mentioned above.  
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",{"B0","B1","B2","B3","B4","B5","B6","B7","B8"},
        "%c",{"Estimate","Std.Err.","z-value","p-value"},
        vbhat~vse~dzvalue~2*tailn(fabs(dzvalue)));      
}