Nonlinear Least-Squares Fitting Help
Reid Nichol
sigma@xander.safepassage.ca
Wed Jun 12 12:11:00 GMT 2002
I have solved a differential equation with the gsl successfully and now am attempting to fit the plots back to the function to extract the constants used in solving the DE. As far as I can tell I have done everything to the specs required by the library. But, the numbers I get back are nonsense (except k) and it quits after attempting the first iteration.
The code is attached along with the input data. Details follow.
The function I have is
y = Ar*cos(k*x) - Ai*sin(k*x) + Br*cos(k*x) + Bi*sin(k*x)
which would mean that my fitting function is
y = (Ar*cos(k*x) - Ai*sin(k*x) + Br*cos(k*x) + Bi*sin(k*x) - y_i)/sigma_i
I am triing to solve for Ar, Ai, Br, Bi and k, so my Jacobian would be made up of
dy/dAr = cos(k*x)/sigma_i
dy/dAi = -sin(k*x)/sigma_i
dy/dBr = cos(k*x)/sigma_i
dy/dBi = sin(k*x)/sigma_i
dy/dk = (-x*Ar*sin(k*x) - x*Ai*cos(k*x) - x*Br*sin(k*x) + x*Bi*cos(k*x))/sigma_i
I get the numbers from x_init but with a +/- on the order of 10^34 except for k which comes out fine.
I have checked my math multiple times and found no mistakes and I can't find an error in my code although it is obvious that one exists. I have spent alot of time on this and would appreciate any help!
--
If you truly love the memory, you must set it free()!
-------------- next part --------------
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_multifit_nlin.h>
#define N 100
#define FIT(i) gsl_vector_get(s->x, i)
#define ERR(i) sqrt(gsl_matrix_get(covar, i, i))
struct data {
size_t n;
double *y;
double *sigma;
} ;
int print_state (size_t iter, gsl_multifit_fdfsolver * s) {
printf ("iter: %3u x = % 15.8f % 15.8f % 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_vector_get (s->x, 3),
gsl_vector_get (s->x, 4),
gsl_blas_dnrm2 (s->f));
}
int func_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;
size_t i;
double Ar = gsl_vector_get(x, 0);
double Ai = gsl_vector_get(x, 1);
double Br = gsl_vector_get(x, 2);
double Bi = gsl_vector_get(x, 3);
double k = gsl_vector_get(x, 4);
for (i = 0; i < n; i++ ) {
gsl_vector_set(f, i, (Ar*cos(k*i) - Ai*sin(k*i) + Br*cos(k*i) + Bi*sin(k*i) - y[i])/sigma[i]);
}
return GSL_SUCCESS;
}
int func_df(const gsl_vector *x, void *params, gsl_matrix *J) {
size_t n = ((struct data *)params)->n;
double *sigma = ((struct data *)params)->sigma;
size_t i;
double Ar = gsl_vector_get(x, 0);
double Ai = gsl_vector_get(x, 1);
double Br = gsl_vector_get(x, 2);
double Bi = gsl_vector_get(x, 3);
double k = gsl_vector_get(x, 4);
for(i = 0; i < n; i++) {
gsl_matrix_set(J, i, 0, cos(k*i)/sigma[i]);
gsl_matrix_set(J, i, 1, -sin(k*i)/sigma[i]);
gsl_matrix_set(J, i, 2, cos(k*i)/sigma[i]);
gsl_matrix_set(J, i, 3, sin(k*i)/sigma[i]);
gsl_matrix_set(J, i, 4, (-i*Ar*sin(k*i) - i*Ai*cos(k*i) - i*Br*sin(k*i) + i*Bi*cos(k*i))/sigma[i]);
}
return GSL_SUCCESS;
}
int func_fdf(const gsl_vector *x, void *params, gsl_vector *f, gsl_matrix *J) {
func_f(x, params, f);
func_df(x, params, J);
return GSL_SUCCESS;
}
int main() {
FILE *fp; float t1, t2, t3, t4, t5, t6, t7, t8, t9;
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 = 5;
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[5] = {1.0, 2.0, 3.0, 4.0, 0.5}; /* actual values */
gsl_vector_view x = gsl_vector_view_array(x_init, p);
f.f = &func_f;
f.df = &func_df;
f.fdf = &func_fdf;
f.n = n;
f.p = p;
f.params = &d;
fp = fopen("out", "r");
for(i = 0; i < n; i++) {
fscanf(fp, "%f %f %f %f %f %f %f %f %f", &t1,&t2,&t3,&t4,&t5,&t6,&t7,&t8,&t9);
y[i] = (double)t2;
sigma[i] = 0.1;
}
fclose(fp);
T = gsl_multifit_fdfsolver_lmsder;
s = gsl_multifit_fdfsolver_alloc(T, n, p);
gsl_multifit_fdfsolver_set(s, &f, &x.vector);
do {
iter++;
status = gsl_multifit_fdfsolver_iterate(s);
print_state (iter, s);
if(status) { fprintf(stderr, "iterate failed -- %s\n", gsl_strerror(status)); break; }
if ((status = gsl_multifit_test_delta(s->dx, s->x, 1e-4, 1e-4))) { fprintf(stderr, "test_delta failed -- %s\n", gsl_strerror(status)); };
} while(status == GSL_CONTINUE && iter < 500);
gsl_multifit_covar(s->J, 0.0, covar);
gsl_matrix_fprintf(stdout, covar, "%g");
printf("Ar = %.5f +/- %.5f\n", FIT(0), ERR(0));
printf("Ai = %.5f +/- %.5f\n", FIT(1), ERR(1));
printf("Br = %.5f +/- %.5f\n", FIT(2), ERR(2));
printf("Bi = %.5f +/- %.5f\n", FIT(3), ERR(3));
printf("k = %.5f +/- %.5f\n", FIT(4), ERR(4));
printf("status == %s\n", gsl_strerror(status));
gsl_multifit_fdfsolver_free(s);
return EXIT_SUCCESS;
}
-------------- next part --------------
1.00000 4.46918 4.46918 -0.08127 -0.08127 4.30664 4.30664 -2.31586 -2.31586
2.00000 3.84415 3.84415 -1.14264 -1.14264 1.55887 1.55887 -3.06472 -3.06472
3.00000 2.27794 2.27794 -1.92425 -1.92425 -1.57057 -1.57057 -3.06322 -3.06322
4.00000 0.15401 0.15401 -2.23474 -2.23474 -4.31548 -4.31548 -2.31175 -2.31175
5.00000 -2.00763 -2.00763 -1.99809 -1.99809 -6.00381 -6.00381 -0.99427 -0.99427
6.00000 -3.67773 -3.67773 -1.27223 -1.27223 -6.22219 -6.22219 0.56663 0.56663
7.00000 -4.44739 -4.44739 -0.23489 -0.23489 -4.91717 -4.91717 1.98881 1.98881
8.00000 -4.12818 -4.12818 0.85996 0.85996 -2.40826 -2.40826 2.92405 2.92405
9.00000 -2.79824 -2.79824 1.74426 1.74426 0.69029 0.69029 3.14339 3.14339
10.00000 -0.78320 -0.78320 2.20151 2.20151 3.61982 3.61982 2.59311 2.59311
11.00000 1.42360 1.42360 2.11975 2.11975 5.66310 5.66310 1.40795 1.40795
12.00000 3.28185 3.28185 1.51900 1.51900 6.31985 6.31985 -0.12192 -0.12192
13.00000 4.33659 4.33659 0.54635 0.54635 5.42929 5.42929 -1.62195 -1.62195
14.00000 4.32958 4.32958 -0.56007 -0.56007 3.20944 3.20944 -2.72486 -2.72486
15.00000 3.26254 3.26254 -1.52936 -1.52936 0.20381 0.20381 -3.16064 -3.16064
16.00000 1.39672 1.39672 -2.12422 -2.12422 -2.85172 -2.85172 -2.82257 -2.82257
17.00000 -0.81107 -0.81107 -2.19899 -2.19899 -5.20905 -5.20905 -1.79345 -1.79345
18.00000 -2.82028 -2.82028 -1.73537 -1.73537 -6.29102 -6.29102 -0.32523 -0.32523
19.00000 -4.13899 -4.13899 -0.84687 -0.84687 -5.83273 -5.83273 1.22263 1.22263
20.00000 -4.44433 -4.44433 0.24897 0.24897 -3.94639 -3.94639 2.47113 2.47113
21.00000 -3.66154 -3.66154 1.28385 1.28385 -1.09383 -1.09383 3.11462 3.11462
22.00000 -1.98228 -1.98228 2.00441 2.00441 2.02653 2.02653 2.99554 2.99554
23.00000 0.18231 0.18231 2.23421 2.23421 4.65073 4.65073 2.14305 2.14305
24.00000 2.30227 2.30227 1.91700 1.91700 6.13627 6.13627 0.76586 0.76586
25.00000 3.85855 3.85855 1.13044 1.13044 6.11943 6.11943 -0.79883 -0.79883
26.00000 4.47012 4.47012 0.06711 0.06711 4.60435 4.60435 -2.16795 -2.16795
27.00000 3.98725 3.98725 -1.01265 -1.01265 1.96196 1.96196 -3.00627 -3.00627
28.00000 2.52816 2.52816 -1.84448 -1.84448 -1.16079 -1.16079 -3.10856 -3.10856
29.00000 0.45009 0.45009 -2.22471 -2.22471 -3.99934 -3.99934 -2.44976 -2.44976
30.00000 -1.73818 -1.73818 -2.06026 -2.06026 -5.85870 -5.85870 -1.19118 -1.19118
31.00000 -3.50088 -3.50088 -1.39139 -1.39139 -6.28366 -6.28366 0.35905 0.35905
32.00000 -4.40644 -4.40644 -0.38185 -0.38185 -5.17015 -5.17015 1.82137 1.82137
33.00000 -4.23316 -4.23316 0.72117 0.72117 -2.79081 -2.79081 2.83775 2.83775
34.00000 -3.02345 -3.02345 1.64763 1.64763 0.27181 0.27181 3.15936 3.15936
35.00000 -1.07349 -1.07349 2.17069 2.17069 3.26789 3.26789 2.70744 2.70744
36.00000 1.13929 1.13929 2.16229 2.16229 5.46387 5.46387 1.59265 1.59265
37.00000 3.07314 3.07314 1.62449 1.62449 6.32211 6.32211 0.08792 0.08792
38.00000 4.25457 4.25457 0.68895 0.68895 5.63247 5.63247 -1.43834 -1.43834
39.00000 4.39434 4.39434 -0.41526 -0.41526 3.56381 3.56381 -2.61243 -2.61243
40.00000 3.45822 3.45822 -1.41781 -1.41781 0.62260 0.62260 -3.14692 -3.14692
41.00000 1.67541 1.67541 -2.07322 -2.07322 -2.47104 -2.47104 -2.91093 -2.91093
42.00000 -0.51761 -0.51761 -2.22104 -2.22104 -4.95969 -4.95969 -1.96224 -1.96224
43.00000 -2.58389 -2.58389 -1.82507 -1.82507 -6.23403 -6.23403 -0.53313 -0.53313
44.00000 -4.01755 -4.01755 -0.98226 -0.98226 -5.98206 -5.98206 1.02651 1.02651
45.00000 -4.46757 -4.46757 0.10104 0.10104 -4.26548 -4.26548 2.33483 2.33483
46.00000 -3.82377 -3.82377 1.15961 1.15961 -1.50456 -1.50456 3.07149 3.07149
47.00000 -2.24379 -2.24379 1.93426 1.93426 1.62473 1.62473 3.05615 3.05615
48.00000 -0.11444 -0.11444 2.23534 2.23534 4.35623 4.35623 2.29256 2.29256
49.00000 2.04292 2.04292 1.98912 1.98912 6.02117 6.02117 0.96766 0.96766
50.00000 3.70011 3.70011 1.25591 1.25591 6.21192 6.21192 -0.59415 -0.59415
51.00000 4.45138 4.45138 0.21520 0.21520 4.88177 4.88177 -2.01049 -2.01049
52.00000 4.11279 4.11279 -0.87820 -0.87820 2.35640 2.35640 -2.93459 -2.93459
53.00000 2.76726 2.76726 -1.75658 -1.75658 -0.74590 -0.74590 -3.14021 -3.14021
54.00000 0.74420 0.74420 -2.20489 -2.20489 -3.66558 -3.66558 -2.57699 -2.57699
55.00000 -1.46107 -1.46107 -2.11337 -2.11337 -5.68780 -5.68780 -1.38283 -1.38283
56.00000 -3.30861 -3.30861 -1.50442 -1.50442 -6.31745 -6.31745 0.14989 0.14989
57.00000 -4.34609 -4.34609 -0.52713 -0.52713 -5.40036 -5.40036 1.64591 1.64591
58.00000 -4.31950 -4.31950 0.57921 0.57921 -3.16108 -3.16108 2.73896 2.73896
59.00000 -3.23534 -3.23534 1.54374 1.54374 -0.14785 -0.14785 3.16141 3.16141
60.00000 -1.35906 -1.35906 2.13031 2.13031 2.90157 2.90157 2.80984 2.80984
61.00000 0.84997 0.84997 2.19531 2.19531 5.24059 5.24059 1.77033 1.77033
62.00000 2.85089 2.85089 1.72282 1.72282 6.29653 6.29653 0.29737 0.29737
63.00000 4.15382 4.15382 0.82852 0.82852 5.81086 5.81086 -1.24839 -1.24839
64.00000 4.43975 4.43975 -0.26863 -0.26863 3.90249 3.90249 -2.48850 -2.48850
65.00000 3.63867 3.63867 -1.30001 -1.30001 1.03865 1.03865 -3.11934 -3.11934
66.00000 1.94672 1.94672 -2.01310 -2.01310 -2.07948 -2.07948 -2.98646 -2.98646
67.00000 -0.22186 -0.22186 -2.23331 -2.23331 -4.68849 -4.68849 -2.12239 -2.12239
68.00000 -2.33612 -2.33612 -1.90674 -1.90674 -6.14959 -6.14959 -0.73868 -0.73868
69.00000 -3.87841 -3.87841 -1.11332 -1.11332 -6.10505 -6.10505 0.82588 0.82588
70.00000 -4.47113 -4.47113 -0.04733 -0.04733 -4.56579 -4.56579 2.18824 2.18824
71.00000 -3.96917 -3.96917 1.03025 1.03025 -1.90866 -1.90866 3.01484 3.01484
72.00000 -2.49541 -2.49541 1.85559 1.85559 1.21578 1.21578 3.10330 3.10330
73.00000 -0.41069 -0.41069 2.22662 2.22662 4.04255 4.04255 2.43197 2.43197
74.00000 1.77458 1.77458 2.05249 2.05249 5.87956 5.87956 1.16520 1.16520
75.00000 3.52537 3.52537 1.37584 1.37584 6.27705 6.27705 -0.38685 -0.38685
76.00000 4.41303 4.41303 0.36234 0.36234 5.13770 5.13770 -1.84418 -1.84418
77.00000 4.22023 4.22023 -0.73988 -0.73988 2.74047 2.74047 -2.84999 -2.84999
78.00000 2.99416 2.99416 -1.66095 -1.66095 -0.32773 -0.32773 -3.15803 -3.15803
79.00000 1.03502 1.03502 -2.17536 -2.17536 -3.31569 -3.31569 -2.69287 -2.69287
80.00000 -1.17753 -1.17753 -2.15716 -2.15716 -5.49185 -5.49185 -1.56840 -1.56840
81.00000 -3.10178 -3.10178 -1.61082 -1.61082 -6.32342 -6.32342 -0.05993 -0.05993
82.00000 -4.26660 -4.26660 -0.67009 -0.67009 -5.60679 -5.60679 1.46321 1.46321
83.00000 -4.38682 -4.38682 0.43470 0.43470 -3.51742 -3.51742 2.62810 2.62810
84.00000 -3.43298 -3.43298 1.43306 1.43306 -0.56687 -0.56687 3.14955 3.14955
85.00000 -1.63864 -1.63864 2.08056 2.08056 2.52248 2.52248 2.89988 2.89988
86.00000 0.55690 0.55690 2.21866 2.21866 4.99423 4.99423 1.94021 1.94021
87.00000 2.61610 2.61610 1.81356 1.81356 6.24322 6.24322 0.50551 0.50551
88.00000 4.03478 4.03478 0.96444 0.96444 5.96366 5.96366 -1.05295 -1.05295
89.00000 4.46560 4.46560 -0.12081 -0.12081 4.22398 4.22398 -2.35361 -2.35361
90.00000 3.80310 3.80310 -1.17649 -1.17649 1.45012 1.45012 -3.07803 -3.07803
91.00000 2.20946 2.20946 -1.94411 -1.94411 -1.67877 -1.67877 -3.04884 -3.04884
92.00000 0.07486 0.07486 -2.23575 -2.23575 -4.39664 -4.39664 -2.27319 -2.27319
93.00000 -2.07806 -2.07806 -1.98001 -1.98001 -6.03807 -6.03807 -0.94098 -0.94098
94.00000 -3.72220 -3.72220 -1.23948 -1.23948 -6.20116 -6.20116 0.62162 0.62162
95.00000 -4.45501 -4.45501 -0.19549 -0.19549 -4.84599 -4.84599 2.03202 2.03202
96.00000 -4.09709 -4.09709 0.89636 0.89636 -2.30436 -2.30436 2.94491 2.94491
97.00000 -2.73605 -2.73605 1.76876 1.76876 0.80147 0.80147 3.13678 3.13678
98.00000 -0.70514 -0.70514 2.20810 2.20810 3.71106 3.71106 2.56067 2.56067
99.00000 1.49842 1.49842 2.10682 2.10682 5.71206 5.71206 1.35761 1.35761
More information about the Gsl-discuss
mailing list