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: the orthogonal matrix, q, produced by the QR facotrization of the final approximate Jacobi matrix, stored column wise. But scipy.optimize.leastsq calls the trusty MINPACK fortran routine lmder, who's comments say: c fjac is an output m by n array. the upper n by n submatrix c of fjac contains an upper triangular matrix r with c diagonal elements of nonincreasing magnitude such that c c p^t *(jac^t *jac)*p = r^t *r, c c where p is a permutation matrix and jac is the final c calculated jacobian. column j of p is column ipvt(j) c (see below) of the identity matrix. the lower trapezoidal c part of fjac contains information generated during c the computation of r. And a recent comment to this list said: On 10-Feb-05, at 11:23 PM, scipy-user-request@scipy.net wrote:
R. Padraic Springuel wrote:
Is there a way to get the uncertainty in the resultant parameters out of the optimize.leastsq fitting function?
optimize.leastsq? reveals that it can return a dictionary of optional values, among which is the Jacobian matrix fjac at the final step. The diagonal elements of fjac times its transpose can be taken as the squares of the errors in the corresponding coefficients.
So fjac is either: 1) the jacobian at the final step, 2) The Q of the QR factorization of the jacobian. 3) or the R of the QR factorization of the jacobian. Why is this important? Because the jacobian at the final step is needed to compute the covariance matrix, which gives the uncertainty in the fitted parameters, which is VERY important. So, any MINPACK experts out there? Can you explain to me what fjac is in clearer terms? Is the scipy.optimize.leastsq docstring accurate? -Brendan ----------------- PS, Here's a bit of linear algebra background: lmder is an implementation of the Levenberg-Marquardt algorithm. The levenberg marquardt method minimizes the sum of the squares of the residual functions by iteratively solving its derivative: ( J^t * J ) * x = J^t * f( x ), for x: x = ( J^t * J )^-1 * J^t * f( x ) where J is the Jacobian matrix (the derivatives of all the functions wrt each parameter x), x is the parameters to solve for, and f(x) is a vector of the residuals. (There's another bit here, where the diagonals of ( J^t * J ) are multiplied by a damping factor, but that's not important for the present problem) That first term on the right side of that last equation, ( J^t * J )^-1, is important, because *that's* the covariance matrix of the parameters x. Here's where I start to get confused. I *think* lmder solves the above function with QR factorization, so that the above equation becomes: x = ( R^t * Q^t * Q * R )^-1* R^t * Q^t * f( x ) = (R^t * R)^-1* R^t * Q^t * f( x ) = R^-1 * Q^T * f( x ) where Q, R are the results of the QR factorization of J. If lmder is indeed doing this (I've looked at the code, but my Fortran skills are poor), than I should be able to get the value of R at the last step, and calculate my covariance from (R^t * R)^-1 But at this point my head is hurting too much to know where to look for this matrix. Am I way off base?