[SciPy-user] optimize.leastsq -> fjac and the covariance redux
oliphant at ee.byu.edu
Sun Mar 13 22:55:32 CST 2005
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:
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,
More information about the SciPy-user