Hi:
I modified the example provided in gsl manual in the section of
"Nonlinear Least-Squares Fitting". The original example have 40 terms
to minimize, now I change it to 2. Though the error produced is not
exactly the error from my program, I think the basically idea is the
gsl behavior when the number of terms is less than the number of
parameters: (N < p).
The error is as follows:
""""""
gsl: covar.c:51: ERROR: Jacobian be rectangular M x N with M >= N
Default GSL error handler invoked.
""""""
I suspect that gsl always expect N>= p, which might be reasonable. But
the code should inform the user earlier, using some checkers. Because
the error such as above is confusing for the users.
Nevertheless, in my own program, I add one more dummy term to raise N
up to p, which works quite well. Thanks for your work.
regards,
yalding
On 3/28/06, Brian Gough <[EMAIL PROTECTED]> wrote:
> Yajun Wang writes:
> > I want to use the Nonlinear Least Squares Fitting. I am working on the
> > sum of 3-dimensional formulas. However, the program quits quietly
> > whenever I have only 2 formulas to solve.
>
> Please can you send a example program which we can compile to
> reproduce the problem -- thanks.
>
> --
> Brian Gough
>
> Network Theory Ltd,
> Publishing the GSL Manual - http://www.network-theory.co.uk/gsl/manual/
>
#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_multifit_nlin.h>
struct data {
size_t n;
double * y;
double * sigma;
};
void
print_state (size_t iter, gsl_multifit_fdfsolver * s)
{
printf ("iter: %3u x = % 15.8f % 15.8f % 15.8f "
"|f(x)| = %g\n",
iter,
gsl_vector_get (s->x, 0),
gsl_vector_get (s->x, 1),
gsl_vector_get (s->x, 2),
gsl_blas_dnrm2 (s->f));
}
int
expb_f (const gsl_vector * x, void *params,
gsl_vector * f)
{
size_t n = ((struct data *)params)->n;
double *y = ((struct data *)params)->y;
double *sigma = ((struct data *) params)->sigma;
double A = gsl_vector_get (x, 0);
double lambda = gsl_vector_get (x, 1);
double b = gsl_vector_get (x, 2);
size_t i;
for (i = 0; i < n; i++)
{
/* Model Yi = A * exp(-lambda * i) + b */
double t = i;
double Yi = A * exp (-lambda * t) + b;
gsl_vector_set (f, i, (Yi - y[i])/sigma[i]);
}
return GSL_SUCCESS;
}
int
expb_df (const gsl_vector * x, void *params,
gsl_matrix * J)
{
size_t n = ((struct data *)params)->n;
double *sigma = ((struct data *) params)->sigma;
double A = gsl_vector_get (x, 0);
double lambda = gsl_vector_get (x, 1);
size_t i;
for (i = 0; i < n; i++)
{
/* Jacobian matrix J(i,j) = dfi / dxj, */
/* where fi = (Yi - yi)/sigma[i], */
/* Yi = A * exp(-lambda * i) + b */
/* and the xj are the parameters (A,lambda,b) */
double t = i;
double s = sigma[i];
double e = exp(-lambda * t);
gsl_matrix_set (J, i, 0, e/s);
gsl_matrix_set (J, i, 1, -t * A * e/s);
gsl_matrix_set (J, i, 2, 1/s);
}
return GSL_SUCCESS;
}
int
expb_fdf (const gsl_vector * x, void *params,
gsl_vector * f, gsl_matrix * J)
{
expb_f (x, params, f);
expb_df (x, params, J);
return GSL_SUCCESS;
}
#define N 2
int
main (void)
{
const gsl_multifit_fdfsolver_type *T;
gsl_multifit_fdfsolver *s;
int status;
size_t i, iter = 0;
const size_t n = N;
const size_t p = 3;
gsl_matrix *covar = gsl_matrix_alloc (p, p);
double y[N], sigma[N];
struct data d = { n, y, sigma};
gsl_multifit_function_fdf f;
double x_init[3] = { 1.0, 0.0, 0.0 };
gsl_vector_view x = gsl_vector_view_array (x_init, p);
const gsl_rng_type * type;
gsl_rng * r;
gsl_rng_env_setup();
type = gsl_rng_default;
r = gsl_rng_alloc (type);
f.f = &expb_f;
f.df = &expb_df;
f.fdf = &expb_fdf;
f.n = n;
f.p = p;
f.params = &d;
/* This is the data to be fitted */
for (i = 0; i < n; i++)
{
double t = i;
y[i] = 1.0 + 5 * exp (-0.1 * t)
+ gsl_ran_gaussian (r, 0.1);
sigma[i] = 0.1;
printf ("data: %d %g %g\n", i, y[i], sigma[i]);
};
T = gsl_multifit_fdfsolver_lmsder;
s = gsl_multifit_fdfsolver_alloc (T, n, p);
gsl_multifit_fdfsolver_set (s, &f, &x.vector);
print_state (iter, s);
do
{
iter++;
status = gsl_multifit_fdfsolver_iterate (s);
printf ("status = %s\n", gsl_strerror (status));
print_state (iter, s);
if (status)
break;
status = gsl_multifit_test_delta (s->dx, s->x,
1e-4, 1e-4);
}
while (status == GSL_CONTINUE && iter < 500);
gsl_multifit_covar (s->J, 0.0, covar);
#define FIT(i) gsl_vector_get(s->x, i)
#define ERR(i) sqrt(gsl_matrix_get(covar,i,i))
printf ("A = %.5f +/- %.5f\n", FIT(0), ERR(0));
printf ("lambda = %.5f +/- %.5f\n", FIT(1), ERR(1));
printf ("b = %.5f +/- %.5f\n", FIT(2), ERR(2));
{
double chi = gsl_blas_dnrm2(s->f);
printf("chisq/dof = %g\n", pow(chi, 2.0)/ (n - p));
}
printf ("status = %s\n", gsl_strerror (status));
gsl_multifit_fdfsolver_free (s);
return 0;
}
_______________________________________________
Bug-gsl mailing list
[email protected]
http://lists.gnu.org/mailman/listinfo/bug-gsl
