Hello,
I've corrected your code a bit (attached). Also included is a plot of
the data and fitted model. The output I get from running the program is
below. Note that I rescaled the x axis by subtracting the mean of all x
values - this makes the problem better conditioned and easier to
converge. You can add back the mean later to get back to the original
scale of the problem. To solve the problem with the original X values
you'll probably need either a very accurate initial guess vector, or an
analytic Jacobian function.
====
$ ./nlfit > dat
summary from method 'lmsder'
number of iterations: 7
function evaluations: 80
Jacobian evaluations: 0
reason for stopping: small step size
initial |f(x)| = 6393.58
final |f(x)| = 2100.08
chisq/dof = 19092.3
Amp = 946.59642 +/- 65.37942
mu = -2.41481 +/- 0.01003
sigma = 0.12576 +/- 0.01003
Amp = 5620.88963 +/- 67.58526
mu = 0.13784 +/- 0.00163
sigma = 0.11744 +/- 0.00163
Amp = 2608.71157 +/- 65.71170
mu = 1.27723 +/- 0.00361
sigma = 0.12409 +/- 0.00361
status = success
====
On 03/21/2016 09:26 AM, Hinxx wrote:
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
/* expfit.c -- model functions for exponential + background */
struct data {
size_t n;
double * t;
double * y;
};
int
gaussian_f (const gsl_vector * x, void *data, gsl_vector * f)
{
size_t n = ((struct data *)data)->n;
double *t = ((struct data *)data)->t;
double *y = ((struct data *)data)->y;
double A1 = gsl_vector_get(x, 0);
double mu1 = gsl_vector_get(x, 1);
double sigma1 = gsl_vector_get(x, 2);
double A2 = gsl_vector_get(x, 3);
double mu2 = gsl_vector_get(x, 4);
double sigma2 = gsl_vector_get(x, 5);
double A3 = gsl_vector_get(x, 6);
double mu3 = gsl_vector_get(x, 7);
double sigma3 = gsl_vector_get(x, 8);
size_t i;
for (i = 0; i < n; ++i)
{
double ti = t[i];
double yi = y[i];
double term1 = (ti - mu1) / sigma1;
double term2 = (ti - mu2) / sigma2;
double term3 = (ti - mu3) / sigma3;
double fi = 0.0;
fi += A1 * exp(-0.5*term1*term1);
fi += A2 * exp(-0.5*term2*term2);
fi += A3 * exp(-0.5*term3*term3);
gsl_vector_set(f, i, fi - yi);
}
return GSL_SUCCESS;
}
#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_multifit_nlin.h>
#include <gsl/gsl_statistics.h>
#include "expfit.c"
/* number of data points to fit */
#define N 240
#define P 9
static double gaussian_Y[N];
static double gaussian_X[N];
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);
struct data d;
gsl_multifit_function_fdf f;
double x_init[P] = { 870.0, -2.5, 0.1, 5530.0, 0.2, 0.1, 2554.0, 1.3, 0.1 };
gsl_vector_view x = gsl_vector_view_array (x_init, p);
gsl_vector *res_f;
double chi, chi0;
const double xtol = 1e-6;
const double gtol = 1e-6;
const double ftol = 0.0;
load("./800W5-95kv_ROI.txt", gaussian_X, gaussian_Y, N);
/* subtract mean from gaussian_X */
{
double mean = gsl_stats_mean(gaussian_X, 1, N);
for (i = 0; i < N; ++i)
gaussian_X[i] -= mean;
}
d.t = gaussian_X;
d.y = gaussian_Y;
d.n = N;
f.f = &gaussian_f;
f.df = NULL;
f.n = n;
f.p = p;
f.params = &d;
s = gsl_multifit_fdfsolver_alloc (T, n, p);
/* initialize solver with starting point */
gsl_multifit_fdfsolver_set (s, &f, &x.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, "mu = %.5f +/- %.5f\n", FIT(i+1), c*ERR(i+1));
fprintf (stderr, "sigma = %.5f +/- %.5f\n", FIT(i+2), c*ERR(i+2));
}
}
fprintf (stderr, "status = %s\n", gsl_strerror (status));
/* print out data and fitted model */
for (i = 0; i < N; ++i)
{
printf("%f %f %f\n",
gaussian_X[i],
gaussian_Y[i],
gsl_vector_get(res_f, i) + gaussian_Y[i]);
}
gsl_multifit_fdfsolver_free (s);
gsl_matrix_free (covar);
gsl_matrix_free (J);
return 0;
}
