cancel
Showing results for 
Show  only  | Search instead for 
Did you mean: 
Choose Language Hide Translation Bar
ragnarl
Level II

How to compute the parameter std error (or standard deviation) from the Hessian returned by the Minimize() function?

Hello,

 

I have a problem with calculating the standard error of the parameter estimates using minimize() with the optional parameter <<showdetails which returns a few extra values including the Hessian Matrix.  

My (rather basic) understanding of the Hessian is that the diagonal elements should be a good approximation of the inverse of the variances. However, I do not get a good agreement between thus extracted values and the ones I get from simple linear regression. 

To illustrate here is an example:

x = [1.309, 1.471, 1.49, 1.565, 1.611, 1.68];
y = [2.138, 3.421, 3.597, 4.34, 4.882, 5.66];
b1 = 1; /*starting value for b1 */
b2 = 5; /*starting value for b2 */
sseExpr = Expr(
       Summation( i = 1, 6, (y[i] - (b1 +b2*x[i] )) ^ 2 )
);
//Use Minimize 
{objVal, iters, gradient, hessian}  = Minimize( sseExpr, {b1, b2} ,<<showdetails(true));
//Use Linear Regression
{Estimates, Std_Error, Diagnostics} = Linear Regression( y, x, <<printToLog );
min_stderror=sqrt(inverse(hessian)/nrows(x));

show("Parameter b1=",Estimates[1],b1);
show("Parameter b2",Estimates[2],b2);
show("Std_error in b1",std_error[1],min_stderror[1,1]);
show("Std_error in b2",std_error[2],min_stderror[2,2]);

The Result from this is shown below 

nParm=2  Newton   ******************************************************

Iter   nFree   Objective  RelGrad    NormGrad2  Ridge   nObj   nGrad   nHess   Parm0      Parm1      
0      2       128.6484   20788.83   4982.355   0       1      1       1       1          5          
1      2       0.074605   7.56e-25   1.88e-25   0       1      1       1       -10.427    9.489346   
Convergence SUCCESS: Gradient
Time: 0.0166666667209938


********************************************************************
************     Parameter Estimates      **************************
********************************************************************
Term           Estimates      Std Error      t Ratio        Prob>|t|       
********************************************************************
Intercept      -10.42696      0.72006        -14.48064      0.00013        
b1             9.48935        0.47199        20.10486       0.00004        
********************************************************************
RSquare:0.990201015051942
RSquare Adj:0.987751268814928

"Parameter b1=";
Estimates[1] = -10.4269614637334;
b1 = -10.4269614637335;
"Parameter b2";
Estimates[2] = 9.48934569169413;
b2 = 9.48934569169416;
"Std_error in b1";
std_error[1] = 0.720062429496937;
min_stderror[1,1] = 1.5220348519284;
"Std_error in b2";
std_error[2] = 0.471992539497802;
min_stderror[2,2] = 0.997676125760012;

While the parameter estimates are in perfect agreement, the errors are not. 

thus the question: How to compute the parameter std error (or standard deviation) from the Hessian returned by the Minimize() function?

 

Thanks for any assistance

 

 

1 ACCEPTED SOLUTION

Accepted Solutions
ragnarl
Level II

Re: How to compute the parameter std error (or standard deviation) from the Hessian returned by the Minimize() function?

A colleague helped me figure out what was wrong. 

The calculation of the standard error is incorrect. It should be:

{objVal, iters, gradient, hessian} = Minimize( sseExpr, {b1, b2} ,<<showdetails(true));
if(!isempty(objVal),
y_est=b1+b2*x;
stderror=sqrt(sum((y-y_est)^2)*inverse(hessian/2)/(nrows(x)-nrows(Estimates)));
);

 

which gives the same std error as the linear regression.

 

View solution in original post

1 REPLY 1
ragnarl
Level II

Re: How to compute the parameter std error (or standard deviation) from the Hessian returned by the Minimize() function?

A colleague helped me figure out what was wrong. 

The calculation of the standard error is incorrect. It should be:

{objVal, iters, gradient, hessian} = Minimize( sseExpr, {b1, b2} ,<<showdetails(true));
if(!isempty(objVal),
y_est=b1+b2*x;
stderror=sqrt(sum((y-y_est)^2)*inverse(hessian/2)/(nrows(x)-nrows(Estimates)));
);

 

which gives the same std error as the linear regression.