* gsl_multifit_linear_est
@ 2005-03-04 16:09 Giulio Bottazzi
0 siblings, 0 replies; only message in thread
From: Giulio Bottazzi @ 2005-03-04 16:09 UTC (permalink / raw)
To: gsl-discuss
[-- 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 --]
^ permalink raw reply [flat|nested] only message in thread
only message in thread, other threads:[~2005-03-04 16:09 UTC | newest]
Thread overview: (only message) (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
2005-03-04 16:09 gsl_multifit_linear_est Giulio Bottazzi
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).