Hi all,
I'm new to GSL and numerical methods so please bear with me. I'm using
GSL-2.1.
I'm trying to fit a 3-peak gaussian curve onto experiment data that I
have (800W5-95kv_ROI.txt).
I started from expfit.c and nlfit.c and modified the code to handle
gauss curve.
For my initial tests I also used code/data from test_gaussian.c. I've
played around with the test set and managed to get proper amplitude,
sigma and peak location(s), also when extending test set to hold three
peaks (copies of first) and code to handle three peaks.
When using the real data of some 240 points, I'm seeing weird behavior.
Solver bails out after 1st iteration and I get my initial guessed
parameters back, unchanged.
If I fiddle with the data set a bit to make all the points between the
peaks near zero I get better results, bat still good enough.
Maybe it is worth noting that I also have some sample Matlab code that
manages to find the fit on the presented data set. It is using code from
http://www2.imm.dtu.dk/projects/immoptibox/.
I'm attaching the modified GSL code and real data set if someone might
be able to help out.
Thanks in advance,
Hinko
646.01843 -41
646.05176 -21
646.08502 0
646.11835 -12
646.15167 -24
646.185 1
646.21826 -22
646.25159 -15
646.28491 -21
646.31818 -10
646.3515 -15
646.38483 -7
646.41809 -12
646.45142 -11
646.48474 -15
646.51801 3
646.55133 -6
646.58459 5
646.61792 -5
646.65118 -17
646.68451 -15
646.71777 -20
646.7511 -21
646.78436 8
646.81769 -5
646.85095 13
646.88428 -1
646.91754 13
646.95081 -2
646.98413 -17
647.0174 -10
647.05066 -4
647.08398 -12
647.11725 -13
647.15051 5
647.18384 -16
647.2171 6
647.25037 10
647.28363 52
647.31696 105
647.35022 149
647.38348 297
647.41675 421
647.45001 589
647.48328 724
647.5166 907
647.54987 902
647.58313 927
647.61639 906
647.64966 768
647.68292 655
647.71619 515
647.74945 367
647.78271 246
647.81598 152
647.84924 125
647.88251 68
647.91577 30
647.94904 36
647.9823 20
648.01556 -4
648.04883 17
648.08203 32
648.1153 23
648.14856 63
648.18182 60
648.21509 78
648.24835 70
648.28156 60
648.31482 45
648.34808 46
648.38135 35
648.41455 59
648.44781 32
648.48108 49
648.51428 39
648.54755 11
648.58081 24
648.61401 28
648.64728 30
648.68048 36
648.71375 18
648.74695 19
648.78021 11
648.81348 23
648.84668 15
648.87994 19
648.91315 21
648.94635 9
648.97961 21
649.01282 -8
649.04608 8
649.07928 -13
649.11255 3
649.14575 30
649.17896 10
649.21222 20
649.24542 -16
649.27863 10
649.31189 15
649.34509 -2
649.3783 17
649.4115 12
649.44476 20
649.47797 5
649.51117 16
649.54437 6
649.57758 5
649.61084 11
649.64404 9
649.67725 46
649.71045 27
649.74365 30
649.77686 57
649.81006 111
649.84326 197
649.87646 433
649.90967 824
649.94287 1459
649.97607 2391
650.00928 3431
650.04248 4442
650.07568 5232
650.10889 5603
650.14209 5522
650.17529 5059
650.2085 4381
650.24164 3440
650.27484 2586
650.30804 1777
650.34125 1141
650.37445 737
650.40765 538
650.4408 407
650.474 404
650.5072 370
650.54034 395
650.57355 368
650.60675 397
650.63995 386
650.6731 372
650.7063 371
650.73944 364
650.77264 324
650.80585 316
650.83899 328
650.87219 307
650.90533 279
650.93854 291
650.97168 277
651.00488 303
651.03802 395
651.07123 602
651.10437 944
651.13757 1439
651.17072 1975
651.20392 2385
651.23706 2657
651.2702 2627
651.30341 2486
651.33655 2193
651.36969 1763
651.40289 1419
651.43604 971
651.46918 668
651.50232 472
651.53552 344
651.56866 330
651.60181 334
651.63495 302
651.66815 297
651.70129 295
651.73444 270
651.76758 277
651.80072 270
651.83386 253
651.867 247
651.90015 236
651.93329 232
651.96643 230
651.99957 200
652.03271 204
652.06586 175
652.099 214
652.13214 184
652.16528 184
652.19843 161
652.23157 195
652.26471 164
652.29785 171
652.33099 161
652.36414 159
652.39728 167
652.43036 140
652.4635 152
652.49664 185
652.52979 169
652.56293 149
652.59601 158
652.62915 173
652.66229 132
652.69537 145
652.72852 165
652.76166 146
652.79474 139
652.82788 132
652.86102 176
652.8941 138
652.92725 148
652.96033 132
652.99347 123
653.02655 153
653.05969 125
653.09277 116
653.12592 106
653.159 130
653.19214 143
653.22522 155
653.25836 148
653.29144 120
653.32458 141
653.35767 163
653.39075 143
653.42389 124
653.45697 145
653.49005 145
653.52319 114
653.55627 110
653.58936 104
653.62244 130
653.65558 123
653.68866 125
653.72174 119
653.75482 141
653.7879 98
653.82098 111
653.85413 141
653.88721 157
653.92029 165
653.95337 191
/* expfit.c -- model functions for exponential + background */
struct data {
size_t n;
double * y;
};
int
expb_f (const gsl_vector * x, void *data, gsl_vector * f)
{
size_t n = ((struct data *)data)->n;
double *y = ((struct data *)data)->y;
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]);
}
return GSL_SUCCESS;
}
int
gaussian_f1peak (const gsl_vector * x, void *data, gsl_vector * f)
{
size_t n = ((struct data *)data)->n;
double *y = ((struct data *)data)->y;
// amplitude
double x1 = gsl_vector_get(x, 0);
// sigma
double x2 = gsl_vector_get(x, 1);
// peak
double x3 = gsl_vector_get(x, 2);
size_t i;
for (i = 0; i < n; ++i)
{
// double ti = (7.0 - i) / 2.0;
double ti = (x3 - i) / 2.0;
double yi = y[i];
// double term = ti - x3;
double term = ti;
double fi = x1 * exp(-x2*term*term/2.0) - yi;
gsl_vector_set(f, i, fi);
}
return GSL_SUCCESS;
}
int
gaussian_f2peak (const gsl_vector * x, void *data, gsl_vector * f)
{
size_t n = ((struct data *)data)->n;
double *y = ((struct data *)data)->y;
// amplitude, sigma, peak
double x11 = gsl_vector_get(x, 0);
double x12 = gsl_vector_get(x, 1);
double x13 = gsl_vector_get(x, 2);
double x21 = gsl_vector_get(x, 3);
double x22 = gsl_vector_get(x, 4);
double x23 = gsl_vector_get(x, 5);
size_t i;
for (i = 0; i < n; ++i)
{
double ti = (x13 - i) / 2.0;
// double ti = i;
double yi = y[i];
// double term = ti - x3;
double term = ti;
double fi = x11 * exp(-x12*term*term/2.0);
ti = (x23 - i) / 2.0;
term = ti;
fi += x21 * exp(-x22*term*term/2.0);
fi -= yi;
gsl_vector_set(f, i, fi);
}
return GSL_SUCCESS;
}
int
gaussian_f (const gsl_vector * x, void *data, gsl_vector * f)
{
size_t n = ((struct data *)data)->n;
double *y = ((struct data *)data)->y;
// amplitude, sigma, peak
double x11 = gsl_vector_get(x, 0);
double x12 = gsl_vector_get(x, 1);
double x13 = gsl_vector_get(x, 2);
double x21 = gsl_vector_get(x, 3);
double x22 = gsl_vector_get(x, 4);
double x23 = gsl_vector_get(x, 5);
double x31 = gsl_vector_get(x, 6);
double x32 = gsl_vector_get(x, 7);
double x33 = gsl_vector_get(x, 8);
size_t i;
for (i = 0; i < n; ++i)
{
double ti = (x13 - i) / 2.0;
// double ti = i;
double yi = y[i];
// double term = ti - x3;
double term = ti;
double fi = x11 * exp(-x12*term*term/2.0);
ti = (x23 - i) / 2.0;
term = ti;
fi += x21 * exp(-x22*term*term/2.0);
ti = (x33 - i) / 2.0;
term = ti;
fi += x31 * exp(-x32*term*term/2.0);
fi -= yi;
gsl_vector_set(f, i, fi);
}
return GSL_SUCCESS;
}
int
expb_df (const gsl_vector * x, void *data,
gsl_matrix * J)
{
size_t n = ((struct data *)data)->n;
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 e = exp(-lambda * t);
gsl_matrix_set (J, i, 0, e);
gsl_matrix_set (J, i, 1, -t * A * e);
gsl_matrix_set (J, i, 2, 1.0);
}
return GSL_SUCCESS;
}
all:
$(CC) -c -Wall -O0 -ggdb3 -I${HOME}/workspace/gsl-2.1/build/include
nlfit.c -o nlfit.o
$(CC) nlfit.o -L${HOME}/workspace/gsl-2.1/build/lib -lgsl -lgslcblas
-lm -o nlfit
clean:
rm nlfit *.o
#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>
#include "expfit.c"
/* number of data points to fit */
//#define N 40
#define N 240
#define P 9
//static double gaussian_Y[N] = {
//0.0009, 0.0044, 0.0175, 0.0540, 0.1295, 0.2420, 0.3521, 0.3989,
//0.3521, 0.2420, 0.1295, 0.0540, 0.0175, 0.0044, 0.0009,
//0.0001, 0.0001, 0.0001,
//0.0009, 0.0044, 0.0175, 0.0540, 0.1295, 0.2420, 0.3521, 0.3989,
//0.3521, 0.2420, 0.1295, 0.0540, 0.0175, 0.0044, 0.0009,
//0.0001, 0.0001, 0.0001,
//0.0009, 0.0044, 0.0175, 0.0540, 0.1295, 0.2420, 0.3521, 0.3989,
//0.3521, 0.2420, 0.1295, 0.0540, 0.0175, 0.0044, 0.0009
//};
static double gaussian_Y[N] = {0};
static double gaussian_X[N] = {0};
int load(char *fname, double *x, double *y, int npoints) {
FILE *fp;
int i;
/* Save data for simple plotting with gnuplot */
fp = fopen(fname, "r");
if (! fp) {
perror("fopen(): ");
return -1;
}
i = 0;
while(!feof(fp)) {
if (i >= npoints) {
break;
}
fscanf(fp, "%lf\t%lf\n", &x[i], &y[i]);
i++;
}
fclose(fp);
printf("read %d / %d points\n", i, npoints);
return 0;
}
int
main (void)
{
const gsl_multifit_fdfsolver_type *T = gsl_multifit_fdfsolver_lmsder;
gsl_multifit_fdfsolver *s;
int status, info;
size_t i;
const size_t n = N;
const size_t p = P;
gsl_matrix *J = gsl_matrix_alloc(n, p);
gsl_matrix *covar = gsl_matrix_alloc (p, p);
double y[N], weights[N];
struct data d = { n, y };
gsl_multifit_function_fdf f;
// double x_init[3] = { 1.0, 0.0, 0.0 };
// amplitude, sigma, peak
// double x_init[3] = { 0.4, 1.0, 0.0 };
// double x_init[P] = { 0.4, 1.0, 7.0, 0.4, 1.0, 22.0 };
// double x_init[P] = { 0.4, 1.0, 7.0, 0.4, 1.0, 22.0, 0.4, 1.0, 37.0 };
// double x_init[P] = { 0.8, 1.0, 7.0, 1.2, 1.0, 25.0, 0.4, 1.0, 40.0 };
double x_init[P] = { 870.0, 5.0, 49.0, 5530.0, 5.0, 125.0, 2554.0, 5.0, 160.0 };
gsl_vector_view x = gsl_vector_view_array (x_init, p);
gsl_vector_view w = gsl_vector_view_array(weights, n);
const gsl_rng_type * type;
gsl_rng * r;
gsl_vector *res_f;
double chi, chi0;
const double xtol = 1e-20;
const double gtol = 1e-20;
const double ftol = 0.0;
load("./800W5-95kv_ROI.txt", gaussian_X, gaussian_Y, N);
gsl_rng_env_setup();
type = gsl_rng_default;
r = gsl_rng_alloc (type);
// f.f = &expb_f;
// f.df = &expb_df; /* set to NULL for finite-difference Jacobian */
// for gaussian
f.f = &gaussian_f;
f.df = NULL;
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;
//double yi = 1.0 + 5 * exp (-0.1 * t);
//double si = 0.1 * yi;
// double si = 0.1 * gaussian_Y[i];
//double dy = gsl_ran_gaussian(r, si);
// weights[i] = 1.0 / (si * si);
//y[i] = yi + dy;
// y[i] = gaussian_Y[i];
double yi = gaussian_Y[i];
// if (i < 15) {
// yi *= 2;
// } else if (i < 30) {
// yi *= 3;
// } else {
// yi *= 1;
// }
// if (yi < 500) {
// yi = 0.0009;
// }
y[i] = yi;
double si = 0.1 * yi;
weights[i] = 1.0 / (si * si);
printf ("data: %zu %g %g\n", i, y[i], si);
};
s = gsl_multifit_fdfsolver_alloc (T, n, p);
/* initialize solver with starting point and weights */
gsl_multifit_fdfsolver_wset (s, &f, &x.vector, &w.vector);
/* compute initial residual norm */
res_f = gsl_multifit_fdfsolver_residual(s);
chi0 = gsl_blas_dnrm2(res_f);
/* solve the system with a maximum of 2000 iterations */
status = gsl_multifit_fdfsolver_driver(s, 2000, xtol, gtol, ftol, &info);
gsl_multifit_fdfsolver_jac(s, J);
gsl_multifit_covar (J, 0.0, covar);
/* compute final residual norm */
chi = gsl_blas_dnrm2(res_f);
#define FIT(i) gsl_vector_get(s->x, i)
#define ERR(i) sqrt(gsl_matrix_get(covar,i,i))
fprintf(stderr, "summary from method '%s'\n",
gsl_multifit_fdfsolver_name(s));
fprintf(stderr, "number of iterations: %zu\n",
gsl_multifit_fdfsolver_niter(s));
fprintf(stderr, "function evaluations: %zu\n", f.nevalf);
fprintf(stderr, "Jacobian evaluations: %zu\n", f.nevaldf);
fprintf(stderr, "reason for stopping: %s\n",
(info == 1) ? "small step size" : "small gradient");
fprintf(stderr, "initial |f(x)| = %g\n", chi0);
fprintf(stderr, "final |f(x)| = %g\n", chi);
{
double dof = n - p;
double c = GSL_MAX_DBL(1, chi / sqrt(dof));
fprintf(stderr, "chisq/dof = %g\n", pow(chi, 2.0) / dof);
for (i = 0; i < p; i +=3) {
fprintf (stderr, "Amp = %.5f +/- %.5f\n", FIT(i), c*ERR(i));
fprintf (stderr, "sigma = %.5f +/- %.5f\n", FIT(i+1), c*ERR(i+1));
fprintf (stderr, "peak = %.5f +/- %.5f\n", FIT(i+2), c*ERR(i+2));
}
}
fprintf (stderr, "status = %s\n", gsl_strerror (status));
gsl_multifit_fdfsolver_free (s);
gsl_matrix_free (covar);
gsl_matrix_free (J);
gsl_rng_free (r);
return 0;
}
~/workspace/gslfit$ ./nlfit
read 240 / 240 points
data: 0 -41 -4.1
data: 1 -21 -2.1
data: 2 0 0
data: 3 -12 -1.2
data: 4 -24 -2.4
data: 5 1 0.1
data: 6 -22 -2.2
data: 7 -15 -1.5
data: 8 -21 -2.1
data: 9 -10 -1
data: 10 -15 -1.5
data: 11 -7 -0.7
data: 12 -12 -1.2
data: 13 -11 -1.1
data: 14 -15 -1.5
data: 15 3 0.3
data: 16 -6 -0.6
data: 17 5 0.5
data: 18 -5 -0.5
data: 19 -17 -1.7
data: 20 -15 -1.5
data: 21 -20 -2
data: 22 -21 -2.1
data: 23 8 0.8
data: 24 -5 -0.5
data: 25 13 1.3
data: 26 -1 -0.1
data: 27 13 1.3
data: 28 -2 -0.2
data: 29 -17 -1.7
data: 30 -10 -1
data: 31 -4 -0.4
data: 32 -12 -1.2
data: 33 -13 -1.3
data: 34 5 0.5
data: 35 -16 -1.6
data: 36 6 0.6
data: 37 10 1
data: 38 52 5.2
data: 39 105 10.5
data: 40 149 14.9
data: 41 297 29.7
data: 42 421 42.1
data: 43 589 58.9
data: 44 724 72.4
data: 45 907 90.7
data: 46 902 90.2
data: 47 927 92.7
data: 48 906 90.6
data: 49 768 76.8
data: 50 655 65.5
data: 51 515 51.5
data: 52 367 36.7
data: 53 246 24.6
data: 54 152 15.2
data: 55 125 12.5
data: 56 68 6.8
data: 57 30 3
data: 58 36 3.6
data: 59 20 2
data: 60 -4 -0.4
data: 61 17 1.7
data: 62 32 3.2
data: 63 23 2.3
data: 64 63 6.3
data: 65 60 6
data: 66 78 7.8
data: 67 70 7
data: 68 60 6
data: 69 45 4.5
data: 70 46 4.6
data: 71 35 3.5
data: 72 59 5.9
data: 73 32 3.2
data: 74 49 4.9
data: 75 39 3.9
data: 76 11 1.1
data: 77 24 2.4
data: 78 28 2.8
data: 79 30 3
data: 80 36 3.6
data: 81 18 1.8
data: 82 19 1.9
data: 83 11 1.1
data: 84 23 2.3
data: 85 15 1.5
data: 86 19 1.9
data: 87 21 2.1
data: 88 9 0.9
data: 89 21 2.1
data: 90 -8 -0.8
data: 91 8 0.8
data: 92 -13 -1.3
data: 93 3 0.3
data: 94 30 3
data: 95 10 1
data: 96 20 2
data: 97 -16 -1.6
data: 98 10 1
data: 99 15 1.5
data: 100 -2 -0.2
data: 101 17 1.7
data: 102 12 1.2
data: 103 20 2
data: 104 5 0.5
data: 105 16 1.6
data: 106 6 0.6
data: 107 5 0.5
data: 108 11 1.1
data: 109 9 0.9
data: 110 46 4.6
data: 111 27 2.7
data: 112 30 3
data: 113 57 5.7
data: 114 111 11.1
data: 115 197 19.7
data: 116 433 43.3
data: 117 824 82.4
data: 118 1459 145.9
data: 119 2391 239.1
data: 120 3431 343.1
data: 121 4442 444.2
data: 122 5232 523.2
data: 123 5603 560.3
data: 124 5522 552.2
data: 125 5059 505.9
data: 126 4381 438.1
data: 127 3440 344
data: 128 2586 258.6
data: 129 1777 177.7
data: 130 1141 114.1
data: 131 737 73.7
data: 132 538 53.8
data: 133 407 40.7
data: 134 404 40.4
data: 135 370 37
data: 136 395 39.5
data: 137 368 36.8
data: 138 397 39.7
data: 139 386 38.6
data: 140 372 37.2
data: 141 371 37.1
data: 142 364 36.4
data: 143 324 32.4
data: 144 316 31.6
data: 145 328 32.8
data: 146 307 30.7
data: 147 279 27.9
data: 148 291 29.1
data: 149 277 27.7
data: 150 303 30.3
data: 151 395 39.5
data: 152 602 60.2
data: 153 944 94.4
data: 154 1439 143.9
data: 155 1975 197.5
data: 156 2385 238.5
data: 157 2657 265.7
data: 158 2627 262.7
data: 159 2486 248.6
data: 160 2193 219.3
data: 161 1763 176.3
data: 162 1419 141.9
data: 163 971 97.1
data: 164 668 66.8
data: 165 472 47.2
data: 166 344 34.4
data: 167 330 33
data: 168 334 33.4
data: 169 302 30.2
data: 170 297 29.7
data: 171 295 29.5
data: 172 270 27
data: 173 277 27.7
data: 174 270 27
data: 175 253 25.3
data: 176 247 24.7
data: 177 236 23.6
data: 178 232 23.2
data: 179 230 23
data: 180 200 20
data: 181 204 20.4
data: 182 175 17.5
data: 183 214 21.4
data: 184 184 18.4
data: 185 184 18.4
data: 186 161 16.1
data: 187 195 19.5
data: 188 164 16.4
data: 189 171 17.1
data: 190 161 16.1
data: 191 159 15.9
data: 192 167 16.7
data: 193 140 14
data: 194 152 15.2
data: 195 185 18.5
data: 196 169 16.9
data: 197 149 14.9
data: 198 158 15.8
data: 199 173 17.3
data: 200 132 13.2
data: 201 145 14.5
data: 202 165 16.5
data: 203 146 14.6
data: 204 139 13.9
data: 205 132 13.2
data: 206 176 17.6
data: 207 138 13.8
data: 208 148 14.8
data: 209 132 13.2
data: 210 123 12.3
data: 211 153 15.3
data: 212 125 12.5
data: 213 116 11.6
data: 214 106 10.6
data: 215 130 13
data: 216 143 14.3
data: 217 155 15.5
data: 218 148 14.8
data: 219 120 12
data: 220 141 14.1
data: 221 163 16.3
data: 222 143 14.3
data: 223 124 12.4
data: 224 145 14.5
data: 225 145 14.5
data: 226 114 11.4
data: 227 110 11
data: 228 104 10.4
data: 229 130 13
data: 230 123 12.3
data: 231 125 12.5
data: 232 119 11.9
data: 233 141 14.1
data: 234 98 9.8
data: 235 111 11.1
data: 236 141 14.1
data: 237 157 15.7
data: 238 165 16.5
data: 239 191 19.1
summary from method 'lmsder'
number of iterations: 1
function evaluations: 20
Jacobian evaluations: 0
reason for stopping: small gradient
initial |f(x)| = nan
final |f(x)| = nan
chisq/dof = nan
Amp = 870.00000 +/- nan
sigma = 5.00000 +/- nan
peak = 49.00000 +/- -nan
Amp = 5530.00000 +/- nan
sigma = 5.00000 +/- nan
peak = 125.00000 +/- nan
Amp = 2554.00000 +/- nan
sigma = 5.00000 +/- nan
peak = 160.00000 +/- nan
status = success
