Estimation of correlation matrix and confidence intervals for given parameter values

[tkeskita]

Hi,

Is there a way in lmfit to force estimation of covariance matrix and confidence intervals for a given set of fixed parameter values? I mean, without conducting parameter estimation? Looking at the docs for calculating confidence intervals, it seems that all of the examples there include a call to minimize() to carry out parameter estimation..? Sorry if I'm missing some obvious way to do this.

I'm aware that error estimates for unconverged parameter values are not reliable. I have two use cases where this would be useful nonetheless:

Case 1: Debugging models. When developing models, I sometimes want to get an estimate of the confidence intervals and correlation coefficient matrix for some given parameter values, when trying to debug cause for e.g. if model calculation fails, or if the minimizer seems to behave oddly.

Case 2: Identify that model has strongly correlated parameters. In some of my models, parameters are highly correlated from the get go, which would be beneficial to see early, before starting parameter estimation. This is also useful when trying to find out good initial values for parameters.

Thanks for any insight on this!

[newville]

@tkeskita

Is there a way in lmfit to force estimation of covariance matrix and confidence intervals for a given set of fixed parameter values? I mean, without conducting parameter estimation? Looking at the docs for calculating confidence intervals, it seems that all of the examples there include a call to minimize() to carry out parameter estimation..? Sorry if I'm missing some obvious way to do this.

Not really. When evaluating covariance and confidence intervals explicitly (see, for example, https://github.com/lmfit/lmfit-py/blob/master/lmfit/confidence.py), one moves each parameter away from its optimal value (typically by its expected 1-sigma uncertainty) and then optimizes all the other parameter values in response to that change. The covariance measures how parameter A responds when parameter B is changed from its optimal value. You have to make those measurements.

The magic of leastsq() is that it generates a very good approximation of the covariance matrix as part of its solution. While finding confidence intervals explicitly as with conf_interval() is clearly exploring parameter space, it is definitely slower than just running leastsq(). And the empirical evidence from the results shows that you have to have a pretty perverse problem, often with some very high correlations for there to be a noticeable difference between the resulting uncertainties. Now, if you are expecting to have more than 1 meaningful significant digit in your uncertainty, want to explore the asymmetry in the uncertainties, or are looking for accurate 3-sigma uncertainties, then this is probably worth running.

[reneeotten]

the code we use to calculate the covariance matrix using numdifftools is:

_covar_ndt = self._calculate_covariance_matrix(result.x)
if _covar_ndt is not None:
    result.covar = self._int2ext_cov_x(_covar_ndt, result.x)
    self._calculate_uncertainties_correlations()

the functions used here all part of the Minimizer class. So I'd imagine (but didn't try) that you could initialize the Minimizer and then use these with your current "best-fit" values to get your covariance matrix.

[tkeskita]

Thanks for answers! I don't know lmfit yet well enough to make the covariance calculation code above to work, but I found out that using Nelder-Mead with a large tolerance stops the optimizer at start values, so this sort of does what I want (using the doc_fitting_withreport.py example) although it's not very elegant solution:

fit_kws={'tol': 1e+308}
out = minimize(residual, fit_params, method='nelder', args=(x,), kws={'data': data}, **fit_kws)

[newville]

@tkeskita I strongly recommend not doing this. The usual interpretation of correlation implies that the parameters are at or near their optimum values. You are deliberately asking to explore correlations away from the optimal value.
The calculation of the uncertainties and correlation requires many calls to your residual function. So, you are already making several function calls and you are deliberately not allowing the parameter values to be optimized. I can believe -- and sympathize with -- those who might view anyone doing that as attempting to bias the "results", and would not be surprised to hear that some people view the fact that lmfit allows this as irresponsible on our part.