[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Bug-gsl] bug in gsl_multifit_wlinear
From: |
Richard Mott |
Subject: |
[Bug-gsl] bug in gsl_multifit_wlinear |
Date: |
Fri, 20 May 2005 14:19:27 +0100 |
User-agent: |
Mozilla Thunderbird 1.0 (Windows/20041206) |
Hi
There is a bug in gsl_multifit_wlinear (and presumably in
gsl_multifit_linear too). The problem manifests itself when the design
matrix X is not of full rank. The documentation says that this case is
handled correctly in the singular value decompostion by setting the tiny
SVD components to 0. In fact his doesn't happen. Instead the function
returns a huge and non-sensical chisq value. The block of code in
question is
for (j = 0; j < p; j++)
{
gsl_vector_view column = gsl_matrix_column (QSI, j);
double alpha = gsl_vector_get (S, j);
if (alpha != 0)
alpha = 1.0 / alpha;
gsl_vector_scale (&column.vector, alpha);
}
and needs to be changed to something like
int rank = 0;
for (j = 0; j < p; j++)
{
gsl_vector_view column = gsl_matrix_column (QSI, j);
double alpha = gsl_vector_get (S, j);
if (alpha >1.0e-7) {
rank ++;
alpha = 1.0 / alpha;
else {
alpha = 0.0;
gsl_vector_set(S,j,0.0);
}
gsl_vector_scale (&column.vector, alpha);
}
Finally, in the case of when the rank of the design matric is not full,
it is useful to return the true rank (given by the number of SVD
components that exceed the test >1.0e-7)
My choice of a cutoff of 1.0e-7 is farily arbitrary, but seems to work
My computer has the following uname:
zeon [289]% uname -a
Linux zeon 2.6.7 #1 SMP Wed Mar 2 15:28:51 GMT 2005 i686 GNU/Linux
Hope this helps
Best Wishes
Richard Mott
--
----------------------------------------------------
Richard Mott | Wellcome Trust Centre
tel 01865 287588 | for Human Genetics
fax 01865 287697 | Roosevelt Drive, Oxford OX3 7BN
- [Bug-gsl] bug in gsl_multifit_wlinear,
Richard Mott <=