最尤法（回帰分析の例）数値計算（スコアリング法、１階導関数＋ヘッシアン行列の期待値必要)

プログラム

// Maximum Likelihood Estimation
// Numerical Result using Scoring method
#include<oxstd.h>
#include<oxfloat.h> // to use M_2PI
#import<maximize>   // to use maximization package

static decl s_mX, s_vY; // static variable to use in the functions
// Loglikelihood Function
fLoglik(const vP, const adFunc, const avScore, const amHess)
{
decl ck,cnobs;decl dsum,dsig2;decl vb,vg,ve;
 cnobs=rows(s_vY);ck=rows(vP);
 vb=vP[0:ck-2];dsig2=vP[ck-1];ve=s_vY-s_mX*vb;
 dsum=-sumsqrc(ve)/(2*dsig2);
 adFunc[0]=-0.5*cnobs*log(M_2PI*dsig2)+dsum;
if(avScore) // provides analytical 1st derivative
 {
 (avScore[0])[0:ck-2]=s_mX'*ve/dsig2;
 (avScore[0])[ck-1]=-0.5*cnobs/dsig2+sumsqrc(ve)/(2*dsig2^2);
 }
if(amHess) // provides expected value of hessian matrix
 {
 (amHess[0])=zeros(ck,ck);
 (amHess[0])[0:ck-2][0:ck-2]=-s_mX'*s_mX/dsig2;
 (amHess[0])[ck-1][ck-1]=-cnobs/(2*dsig2^2);
 }
return 1;
}
main()
{
decl ck,cnobs; // scalar (count, number of ...)
decl dfunc,ds2,dssr; // scalar (double)
decl minvhess,mi,msigma; // matrix
decl vb,ve,vp,vpvalue,vse,vyhat,vzvalue; // vector
 println("Maximum Likelihood Esimation using scoring method.");
 s_mX=<0;1;2;3;4>;  // X matrix: independent variable
 s_mX=ones(5,1)~s_mX; //          + constant term
 s_vY=<10;20;20;30;40>; // Y vector: dependent variable
 cnobs=rows(s_mX); // number of observations
 ck=columns(s_mX)+1; // number of parameters (including constant term & sigma^2)
 vp=<0;0;10>; // initial values
//
 println("initial values:",vp');
// Scoring Method 
 vp=<0;0;10>; // initial values
 MaxNewton(fLoglik,&vp,&dfunc,&minvhess,0);
// minvhess=invert(Hessian), Information matrix=-Hessian
// -mhess=invert(Information matrix)
 msigma=-minvhess;
 vse=sqrt(diagonal(msigma))';
 vzvalue=vp./vse; // z-value for H0: parameter=0
 vpvalue=tailn(fabs(vzvalue)); 
 println("%r",{"beta0","beta1","sigma^2"},
         "%c",{"est. coeff.", "std. err.","z-value", "p-value"},
         "%10.6",vp~vse~vzvalue~vpvalue);
}

計算結果

Maximum Likelihood Esimation using scoring method.
initial values:
      0.00000      0.00000       10.000

               est. coeff.    std. err.      z-value      p-value
beta0               10.000       1.8974       5.2705  6.8040e-008
beta1               7.0000      0.77460       9.0370  8.0544e-020
sigma^2             6.0000       3.7947       1.5811     0.056923