最尤法（回帰分析の例）　解析的にもとめられる場合

プログラム

// Maximum Likelihood Estimation
// Analytical Result
#include<oxstd.h>
main()
{
decl ck,cnobs; // scalar (count, nnumber of...)
decl ds2,dssr; // scalar (double)
decl mi,msigma,mx; // matrix
decl vb,ve,vp,vpvalue,vse,vy,vyhat,vzvalue; // vector
 mx=<0;1;2;3;4>;  // x matrix: independent variable
 mx=ones(5,1)~mx; //          + constant term
 vy=<10;20;20;30;40>; // y vector: dependent variable
 cnobs=rows(mx); // number of observations
 ck=columns(mx)+1; // number of parameters (including constant term & sigma^2)
// use analytically obtained result
 vb=invert(mx'*mx)*mx'*vy; // beta
 vyhat=mx*vb; // y-hat
 ve=vy-vyhat; // residual
 dssr=ve'*ve; // Sum of Squared Residuals
 ds2=dssr/cnobs; // MLE:s^2
 vp=vb|ds2;   // parameter vector 
 // mi:information matrix
 mi=zeros(ck,ck);mi[0:ck-2][0:ck-2]=(1/ds2)*mx'*mx;
 mi[ck-1][ck-1]=cnobs/(2*ds2^2);
 msigma=invert(mi); // asymptotic covariance matrix
 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);
}

計算結果

               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