gsl_multifit_linear_est

[Giulio Bottazzi <bottazzi@sssup.it>]

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

hello

I've noticed that the function in subject is not currently
implemented in GSL. The following code is a possible (highly
non-optimal) candidate. I can send patch w.r.t. CVS if preferred.

	Giulio.


int
gsl_multifit_linear_est (const gsl_vector *x,
			 const gsl_vector *c,
			 const gsl_matrix *cov,
			 double *y, double *y_err)
{

  if (x->size != c->size)
    {
      GSL_ERROR ("number of parameters c does not match number of
observations x",                 GSL_EBADLEN);
    }
  else if (cov->size1 != cov->size2)
    {
      GSL_ERROR ("covariance matrix is not square", GSL_ENOTSQR);
    }
  else if (c->size != cov->size1)
    {
      GSL_ERROR
        ("number of parameters does not match size of covariance
matrix",         GSL_EBADLEN);
    }
  else{

    size_t dim = x->size;
    size_t i,j;
    double estimate,estimate_err;

    estimate=0.0;
    for(i=0;i<x->size;i++)
      estimate += gsl_vector_get (x,i)*gsl_vector_get (c,i);

    estimate_err=0.0;
    for(i=0;i<x->size;i++){
      const double dtmp1=gsl_vector_get (x,i);
      estimate_err += dtmp1*dtmp1*gsl_matrix_get(cov,i,i);
      for(j=0;j<i;j++)
	estimate_err += 2.*dtmp1*gsl_vector_get (x,j)*gsl_matrix_get(cov,i,j);
    }
    
    *y = estimate;
    *y_err =  sqrt(estimate_err);
    
    return GSL_SUCCESS;
  }
}





-- 
Giulio Bottazzi                 PGP Key ID:BAB0A33F 

CAFiM  Center for the Analysis of Financial Markets
LEM          Laboratory of Economics and Management
Sant'Anna School for Advanced Studies, Pisa, Italy
Phone: (+39)-050-883365       Fax: (+39)-050-883344
email: bottazzi@sssup.it    www.sssup.it/~bottazzi/

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