[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Help-gsl] Re: once more fitting
From: |
Thad Harroun |
Subject: |
[Help-gsl] Re: once more fitting |
Date: |
Wed, 01 Nov 2006 13:26:39 -0500 |
User-agent: |
Thunderbird 1.5.0.7 (Windows/20060909) |
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
- [Help-gsl] Re: once more fitting,
Thad Harroun <=