[SciPy-user] negative values in diagonal of covariance matrix

[josef.pktd at gmail.com]

On Fri, Dec 12, 2008 at 4:25 AM, Pauli Virtanen <pav at iki.fi> wrote:
> Fri, 12 Dec 2008 13:49:55 +0900, David Cournapeau wrote:
>> Do you mean the python wrapper miss some diagnostic information
>> available in fortran ? Otherwise, would it make sense to check the ier
>> value from fortran and at least generate a warning about failed
>> convergence ?
>
> The Python wrapper to minpack.lmder is quite thin, I'd expect any
> problems to be in the MINPACK code itself (which is actually LMDIF and
> not LMDER as I wrote earlier). The algorithm itself is something like
> Levenberg-Marquardt with a trust region.
>
> Unless the problem is in the minpack.leastsq code that forms the cov_x
> matrix from return values of LMDIF:
>
>        perm = take(eye(n),retval[1]['ipvt']-1,0)
>        r = triu(transpose(retval[1]['fjac'])[:n,:])
>        R = dot(r, perm)
>        cov_x = inv(dot(transpose(R),R))
>
> I don't have the time to check this now, though.
>
> --
> Pauli Virtanen
>
> _______________________________________________
> SciPy-user mailing list
> SciPy-user at scipy.org
> http://projects.scipy.org/mailman/listinfo/scipy-user
>

I compared the covariance that is returned by optimize.leastsq with
the covariance from the ols estimation using scipy.linalg.leastsq with
the linear in parameter model from the test suite.

The covariance in optimize.leastsq is in this case equal to:
(X'X)**(-1)  i.e. linalg.inv(dot(self.xx.T,self.xx))
which is not equal to the covariance of the estimated parameters,
which is  SSE/dof * (X'X)**(-1).

So at least the docstring of optimize.leastsq is ambiguous or misleading,
but except for the missing scale factor, the returned numbers are
essentially the same, in my example:

tls.inv_xx - optls.res_full[1]
array([[ -1.76561548e-11,   2.08572459e-10,  -3.03803875e-10],
       [  2.08572459e-10,  -2.37702299e-09,   3.39151827e-09],
       [ -3.03803876e-10,   3.39151827e-09,  -4.22841509e-09]])

estimated parameters are also ok:

>>> res_full[0] - tls.p_ols.T
array([[  1.01658158e-07,  -1.46458449e-06,   3.41746076e-06]])

I'm a bit rusty on non-linear least squares, and didn't check for
non-linear cases. Numpy/scipy seems to be missing a numerical hessian
calculation to verify non-linear optimization results (at least I
didn't find any).

For the original problem of this thread, only optimize.fmin is able to
get a reasonable estimate, and, as I expected, optimize.fmin_bfgs is
doing even worse than optimize.leastsq.

Josef