[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* **<=**