[SciPy-user] Re: optimize.leastsq -> fjac and the covariance redux
Brendan Simons
brendansimons at yahoo.ca
Mon Mar 14 13:07:21 CST 2005
>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
Travis,
I too would like to see the covariance in the leastsq
optional outputs. Using your suggestion, I whipped up
the following script, which fits a simple quadratic
(chosen because I can calculate the jacobian easily by
hand):
#---begin script--------
"""script to test the fjac output matrix of
scipy.optimize.leastsq
Placed in public domain by Brendan Simons
March 2005
"""
from scipy import *
import MLab
from scipy.optimize import leastsq
x = array([0., 1., 2., 3., 4., 5.], typecode='f')
coeffs = [4., 3., 5.]
yErrs = array([0.1, 0.2, -0.1, 0.05, 0, -.02],
typecode='f')
m = len(x)
n = len(coeffs)
def evalY(p, x):
"""a simple quadratic"""
return p[2]*x**2 + p[1]*x + p[0]
def J(p, x):
"""the jacobian of a quadratic"""
result = zeros([m,n], typecode='f')
result[:, 0] = ones(m, typecode='f')
result[:, 1] = x
result[:, 2] = x*x
return result
def resid(p, y, x):
return y - evalY(p, x)
y_true = y(coeffs, x)
y_meas = y_true + yErrs
p0 = [3., 2., 4.]
pMin, infodict, ier, mesg = leastsq(resid, p0,
args=(y_meas, x), full_output=True)
fjac = infodict['fjac']
ipvt = infodict['ipvt']
#here's travis' suggestion
perm = mat(take(eye(n),ipvt-1))
r = mat(MLab.triu(fjac[:,:n]))
R = r * perm.T
C = (R.T * R).I
#check against know algorithms for computing C
Q_true, R_true = linalg.qr(J(pMin, x))
C_true = linalg.inv(dot(transpose(J(pMin, x)),
J(pMin, x)))
C_true2 = linalg.inv(dot(transpose(R_true), R_true))
print 'pMin'
print pMin
print '\nfjac'
print fjac
print'\nipvt'
print ipvt
print '\nperm'
print perm
print '\nr'
print r
print '\nR'
print R
print '\nR_true'
print R_true
print '\nC'
print C
print '\nC_true'
print C_true
print '\nC_true2'
print C_true2
#---end script--------
Note that I had to make the following modifications to
your algorithm:
a) The Fortran docstring says r is in the first nxn
submatrix of fjac, so I truncated fjac with [:, :n]
before applying MLab.triu
b) the identity matrix columns returned by ipvt are
1-based, not 0-based (Fortran is 1-based I guess). So
I subtracted 1 in
which results in perm = mat(take(eye(n),ipvt-1))
The results from the script are as follows:
#--- begin script output----
pMin
[ 4.13714286 2.93428571 5.00714286]
fjac
[[-31.28897539 -0.03196014 -0.12784056 -0.28764125
-0.51136222 -0.79900346]
[ -7.19103119 1.81357942 0.59589038 0.51365984
0.17797854 -0.41115309]
[ -1.75780755 1.30104626 -1.10335453 0.03523648
-0.29732131 -0.95308465]]
ipvt
[3 2 1]
perm
Matrix([[ 0., 0., 1.],
[ 0., 1., 0.],
[ 1., 0., 0.]])
r
Matrix([[-31.28897539, -0.03196014, -0.12784056],
[ 0. , 1.81357942, 0.59589038],
[ 0. , 0. , -1.10335453]])
R
Matrix([[ -0.12784056, -0.03196014, -31.28897539],
[ 0.59589038, 1.81357942, 0. ],
[ -1.10335453, 0. , 0. ]])
R_true
[[ -2.44948983 -6.12372446 -22.45365715]
[ 0. 4.18329954 20.916502 ]
[ 0. 0. 6.11010218]
[ 0. 0. 0. ]
[ 0. 0. 0. ]
[ 0. 0. 0. ]]
C
Matrix([[ 8.21428627e-001, -2.69897976e-001,
-3.08050728e-003],
[-2.69897976e-001, 3.92718047e-001,
7.01607645e-004],
[-3.08050728e-003, 7.01607645e-004,
1.03332016e-003]])
C_true
[[ 0.8214277 -0.58928472 0.08928555]
[-0.58928496 0.72678488 -0.13392843]
[ 0.08928559 -0.13392843 0.02678569]]
C_true2
[[ 0.82142925 -0.58928621 0.08928578]
[-0.58928627 0.72678584 -0.13392855]
[ 0.08928579 -0.13392855 0.0267857 ]]
#--- end script output----
So you can see it doesn't work as expected. C does not
equal C_true. More troublesome is that *no* columns
of fjac resemble R_true, no matter what permutation
matrix is applied. At this point I'm lost as to what
is what. Any suggestions?
-Brendan
--
PS, If and when we figure out the answer, I'd like to
submit the change to CVS. However, never having done
such a thing, I will need some instruction.
______________________________________________________________________
Post your free ad now! http://personals.yahoo.ca
More information about the SciPy-user
mailing list