[Help-gsl] Re: once more fitting

[Thad Harroun]

> From: Petr Ent<address@hidden
> Subject: [Help-gsl] fitting hi
> everyone, is it possible to do following thing with non-linear least
> squares fitting...?
>
> double x_init[3] = { 1.0, 0.0, 0.0 }; gsl_vector_view x =
> gsl_vector_view_array (x_init, p); T = gsl_multifit_fdfsolver_lmsder;
>  s = gsl_multifit_fdfsolver_alloc (T, n, p);
>
> for ( int i = 0; i < 100 000; i++ )
>
> gsl_multifit_fdfsolver_set (s, &f, &x.vector);
>
> do { iter++; status = gsl_multifit_fdfsolver_iterate (s);
>
> if (status) break;
>
> status = gsl_multifit_test_delta (s->dx, s->x,1e-4, 1e-4); } while
> (status == GSL_CONTINUE && iter < 500);

In this, "x.vector" is your initial guess for the solution, and I'm 
pretty certain is unchanged during the fitting. The vector containing 
the parameters that are being changed is "s->x".

Why the gsl_vector_view? I would:

gsl_vector *guess = gsl_vector_alloc(p)
for(i=0, i<p, i++) gsl_vector_set(guess, i, x_init[i])
...
gsl_multifit_fdfsolver_set (s, &f, guess);

> end of for cycle

I can't understand the purpose of the outer "for" loop. As you have it 
now, if you restart the simulation with the same initial guess, it will 
reach the same solution 100,000 times. Are you trying to restart the 
fitting with new, different values?


> I have one more question about nonlinear least squares fitting, I
> will use the example from GSL online manual. I copied it to my
> program and just changed the equation to
>
> y = Ai(1-exp(-t/ti)) + Ap(1-exp(-t/tp)) with parameters Ai, Ap, ti,
> tp
>
> I set correct partial derivations too, there is no problem with this.
>
>
> if I let it be this way, everything works ok (I try few thousand
> iterations), each one converges to parameters I used to generate my
> data. But when I add some random noise to the source data, lets say
> only 1 percent of the value (values are from 400 to 600), solver
> cannot find solution (10 cases from 1000 tries) even when I give him
> nearly desired values as a starting guess. I traced the fitting
> proces and I have found that it is due its first iteration steps,
> when solver from near-ok initial guess makes some nonsense (parameter
> is set to 400 when generating equation, guess is 390 and after first
> iteration it is 1.xxE10). From this point, solver never returns to
> "way that converges to 400".

The LMDR and LMSDR routines assume that the Jacobian contains
the derivatives or the _residual_ with respect to the parameter. Wrong
values here and the inital step direction may be wrong.

My first guess is that there is a bug in the partial derivatives, or 
more likely, a bug in putting them into the correct matrix elements for 
the Jacobian. Make sure that the Jacobian values you calculate are not 
NaN or other unreasonably large number during the first or second step, 
and that each element has its correct value.

Second might be the uncertainty value of the data going into
the fit. I had a similar problem because I thought setting sigma_y = 1 
would mean the error in (x,y) would be ignored, but sigma_y was 
comparable to y or the noise, so the Jacobian was sending the fit 
solution in the wrong direction.


Thad
----------------------------------
Brock University
Department of Physics
St.Catharines, Ontario, L2S 3A1
  Canada
Phone: 905-688-5550 x4294