[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re:[Help-gsl] Re: once more fitting
|
From: |
Petr Ent |
|
Subject: |
Re:[Help-gsl] Re: once more fitting |
|
Date: |
Thu, 02 Nov 2006 10:03:33 +0100 (CET) |
hello,
to the "outer for loop" - I am sorry, I forgot to add a line with input data change,
it was meant to be as simple as possible :). I just wanted to tell that "when I run this
loop for the first time, it works as it is supposed to be. But when I run it (once) more at
same program run, I got NaNs in x vector.
I am quite sure with partial derivations I use, but I would better place some
code here to show you how I do that..... as I said, I took the sample program
from GSL manual and just replace equation to fit and its partial derivations
int TLogFit::fit(const gsl_vector *x,gsl_vector *f)
{
double Ai = gsl_vector_get(x,0);
double ti = gsl_vector_get(x,1);
double Ap = gsl_vector_get(x,2);
double tp = gsl_vector_get(x,3);
int n = numberPoint();
for ( int i = 0; i < n; i++ )
{
double t = i;
double v = 100 + Ai*(1-exp(-t/ti)) + Ap*(1-exp(-t/tp));
double y = Y->at(i);
double r = v - y;
gsl_vector_set(f,i, r);
}
return 0;
}
this is method where the value of function is computed. I call it indirectly,
so there is one parameter missing, but thats fine.
int TLogFit::fitd(const gsl_vector *x,gsl_matrix *m)
{
double Ai = gsl_vector_get(x,0);
double ti = gsl_vector_get(x,1);
double Ap = gsl_vector_get(x,2);
double tp = gsl_vector_get(x,3);
int n = numberPoint();
for ( int i = 0; i < n; i++ )
{
double t = i;
double ei = exp(-t/ti);
double ep = exp(-t/tp);
gsl_matrix_set (m, i, 0, 1-ei);
gsl_matrix_set (m, i, 1, -Ai*ei*t/(ti*ti));
gsl_matrix_set (m, i, 2, 1-ep);
gsl_matrix_set (m, i, 3, -Ap*ep*t/(tp*tp));
}
return 0;
}
here is partial derivations method. I asked many people to check these
equations, so I hope they are correct.
f = 100 + Ai*(1-exp(-t/ti)) + Ap*(1-exp(-t/tp))
(100 is just some offset, I add it when I generate data, so I have to add it
here too. But I got same results with/without it.)
f'(Ai) = 1 - exp(-t/ti)
f'(ti) = -Ai * exp(-t/ti) * t/ti^2
f'(Ap) = 1 - exp(-t/tp)
f'(tp) = -Ap * exp(-t/tp) * t/tp^2
I checked partial derivations in first step and they range from +-0.02 to
+-0.9, and here is the log with computation steps
Ai: 90.000000 ti: 40.000000 Ap: 60.000000 tp: 30.000000
Ai: 8831.107932 ti: -112.688446 Ap: -8510.709743 tp: -967.486498
Ai: 854.019667 ti: 527.461525 Ap: -482.653284 tp: -484.379713
Ai: 148.817022 ti: 328.201707 Ap: 90.790375 tp: -283.839888
Ai: 76.927852 ti: 14.208275 Ap: 57.712913 tp: -22.750787
Ai: 93.710786 ti: 37.245677 Ap: 63.181652 tp: 26.916838
Ai: 100.409324 ti: 31.566805 Ap: 69.085187 tp: 20.118399
Ai: 90.873797 ti: 27.168735 Ap: 64.522027 tp: -2.468986
Ai: 101.502109 ti: 30.425712 Ap: 70.113333 tp: 18.552237
Ai: 103.170378 ti: 28.247236 Ap: 71.835084 tp: 14.989352
Ai: 98.827858 ti: 28.155186 Ap: 69.477609 tp: 6.000933
Ai: 83.070414 ti: 13.166634 Ap: 62.498510 tp: 7.030133
Ai: 94.607768 ti: 25.471029 Ap: 69.554715 tp: 6.563325
Ai: 86.854841 ti: 19.664240 Ap: 68.123561 tp: 7.203937
Ai: 98.480305 ti: 14.671131 Ap: 52.022971 tp: 7.286642
Ai: 85.468642 ti: 17.898347 Ap: 67.357410 tp: 7.859267
Ai: 88.879871 ti: 16.719137 Ap: 62.951567 tp: 7.981395
Ai: 92.928240 ti: 16.438553 Ap: 58.843418 tp: 7.830807
Ai: 100.930343 ti: 15.755671 Ap: 50.579582 tp: 7.448265
Ai: 116.962902 ti: 14.585029 Ap: 34.057786 tp: 6.563051
Ai: 132.976250 ti: 13.683968 Ap: 17.553396 tp: 5.110763
Ai: 121.393902 ti: 14.440941 Ap: 29.542000 tp: 5.978224
Ai: 130.156638 ti: 13.832772 Ap: 20.413111 tp: 5.124316
Ai: 146.713814 ti: 12.691763 Ap: 2.998488 tp: 2.688985
Ai: 131.882607 ti: 13.743564 Ap: 18.623268 tp: 4.633074
Ai: 135.362803 ti: 13.485765 Ap: 14.941807 tp: 4.054630
Ai: 141.999370 ti: 12.952887 Ap: 7.815718 tp: 2.661430
Ai: 146.602055 ti: 12.449263 Ap: 2.564257 tp: 0.095629
Ai: 143.995645 ti: 12.712751 Ap: 5.486675 tp: 1.068644
Ai: 144.375252 ti: 12.592269 Ap: 4.891963 tp: 0.420256
Ai: 144.320518 ti: 12.591193 Ap: 4.938670 tp: 0.320396
Ai: 144.361350 ti: 12.584322 Ap: 4.886750 tp: 0.246665
Ai: 144.386895 ti: 12.579999 Ap: 4.854196 tp: 0.166695
Ai: 144.397113 ti: 12.578251 Ap: 4.841125 tp: -0.005567
Ai: 144.381139 ti: 12.581007 Ap: 4.861781 tp: 0.151074
Ai: 144.379198 ti: 12.581273 Ap: 4.863950 tp: 0.116240
Ai: 144.378210 ti: 12.581434 Ap: 4.865180 tp: 0.081526
Ai: 144.377750 ti: 12.581511 Ap: 4.865764 tp: 0.046855
Ai: 144.377729 ti: 12.581515 Ap: 4.865791 tp: 0.078058
Ai: 144.377718 ti: 12.581517 Ap: 4.865805 tp: 0.071448
Ai: 144.377718 ti: 12.581517 Ap: 4.865805 tp: 0.064501
Ai: 144.377717 ti: 12.581517 Ap: 4.865807 tp: 0.070753
Ai: 144.377717 ti: 12.581517 Ap: 4.865806 tp: 0.069364
Ai: 144.377717 ti: 12.581517 Ap: 4.865806 tp: 0.066585
Ai: 144.377717 ti: 12.581517 Ap: 4.865806 tp: 0.069086
Ai: 144.377717 ti: 12.581517 Ap: 4.865807 tp: 0.069336
Ai: 144.377717 ti: 12.581517 Ap: 4.865807 tp: 0.069280
Ai: 144.377717 ti: 12.581517 Ap: 4.865806 tp: 0.069169
Ai: 144.377717 ti: 12.581517 Ap: 4.865807 tp: 0.069269
Ai: 144.377717 ti: 12.581517 Ap: 4.865807 tp: 0.069247
Ai: 144.377717 ti: 12.581517 Ap: 4.865807 tp: 0.069203
Ai: 144.377717 ti: 12.581517 Ap: 4.865806 tp: 0.069114
Ai: 144.377717 ti: 12.581517 Ap: 4.865807 tp: 0.069194
Ai: 144.377717 ti: 12.581517 Ap: 4.865806 tp: 0.069176
Ai: 144.377717 ti: 12.581517 Ap: 4.865806 tp: 0.069140
Ai: 144.377717 ti: 12.581517 Ap: 4.865806 tp: 0.069172
Ai: 144.377717 ti: 12.581517 Ap: 4.865806 tp: 0.069165
Ai: 144.377717 ti: 12.581517 Ap: 4.865806 tp: 0.069172
here I got exp: Overflow error
the source data was generated with parameters
Ai = 100
ti = 10
Ap = 50
tp = 20
ofset = 100
and the noise was in range (-5,5);
Interesting thing is that this was the 23rd try, 22 tries before ended
succesfully and their result was very good (when ploted, they didnot differ
much from original curve).
If I add more noise to source data, this behaviour becomes more often. All is in first
iteration step, if solver "listen more" to my advice (initial guess), I am
pretty sure that it would converge.
I add just one more example, I used same example program from GSL manual, I
copied it line by line. I only changed equation this way
f = a*exp(lambda*t)+b; (I just removed minus from expression)
derivations follows:
gsl_matrix_set (m, i, 0, exp(lambda*t));
gsl_matrix_set (m, i, 1, t * a * exp(lambda*t));
gsl_matrix_set (m, i, 2, 1);
parameters used by generator
a = 5;
lambda = 0.5;
b = 50;
no noise at all, starting guess
a = 5;
lambda = 1;
b = 20;
and here is the log. Notice behaviour of b variable
a: 5.000000 lambda: 1.000000 b: 20.000000
a: 0.000000 lambda: 1.000000 b: 1146731183.144235
a: 0.000000 lambda: 0.980099 b: 1146710605.800375
a: 0.000000 lambda: 0.940297 b: 1580849253.130197
a: 0.000000 lambda: 0.860693 b: 2113733309.810215
a: 0.000000 lambda: 0.701485 b: 2649974321.661100
a: 0.000000 lambda: 0.810630 b: 2810395464.666553
a: 0.000000 lambda: 0.841103 b: 2855184146.770248
a: 0.000000 lambda: 0.821512 b: 2728600930.934087
a: 0.000001 lambda: 0.801922 b: 2607082703.329340
a: 0.000003 lambda: 0.792127 b: 2514172932.704934
a: 0.000006 lambda: 0.772536 b: 2402488753.208288
a: 0.000004 lambda: 0.786289 b: 2459147111.154609
a: 0.000007 lambda: 0.774614 b: 2390137152.747778
a: 0.000005 lambda: 0.782880 b: 2430175024.612602
a: 0.000007 lambda: 0.776061 b: 2389927577.762827
a: 0.000006 lambda: 0.780008 b: 2410334359.097036
a: 0.000008 lambda: 0.774265 b: 2376709990.614474
a: 0.000009 lambda: 0.771455 b: 2355576602.765717
a: 0.000011 lambda: 0.765834 b: 2321923850.078345
a: 0.000013 lambda: 0.763079 b: 2300896835.368453
a: 0.000017 lambda: 0.757568 b: 2267314649.661103
a: 0.000020 lambda: 0.754865 b: 2246423317.138189
a: 0.000026 lambda: 0.749460 b: 2212905540.207484
a: 0.000030 lambda: 0.746808 b: 2192145393.455541
a: 0.000038 lambda: 0.741503 b: 2158685016.720283
a: 0.000044 lambda: 0.738900 b: 2138051720.468261
a: 0.000055 lambda: 0.733693 b: 2104642830.750894
a: 0.000071 lambda: 0.728624 b: 2069618974.062552
a: 0.000090 lambda: 0.723609 b: 2034221399.314935
a: 0.000115 lambda: 0.718638 b: 1998724971.209173
a: 0.000147 lambda: 0.713710 b: 1963180931.539362
a: 0.000186 lambda: 0.708827 b: 1927609541.393495
a: 0.000236 lambda: 0.703989 b: 1892027086.189487
a: 0.000298 lambda: 0.699199 b: 1856449669.610041
a: 0.000376 lambda: 0.694458 b: 1820893728.648119
a: 0.000473 lambda: 0.689768 b: 1785376083.901281
a: 0.000607 lambda: 0.684561 b: 1746243314.607982
a: 0.000779 lambda: 0.679450 b: 1706888412.132023
a: 0.000996 lambda: 0.674412 b: 1667534644.748884
a: 0.001271 lambda: 0.669443 b: 1628302230.921109
a: 0.001615 lambda: 0.664546 b: 1589230244.635055
a: 0.002046 lambda: 0.659722 b: 1550346475.142906
a: 0.002582 lambda: 0.654972 b: 1511677285.562059
a: 0.003312 lambda: 0.649765 b: 1469473636.444143
a: 0.004246 lambda: 0.644683 b: 1427352052.505335
a: 0.005420 lambda: 0.639701 b: 1385467284.947267
a: 0.006885 lambda: 0.634819 b: 1343959297.684775
a: 0.008703 lambda: 0.630037 b: 1302878258.316209
a: 0.011170 lambda: 0.624830 b: 1258271455.694746
a: 0.014300 lambda: 0.619778 b: 1214074427.639879
a: 0.018201 lambda: 0.614858 b: 1170390818.316526
a: 0.023023 lambda: 0.610068 b: 1127381059.250725
a: 0.029549 lambda: 0.604861 b: 1080722738.059405
a: 0.037771 lambda: 0.599846 b: 1034855173.940194
a: 0.047924 lambda: 0.594998 b: 989834258.222745
a: 0.061540 lambda: 0.589792 b: 941480416.181589
a: 0.078555 lambda: 0.584814 b: 894338191.065179
a: 0.099352 lambda: 0.580038 b: 848423015.814896
a: 0.127478 lambda: 0.574833 b: 798433000.679436
a: 0.147542 lambda: 0.572298 b: 771314147.830819
a: 0.184154 lambda: 0.567415 b: 724555293.723404
a: 0.164824 lambda: 0.570111 b: 749125718.535294
a: 0.199270 lambda: 0.565943 b: 709258454.498209
a: 0.222672 lambda: 0.563929 b: 688044045.265559
a: 0.265892 lambda: 0.560064 b: 650320642.317852
a: 0.317345 lambda: 0.556431 b: 613843003.410136
a: 0.384690 lambda: 0.552407 b: 573130198.889405
a: 0.462657 lambda: 0.548646 b: 534390674.203019
a: 0.561260 lambda: 0.544634 b: 492876820.972454
a: 0.673436 lambda: 0.540932 b: 453984031.833932
a: 0.815871 lambda: 0.536941 b: 411934487.560946
a: 0.974782 lambda: 0.533338 b: 373364457.261235
a: 1.177969 lambda: 0.529391 b: 331025986.734599
a: 1.304584 lambda: 0.527573 b: 310318099.194578
a: 1.542785 lambda: 0.523921 b: 270929848.860883
a: 1.405788 lambda: 0.526084 b: 293751872.290390
a: 1.592533 lambda: 0.523406 b: 264864557.715904
a: 1.832252 lambda: 0.520474 b: 232884262.509994
a: 2.107358 lambda: 0.517594 b: 201047972.572021
a: 2.417435 lambda: 0.514774 b: 169573613.462840
a: 2.763754 lambda: 0.512026 b: 138655778.468354
a: 3.147084 lambda: 0.509361 b: 108443388.282738
a: 3.567305 lambda: 0.506791 b: 79074721.109322
a: 4.023374 lambda: 0.504325 b: 50671714.970997
a: 4.513280 lambda: 0.501970 b: 23340806.631467
a: 4.979493 lambda: 0.499982 b: 47286.200140
a: 4.999996 lambda: 0.500000 b: 119.237751
a: 4.999998 lambda: 0.500000 b: 115.400448
a: 4.999998 lambda: 0.500000 b: 115.400339
error at check of status after solver iteration
I am really out of ideas, I would appreciate any help
Best regards
Petr Ent
> 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 mailing list
address@hidden
http://lists.gnu.org/mailman/listinfo/help-gsl