最尤法（回帰分析の例）　数値計算（BFGS法、導関数不要)

プログラム

// Maximum Likelihood Estimation
// Numerical Result using BFGS
#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;
cnobs=rows(s_vY);ck=rows(vP);
vb=vP[0:ck-2];dsig2=vP[ck-1];
dsum=(sumsqrc(s_vY-s_mX*vb))/(2*dsig2);
adFunc[0]=-0.5*cnobs*log(M_2PI*dsig2)-dsum;
return 1;
}
main()
{
decl ck,cnobs; // scalar (count, number of ...)
decl dfunc,ds2,dssr; // scalar (double)
decl mhess,mi,msigma; // matrix
decl vb,ve,vp,vpvalue,vse,vyhat,vzvalue; // vector
 println("Maximum Likelihood Esimation using BFGS 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');
 MaxBFGS(fLoglik,&vp,&dfunc,0,TRUE);
 Num2Derivative(fLoglik,vp,&mhess); // Numerical Hessian Matrix
 mi=-mhess;msigma=invert(mi);
 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 BFGS 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.8037e-008
beta1               7.0000      0.77459       9.0370  8.0519e-020
sigma^2             6.0000       3.7947       1.5811     0.056922