I'm sorry to bother you folks with what should be simple, but I have been unable to get the sample code that executes the routine gsl_multilinear_wfit.c (which I literally cut and pasted from the on-line HTML documentation) to give an answer that is the same or even close to what is given in the documentation. Best as I can see this is not reported in the BUGS list or in the on-line bug reports, but I must admit to some difficulty navigating them.
The input data for the test is created using the program generate.c which
computes
y = exp(x) + noise
where 0.1 <= x <= 2.0. As the noise -> 0 one should recover the
Taylor expansion of exp(x), namely
exp(x) = 1.0 + x/1! + x^2/2! + ...
which is in fact what is shown in the documentation for this case
(allowing for the fact some considerable amount of noise has been
added).
My compiled code gives the result
# best fit: Y = 1.18246 + 0.184573 X + 1.3031 X^2
If I cut the noise size down in generate by setting
double sigma = 0.00001*y0 then I should get close to the analytic
values, but instead get
# best fit: Y = 1.08437 + 0.435157 X + 1.24313 X^2
I am running the code on a Ubuntu 10.03 LTS platform with an
AMD Phenom processor with GSL version 1.15. The same behavior (and
exactly the same answers are found on a 32 bit AMD XP2400 processor
with GSL 1.9 or GSL 1.15)
I have recompiled the library turning off all optimization in all
modules to no avail, same results.
So, having run out of ideas and having tested all the "standard"
kind of stuff I'm submitting this report. (BTW, I tried the test
example immediately preceding this one in the documentation and
it reproduces the result in the manual perfectly.) I'm attaching to
this report the following files:
1) generate.c - data generate (cut/pasted from the web)
2) exp.dat - output from above, input to fit
3) fit.c - data fitter (cut/pasted from the web)
4) fit.output - output from fit.c and exp.dat
Please let me know if there is any other information I can provide to
help explain this rather strange result.
--
Jim Heasley
Institute for Astronomy [email protected]
University of Hawaii phone: 808-956-6826
2680 Woodlawn Drive fax: 808-956-4532
Honolulu, HI 96822
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_randist.h>
int
main (void)
{
double x;
const gsl_rng_type * T;
gsl_rng * r;
gsl_rng_env_setup ();
T = gsl_rng_default;
r = gsl_rng_alloc (T);
for (x = 0.1; x < 2; x+= 0.1)
{
double y0 = exp (x);
double sigma = 0.1 * y0;
double dy = gsl_ran_gaussian (r, sigma);
printf ("%g %g %g\n", x, y0 + dy, sigma);
}
gsl_rng_free(r);
return 0;
}
exp.dat
Description: Netscape Proxy Auto Config
#include <stdio.h>
#include <gsl/gsl_multifit.h>
int
main (int argc, char **argv)
{
int i, n;
double xi, yi, ei, chisq;
gsl_matrix *X, *cov;
gsl_vector *y, *w, *c;
if (argc != 2)
{
fprintf (stderr,"usage: fit n < data\n");
exit (-1);
}
n = atoi (argv[1]);
X = gsl_matrix_alloc (n, 3);
y = gsl_vector_alloc (n);
w = gsl_vector_alloc (n);
c = gsl_vector_alloc (3);
cov = gsl_matrix_alloc (3, 3);
for (i = 0; i < n; i++)
{
int count = fscanf (stdin, "%lg %lg %lg",
&xi, &yi, &ei);
if (count != 3)
{
fprintf (stderr, "error reading file\n");
exit (-1);
}
printf ("%g %g +/- %g\n", xi, yi, ei);
gsl_matrix_set (X, i, 0, 1.0);
gsl_matrix_set (X, i, 1, xi);
gsl_matrix_set (X, i, 2, xi*xi);
gsl_vector_set (y, i, yi);
gsl_vector_set (w, i, 1.0/(ei*ei));
}
{
gsl_multifit_linear_workspace * work
= gsl_multifit_linear_alloc (n, 3);
gsl_multifit_wlinear (X, w, y, c, cov,
&chisq, work);
gsl_multifit_linear_free (work);
}
#define C(i) (gsl_vector_get(c,(i)))
#define COV(i,j) (gsl_matrix_get(cov,(i),(j)))
{
printf ("# best fit: Y = %g + %g X + %g X^2\n",
C(0), C(1), C(2));
printf ("# covariance matrix:\n");
printf ("[ %+.5e, %+.5e, %+.5e \n",
COV(0,0), COV(0,1), COV(0,2));
printf (" %+.5e, %+.5e, %+.5e \n",
COV(1,0), COV(1,1), COV(1,2));
printf (" %+.5e, %+.5e, %+.5e ]\n",
COV(2,0), COV(2,1), COV(2,2));
printf ("# chisq = %g\n", chisq);
}
gsl_matrix_free (X);
gsl_vector_free (y);
gsl_vector_free (w);
gsl_vector_free (c);
gsl_matrix_free (cov);
return 0;
}
0.1 1.11997 +/- 0.110517 0.2 1.21064 +/- 0.12214 0.3 1.57588 +/- 0.134986 0.4 1.60127 +/- 0.149182 0.5 1.81319 +/- 0.164872 0.6 1.58934 +/- 0.182212 0.7 1.53111 +/- 0.201375 0.8 2.07436 +/- 0.222554 0.9 2.44999 +/- 0.24596 1 2.96118 +/- 0.271828 1.1 2.99886 +/- 0.300417 1.2 2.88964 +/- 0.332012 1.3 3.42419 +/- 0.36693 1.4 4.12889 +/- 0.40552 1.5 4.85414 +/- 0.448169 1.6 4.68148 +/- 0.495303 1.7 5.12469 +/- 0.547395 1.8 6.05394 +/- 0.604965 1.9 6.23804 +/- 0.668589 # best fit: Y = 1.18246 + 0.184573 X + 1.3031 X^2 # covariance matrix: [ +1.25612e-02, -3.64387e-02, +1.94389e-02 -3.64387e-02, +1.42339e-01, -8.48761e-02 +1.94389e-02, -8.48761e-02, +5.60243e-02 ] # chisq = 16.0065
