[SciPy-user] optimize.leastsq -> fjac and the covariance redux

[Travis Oliphant]

Brendan Simons wrote:

> I am scratching my brain trying to find the quickest way to construct 
> the covariance matrix corresponding to a least squares fit, and I've 
> come across some documentation quirks that are confusing me.
>
> The method scipy.optimize.leastsq returns, in its optional infodict 
> parameter, a value called fjac which according to the docstring is:

Brendan,

Thanks for your detailed comments.  I would love to see the covariance 
matrix as a standard optional output.  If you can help us get there from 
the code that would be great. 

As I look at the Fortran docstring, I think it's pretty clear that the 
leastsq Python docstring is wrong (sorry about that since I probably 
wrote it...)  fjac is basically the R matrix (once you apply some 
permutations), thus, I'm pretty sure you can use the upper triangular 
part of the fjac matrix with the ipvt vector (also in infodict) to 
caculate the covariance matrix as you describe.

So, I would construct a permutation matrix P (using ipvt) like this:

n = len(ipvt)
p = mat(take(eye(n),ipvt))

and then get R using some thing like

r = mat(MLab.triu(fjac))
R = r * p.T

Then covariance matrix estimate is

C = (R.T * R).I


This is totally untested, but I think you are on the right track to use 
fjac as the R matrix of the QR factorization.

Thanks for your help,

-Travis