最尤法（回帰分析の例）数値計算（ニュートン法、１階導関数必要)

プログラム

// Maximum Likelihood Estimation
// Numerical Result using Newton method (1st derivative)
#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);
 }
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 Newton method with 1st derivative.");
 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');
// Newton Method with 1st derivative only.
 MaxNewton(fLoglik,&vp,&dfunc,&minvhess,1);
// Numerical 2nd derivatives are used.
// 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);
// Use BHHH estimator instead
decl mg;
 println("Use BHHH estimates instead");
 mg=zeros(ck,cnobs);
 mg[0:ck-2][]=s_mX'.*(s_vY-s_mX*vp[0:1])'/vp[2];
 mg[ck-1][]=-0.5/vp[2]+((s_vY-s_mX*vp[0:1]).^2)'/(2*vp[2]^2);
 msigma=invert(mg*mg');
 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 Newton method with 1st derivative.
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

Use BHHH estimates instead

               est. coeff.    std. err.      z-value      p-value
beta0               10.000       2.5962       3.8519  5.8613e-005
beta1               7.0000       1.2000       5.8333  2.7165e-009
sigma^2             6.0000       6.1188      0.98058      0.16340