[SciPy-user] negative values in diagonal of covariance matrix
Pavlo Shchelokovskyy
shchelokovskyy at gmail.com
Thu Dec 11 08:09:14 EST 2008
Hi all,
I'm a moderately new user of scipy, trying to make some curve-fitting
with it. I wanted to use cov matrix output of leastsq to estimate
errors of fitted parameters, but stumbled upon strange discrepancy
(for one particular dataset):
On Linux (Fedora 8), using Python 2.5.1, numpy 1.2.0 and scipy 0.6.0
>>> from scipy import *
>>> from scipy import optimize
>>> y = asarray([217, 182, 162, 170, 255])
>>> x = linspace(0, y.size - 1, y.size)
>>> gauss = lambda p, x: p[0] + p[1] * exp(-(x-p[2])**2/(2*p[3]**2))
>>> errfunc = lambda p, x, y: y - gauss(p, x)
>>> pinit = (max(y), -y.ptp(), argmin(y), y.size/4)
>>> fit,cov,info,mesg,success = optimize.leastsq(errfunc, pinit, args=(x,y), full_output = 1)
>>> fit
array([ 3.01602865e+05, -3.01444487e+05, 1.83283239e+00, 8.87273494e+01])
>>> cov
array([[ -2.27903048e+16, 2.27903047e+16, -4.72378743e+05, -3.35454486e+12],
[ 2.27903047e+16, -2.27903046e+16, 4.72378733e+05, 3.35454485e+12],
[ -4.72378950e+05, 4.72378947e+05, 6.38886491e-05, -6.95317302e+01],
[ -3.35454486e+12, 3.35454485e+12, -6.95316881e+01, -4.93761332e+08]])
I know that the fit is very bad (I can reject it afterwards in my
procedure), but what bothers me the most are the negative numbers on
the diagonal of the covariance matrix - as far as I understand, there
shouldn't be any... especially when the same script run on Windows XP,
using Python 2.5.2, numpy 1.2.0 and scipy 0.6.0 gives no such values
>>> fit
array([ 3.58426106e+05, -3.58267727e+05, 1.83283417e+00, 9.67300286e+01])
>>> cov
array([[ 6.40320190e+14, -6.40320188e+14, 4.24641077e+04, 8.64500286e+10],
[ -6.40320188e+14, 6.40320186e+14, -4.24641083e+04, -8.64500284e+10],
[ 4.24641250e+04, -4.24641094e+04, 7.64931901e-05, 5.73152924e+00],
[ 8.64500286e+10, -8.64500284e+10, 5.73152989e+00, 1.16716728e+07]])
Is it my error or misunderstanding somewhere, or is it really a bug in
Linux implementation?
Thanks in advance,
More information about the SciPy-User
mailing list