# How to use an optimization solver to get t-stats and p-values for the estimates?

I calculate a data log likelihood (evaluated at a set of parameters to be estimated), and my task is to find the set of parameters that maximize my log likelihood.

My problem is: thought there are a bunch of matlab optimization algorithms that can give me the estimated parameter (like fminserch, patternsearch), they don't give the t-stats (or standard error) of each estimate....So if I want to check the robustness of the estimation, what should I do?

Thank you :)

• These algorithms are optimization algorithms. They usually stop when they reach some local-minima. They do not really care about "robustness of the estimation" procedure because they are quite agnostic to it. I think what you want is actually the Hessian of your optimization function at your minimum... Apr 23 '15 at 5:39
• Thank you. But hessian matrix only tells me the 2nd derivatives at the optimum, how can I get t-stats from that?
– Ruby
Apr 23 '15 at 5:53
• I changed your title because I think you are actually asking a much more generic question than what your title currently conveys. Feel free to edit it back. Apr 23 '15 at 8:14
• Please see my answer below on how to get standard errors for your estimates using the Hessian matrix. From there getting $t$-values is just a fraction. Apr 23 '15 at 8:25
• Because you cross-posted this at stackoverflow.com/questions/29814239, we hope and expect that you will post a summary of the answers you receive there, or that you will flag one of your posts for migration to the other site, where the two sets of answers can be merged.
– whuber
Apr 23 '15 at 15:13

Assuming you have a function $f$ the Hessian matrix $H$ of $f$ at a particular point $\theta$ is is essentially the second-derivatives of that function at $\theta$. More importantly if now the function you optimize happens to be a log-likelihood the Hessian matrix $H$ is equal to the inverse of the covariance matrix $K$. This is why the inverse of the covariance matrix is often call precision matrix. Now if you want the standard errors you can get them immediately of each variable in $\theta$ you can get them directly by taking the square root of the diagonal elements in $K$ (where the diagonal elements $K_{i,i}$ are practically the variances of each $\theta_i$.

OK, so lets make this a bit more obvious using code. I base this example of linear regression for illustration purposes.

First we generate some data:

% Generate some data
N = 100;
rng(1234);
X = rand(N,2);
beta = [10,20];
epsilon = randn(100,1) *2;
y = epsilon + X*beta';


Then we use MATLAB's in-built function to fit our model so have some idea if we right or just completely off (always a good thing). We also quickly check our coefficient covariance:

goodLM = fitlm(X,y,'Intercept',false);
goodLM.CoefficientCovariance
% ans =
% 0.3164   -0.2338
% -0.2338   0.2910


Then we define the cost function we optimize against. Here as we have as simple linear model we optimize against the normal log-likelihood. Luckily for us MATLAB has got an in-built function for that (normlike). At this point we stop and remember that actually our linear model means that $y \sim N(X\beta, \sigma^2)$ and that similarly we assume that $\epsilon \sim N(0,\sigma^2)$. So we can ahead and optimize against the likelihood of our residuals coming from a normal with mean 0 and standard deviation $\sigma$.

myNlogLik =@(beta) normlike([0, std(y-X*beta')], y-X*beta');


We move ahead and optimize using a solver that actually offers us information about the Hessian (eg. fminunc). I set some starting values and I start my optimizer:

[x,fval,exitflag,output,grad,hessian] = fminunc(myNlogLik, [19,25]);


So what is the inverse of our Hessian?

inv(hessian)
% ans =
% 0.3103   -0.2289
% -0.2289   0.2854


Something quite close to our original matrix! :) If you compute the square roots of the diagonal of this matrix you will see that actually they are very close to the standard errors MATLAB gave originally (or in general the standard errors one would get using close form matrix equations).

[ sqrt(diag(inv(hessian))) goodLM.Coefficients.SE]
% ans =
% 0.5570    0.5625  % first column "our fits", second "MATLAB's"
% 0.5342    0.5395


Clearly as you are not going to get the same exact results for a couple of reasons. For starters one might not have converged to the same minimum. Second and more importantly one approximates the Hessian in some way; a stencil approach, some form of spline fitting, Richardson extrapolation? Pick your numerical poison...

• Thank you~! I tried your method, I probably don't quite understand what kind of information the hessian can deliver to us... Now I have a problem: my hessian is near singular, and calculated standard error for some paramters are zeros (some columns in the hessian are zeros acctually)... Does that mean I can NOT trust the estimation results? And the reasons may be that there is no variation in my data?
– Ruby
Apr 23 '15 at 14:36
• The information delivered by the Hessian is the local curvature of a function. In your particular case now as the local curvature is "flat" (so to speak) it means that your point estimate is more variable. This then translates to higher standard errors. You should trust your estimated results but you should present them for that they are: local optima obtained by a optimization procedure. As a result the standard errors of them will be possibly large. Try to alter the way the Hessian is computed (or even better compute yourself manually) this might give you better estimates. (cont.) Apr 23 '15 at 23:47
• The fact that a matrix (the Hessian in this case) has zeros all but guarantees that the inverse of it has the same sparsity pattern. Maybe you can regularize your Hessian (roughly speaking set all the negative or small eigenvalues of it to some arbitrary small $\lambda$). Check this thread for information on how to deal with an ill-condition matrix. Apr 23 '15 at 23:56
• MATLAB also has a relatively new equilibrate function for matrices Oct 23 '20 at 4:23

The generally preferred way is probably the bootstrap, or related methods like the jackknife. The basic idea is to resample your dataset a bunch of times and see how much the answer varies. This usually gives pretty good estimates, but can be computationally expensive. The matlab function bootci will compute this for you automatically.

If your optimization is convex and reasonably well-conditioned, you can probably get a good speedup on the bootstrap by doing it manually and warm-starting the optimization from previous solutions. If there are optimization problems, though, this will probably underestimate the true variation.

An alternative, when you can compute the gradient of your likelihood and the estimator is asymptotically normal, is to use the delta method.

There may be more specific methods if you tell us more about your model.

• Thank you Dougal. I will try the methods you suggest. The reason I want to calculate p-value is that, I guess my lilelihood function is ill-behaved, and there may be a lot of local optima, so the results are sensitive to a starting point. I just wannted to check whether the estimates are robust...
– Ruby
Apr 23 '15 at 5:55
• If there are lots of local optima, definitely don't use the bootstrapping shortcut I mentioned, and the delta method is also unlikely to be a good approximation. Apr 23 '15 at 5:56

One approach would be to use the observed Fisher information matrix (F), that is the Hessian (approximate curvature matrix) of the loglikelihood function at the estimated parameters.

The inverse of this gives you the variance-covariance matrix of the estimated parameters: V=inv(F).

The standard errors can be extracted from this: sd = sqrt(diag(V)).

Then you can use the asymptotic normality property of the maximum likelihood estimator to construct confidence intervals, to quantify the uncertainty of parameter estimation.

If you have a simple problem/model, e.g. estimating the parameters of some common distribution then there are matlab functions which provide you the variance-covariance matrix, e.g. [nlogL,V] = lognlike(params,data), similar functions are available for many other distributions.

If you have a more complex model you can numerically or symbolically calculate the Hessian and follow the steps outlined above.

• Hi Arpi, thank you! If my hessian is near singular (in matlab I can use pinv to get the inverse though).....does it mean the the standard errors calculated this way will be meaingless?
– Ruby
Apr 23 '15 at 13:26
• @Ruby, well I think it means something like that the standard errors are very large (~ infinite). You could construct profile likelihood curves to see how flat is the curve in the vicinity of the estimated parameters. Almost flat curve should correspond to the near singular Hessian. It means that the estimated parameters are very uncertain, I do not know whether it makes sense with your model and data. proflike Apr 23 '15 at 15:02
• Thank you@ Arpi. Well actually some standard errors are zeros...I will try the method you suggest, it does semm that my loglikelihood function is flat hence the results are unstable...
– Ruby
Apr 23 '15 at 15:08