Possible inconsistency in gsl_multifit_covar

Possible inconsistency in gsl_multifit_covar

[Giulio Bottazzi]

[-- Attachment #1: Type: text/plain, Size: 785 bytes --]

Hello,
I've noticed that the matrix returned by the
function "gsl_multifit_covar" is defined as

(J^T J)^{-1}


Now the problem is that this matrix IS NOT the
(approximated) estimated variance-covariance matrix of the asymptotic
normal distribution of the estimates. This should indeed be defined as

\sigma^2 (J^T J)^{-1}

where

\sigma^2 = \sum_j res_j / (N-K)

res being the residuals, N-K the degrees of freedom.

This behavior of gsl_multifit_covar is strange, first because calling a
"covariance matrix" what the function returns sound a little bit
deceiving and second because it is inconsistent with the covariance
matrix returned by the GSL linear routines (which contain the proper
factor \sigma^2). 

Is this a feature or a bug?

Best,
	Giulio.

[-- Attachment #2: Type: application/pgp-signature, Size: 189 bytes --]

Re: Possible inconsistency in gsl_multifit_covar

[Brian Gough]

Giulio Bottazzi writes:
 > I've noticed that the matrix returned by the
 > function "gsl_multifit_covar" is defined as
 > 
 > (J^T J)^{-1}
 > 
 > 
 > Now the problem is that this matrix IS NOT the
 > (approximated) estimated variance-covariance matrix of the asymptotic
 > normal distribution of the estimates. This should indeed be defined as
 > 
 > \sigma^2 (J^T J)^{-1}
 > 

Hello,

The multifit functions are a translation of the corresponding MINPACK
routines, which defined 'covar' as in GSL.  In MINPACK the definition
of J includes a factor of sigma 

  J_ij = dF_i/dx_j
  F_i = (y_i-y(x))/sigma_i

Does this account for the difference above?

-- 
Brian Gough

Re: Possible inconsistency in gsl_multifit_covar

[Giulio Bottazzi]

[-- Attachment #1: Type: text/plain, Size: 2464 bytes --]

On Tue, 23 Aug 2005 15:29:29 +0100
Brian Gough <bjg@network-theory.co.uk> wrote:

> Hello,
> 
> The multifit functions are a translation of the corresponding MINPACK
> routines, which defined 'covar' as in GSL.  In MINPACK the definition
> of J includes a factor of sigma 
> 
>   J_ij = dF_i/dx_j
>   F_i = (y_i-y(x))/sigma_i
> 
> Does this account for the difference above?
> 
> -- 

Hello,
as far as I understand, given that I do not know MINPACK, I would say
yes and no. The point is subtle and has to do with the double naure of
the least-quare procedure: a generic tools for "fitting things" and a
maximum likelihood estimator of a well defined class of statistical
models. I hope you'll forgive me if I try to elaborate the issue a
little bit below.

From a statistical point of view, you can assume that your data follow
the model

(I)   y_i ~ f(x_i,B) + \w_i * \epsilon_i

where \epsilon_i ~ N(0,1) (i.i.d. normally distributed with mean zero
and variance 1) or instead the model

(II)  y_i ~ f(x_i,B) + \w_i * \epsilon_i

where \epsilon_i ~ N(0,\sigma) (i.i.d. normally distributed with mean
zero and unknon variance). Both prescriptions allow for the
specification of weights w_i, but differ on the number of the unknowns.
In particular, the variance of the distribution \sigma is estimated from
data in Model II while it is assumed known in Model I.

Now in GSL the model assumption for the different linear routines
(f(x_i,B) = x_i . B) are:


function:            assumption:


gsl_fit_linear         Model II

gsl_fit_wlinear        Model I

gsl_multifit_linear    Model II

gsl_multifit_wlinear   Model I


which is consistent (even if, AFAIK, undocumented): if weights are
provided, than the variance is assumed know, if they are not, it is
estimated.

The possible inconsistency is with the non-linear function

gsl_multifit_covar   Model I

now weights are NOT explicitly required, but nonetheless the variance is
assumed known. For comparison, consider the analogous function in NAG C
library, "e04ycc", which instead assume Model II.

Probably it would help to explicitly mention the formula used to compute
variance-covariance matrix in ALL the routines, so that the average
dumb user (like me), by comparing the different formulas, can immediatly
understand were differences can possibly arise. What do you think?

Apologies for the lenght of the message.

	Giulio.

[-- Attachment #2: Type: application/pgp-signature, Size: 189 bytes --]

Re: [Bug-gsl] Covariance estimate in weighted regression

[Brian Gough]

Alex Tartakovsky writes:
 > >From GSL manual (pp. 361-362), standard texts, and just common
 > >sense, one expects that the output produced by a weighted
 > >regression with all the weights set to 1 should be the same as
 > >from unweighted regression.  This is not the case for the
 > >covariance estimates produced by "fit" and "multifit"
 > >least-squares GSL functions.  The reason is that the cov estimates
 > >in the straight versions of the functions include s2 (an estimate
 > >of the error variance):

Giulio Bottazzi writes:
 > Probably it would help to explicitly mention the formula used to
 > compute variance-covariance matrix in ALL the routines, so that the
 > average dumb user (like me), by comparing the different formulas,
 > can immediatly understand were differences can possibly arise. What
 > do you think?

Thanks for the comments, I have added some longer explanations in the
manual about how the covariance matrices are computed and their
definitions for the different cases.  The new chapters are available
at http://www.network-theory.co.uk/download/gsl/newchaps.ps.gz

-- 
Brian Gough

Network Theory Ltd,
Publishing Free Software Manuals --- http://www.network-theory.co.uk/