Hello, there are some tips in the Troubleshooting section of the
Nonlinear least squares manual:
http://www.gnu.org/software/gsl/doc/html/nls.html#troubleshooting
In particular, your problem appears to be an error in your Jacobian
matrix. I've attached a modified version of your program where I
switched to the finite difference Jacobian calculation. I also lowered
the solution tolerance to get a more accurate solution. You'll see that
the numerical solution now agrees with your model to about 1e-7,
although the computed parameters are not exactly what you started with.
This could likely be improved further when you find the error in your
Jacobian function (which I have not done).
Here are the differences with your original program:
> diff main.c main.c.orig
> 188,189c188,189
> < const double xtol = 1.0e-6;
> < const double gtol = 1.0e-6;
> ---
> > const double xtol = 1.0e-2;
> > const double gtol = 1.0e-2;
> 276d275
> < #if 0
> 279,282d277
> < #else
> < fdf.df = NULL;
> < fdf.fvv = NULL;
> < #endif
> 295c290
> < fdf_params.trs = gsl_multifit_nlinear_trs_dogleg;
> ---
> > fdf_params.trs = gsl_multifit_nlinear_trs_lmaccel;
Note I also switched to the dogleg algorithm - in order to use lmaccel
you need an accurate Jacobian routine and fvv routine.
On 10/16/2017 07:56 PM, Jiawei Zhao Zhao wrote:
> Dear GSL mailing list.
>
>
> I am new to GSL and trying to fit a Non linear function. the function is
> given below.
>
>
> y = A*( ( b / (x-r0) ) ^ m- B* ( b / (x-r0) ) ^ n)
>
>
> I modified an example of non linear fitting (code attached). The function
> successfully complied and first deviation gives correct value.
>
>
> However, when I test with generating data from:
>
>
> A=1, b= 2.5, m=12, n=6, r0=0 and B=1
>
>
> and try to fit it start from
>
>
> A=1, b= 2.4, m=12, n=6, r0=0 and B=1
>
>
> The iteration seems stopped at the first steps with info = 27.
>
>
> Can anyone help me on this?
>
>
> Regards.
>
>
> Jiawei Zhao
>
>
#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_multifit_nlinear.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
struct data
{
double *t;
double *y;
size_t n;
};
double
gaussian(const double a, const double b, const double c, const double t)
{
const double z = (t - b) / c;
return (a * exp(-0.5 * z * z));
}
//const double a = 1; //energy
//const double b = 2.5; //sigema,distance
//const double m = 12; /* m */
//const double n = 6; /* n */
//const double r0 = 0; /* cutoff */
//const double B = 1; // Extra attractive
int
func_f (const gsl_vector * x, void *params, gsl_vector * f)
{
struct data *d = (struct data *) params;
double a = gsl_vector_get(x, 0);
double b = gsl_vector_get(x, 1);
double m = gsl_vector_get(x, 2);
double n = gsl_vector_get(x, 3);
double r0 = gsl_vector_get(x, 4);
double B = gsl_vector_get(x, 5);
size_t i;
for (i = 0; i < d->n; ++i)
{
double ti = d->t[i];
double yi = d->y[i];
double y = a*(pow(b/(ti-r0),m)-B*pow(b/(ti-r0),n));//gaussian(a, b, c, ti);
gsl_vector_set(f, i, yi - y);
}
return GSL_SUCCESS;
}
int
func_df (const gsl_vector * x, void *params, gsl_matrix * J)
{
struct data *d = (struct data *) params;
double a = gsl_vector_get(x, 0);
double b = gsl_vector_get(x, 1);
double m = gsl_vector_get(x, 2);
double n = gsl_vector_get(x, 3);
double r0 = gsl_vector_get(x, 4);
double B = gsl_vector_get(x, 5);
size_t i;
for (i = 0; i < d->n; ++i)
{
double ti = d->t[i];
gsl_matrix_set(J, i, 0, pow((-b/(r0 - ti)),m) - B*pow((-b/(r0 - ti)),n));
gsl_matrix_set(J, i, 1, -a*((m*pow((-b/(r0 - ti)),(m - 1)))/(r0 - ti) - (B*n*pow((-b/(r0 - ti)),(n - 1)))/(r0 - ti)) );
gsl_matrix_set(J, i, 2, a*log(-b/(r0 - ti))*pow((-b/(r0 - ti)),m));
gsl_matrix_set(J, i, 3, -B*a*log(-b/(r0 - ti))*pow((-b/(r0 - ti)),n));
gsl_matrix_set(J, i, 4, a*((b*m*pow((-b/(r0 - ti)),(m - 1)))/(r0 - ti)/(r0 - ti) - (B*b*n*pow((-b/(r0 - ti)),(n - 1)))/(r0 - ti)/(r0 - ti)));
gsl_matrix_set(J, i, 5, -a*pow((-b/(r0 - ti)),n));
}
return GSL_SUCCESS;
}
int
func_fvv (const gsl_vector * x, const gsl_vector * v,
void *params, gsl_vector * fvv)
{
struct data *d = (struct data *) params;
double a = gsl_vector_get(x, 0);
double b = gsl_vector_get(x, 1);
double m = gsl_vector_get(x, 2);
double n = gsl_vector_get(x, 3);
double r0 = gsl_vector_get(x, 4);
double B = gsl_vector_get(x, 5);
double va = gsl_vector_get(v, 0);
double vb = gsl_vector_get(v, 1);
double vm = gsl_vector_get(v, 2);
double vn = gsl_vector_get(v, 3);
double vr0 = gsl_vector_get(v, 4);
double vB = gsl_vector_get(v, 5);
size_t i;
for (i = 0; i < d->n; ++i)
{
double ti = d->t[i];
//double Daa = 0;
double Dab = (B*n*pow((-b/(r0 - ti)),(n - 1)))/(r0 - ti) - (m*pow((-b/(r0 - ti)),(m - 1)))/(r0 - ti);
double Dam = log(-b/(r0 - ti))*pow((-b/(r0 - ti)),m);
double Dan = -B*log(-b/(r0 - ti))*pow((-b/(r0 - ti)),n);
double Dar0 = (b*m*pow((-b/(r0 - ti)),(m - 1)))/(r0 - ti)/(r0 - ti) - (B*b*n*pow((-b/(r0 - ti)),(n - 1)))/(r0 - ti)/(r0 - ti);
double DaB = -pow((-b/(r0 - ti)),n);
double Dbb = a*((m*pow((-b/(r0 - ti)),(m - 2))*(m - 1))/(r0 - ti)/(r0 - ti) - (B*n*pow((-b/(r0 - ti)),(n - 2))*(n - 1))/(r0 - ti)/(r0 - ti));
double Dbm = -a*(pow((-b/(r0 - ti)),(m - 1))/(r0 - ti) + (m*log(-b/(r0 - ti))*pow((-b/(r0 - ti)),(m - 1)))/(r0 - ti));
double Dbn = a*((B*pow((-b/(r0 - ti)),(n - 1)))/(r0 - ti) + (B*n*log(-b/(r0 - ti))*pow((-b/(r0 - ti)),(n - 1)))/(r0 - ti));
double Dbr0 = a*((m*pow((-b/(r0 - ti)),(m - 1)))/(r0 - ti)/(r0 - ti) - (B*n*pow((-b/(r0 - ti)),(n - 1)))/(r0 - ti)/(r0 - ti) -
(b*m*pow((-b/(r0 - ti)),(m - 2))*(m - 1))/(r0 - ti)/(r0 - ti)/(r0 - ti) + (B*b*n*pow((-b/(r0 - ti)),(n - 2))*(n - 1))/(r0 - ti)/(r0 - ti)/(r0 - ti));
double DbB = (a*n*pow((-b/(r0 - ti)),(n - 1)))/(r0 - ti);
double Dmm = a*log(-b/(r0 - ti))*log(-b/(r0 - ti))*pow((-b/(r0 - ti)),m);
//double Dmn = 0;
double Dmr0 = (a*b*m*log(-b/(r0 - ti))*pow((-b/(r0 - ti)),(m - 1)))/(r0 - ti)/(r0 - ti) - (a*pow((-b/(r0 - ti)),m))/(r0 - ti);
//double DmB = 0;
double Dnn = -B*a*log(-b/(r0 - ti))*log(-b/(r0 - ti))*pow((-b/(r0 - ti)),n);
double Dnr0 = (B*a*pow((-b/(r0 - ti)),n))/(r0 - ti) - (B*a*b*n*log(-b/(r0 - ti))*pow((-b/(r0 - ti)),(n - 1)))/(r0 - ti)/(r0 - ti);
double DnB = -a*log(-b/(r0 - ti))*pow((-b/(r0 - ti)),n);
double Dr0r0 = -a*((2*b*m*pow((-b/(r0 - ti)),(m - 1)))/(r0 - ti)/(r0 - ti)/(r0 - ti) -
(b*b*m*pow((-b/(r0 - ti)),(m - 2))*(m - 1))/(r0 - ti)/(r0 - ti)/(r0 - ti)/(r0 - ti) -
(2*B*b*n*pow((-b/(r0 - ti)),(n - 1)))/(r0 - ti)/(r0 - ti)/(r0 - ti)
+ (B*b*b*n*pow((-b/(r0 - ti)),(n - 2))*(n - 1))/(r0 - ti)/(r0 - ti)/(r0 - ti)/(r0 - ti));
double Dr0B = -(a*b*n*pow((-b/(r0 - ti)),(n - 1)))/(r0 - ti)/(r0 - ti);
//double DBB =0
double sum;
sum = 2.0 * (va * vb * Dab + va * vm * Dam + va * vn * Dan + va * vr0 * Dar0 + va * vB * DaB) +
2.0 * (vb * vm * Dbm + vb * vn * Dbn + vb * vr0 * Dbr0 + vb * vB * DbB) +
vb * vb * Dbb +
vm * vm * Dmm +
2.0 * vm * vr0 * Dmr0 +
vn * vn * Dnn +
2.0 * (vn * vr0 * Dnr0 + vn * vB * DnB + vr0 * vB * Dr0B) +
vr0 * vr0 * Dr0r0;
gsl_vector_set(fvv, i, sum);
}
return GSL_SUCCESS;
}
void
callback(const size_t iter, void *params,
const gsl_multifit_nlinear_workspace *w)
{
gsl_vector *f = gsl_multifit_nlinear_residual(w);
gsl_vector *x = gsl_multifit_nlinear_position(w);
double avratio = gsl_multifit_nlinear_avratio(w);
double rcond;
(void) params; /* not used */
/* compute reciprocal condition number of J(x) */
gsl_multifit_nlinear_rcond(&rcond, w);
printf("calculating rcond, rcond=%f\n",rcond );
fprintf(stderr, "iter %2zu: a = %.4f, b = %.4f, m = %.4f ,n = %.4f, r0 = %.4f, B = %.4f, |a|/|v| = %.4f cond(J) = %8.4f, |f(x)| = %.4f\n",
iter,
gsl_vector_get(x, 0),
gsl_vector_get(x, 1),
gsl_vector_get(x, 2),
gsl_vector_get(x, 3),
gsl_vector_get(x, 4),
gsl_vector_get(x, 5),
avratio,
1.0 / rcond,
gsl_blas_dnrm2(f));
}
void
solve_system(gsl_vector *x, gsl_multifit_nlinear_fdf *fdf,
gsl_multifit_nlinear_parameters *params)
{
const gsl_multifit_nlinear_type *T = gsl_multifit_nlinear_trust;
const size_t max_iter = 200;
const double xtol = 1.0e-6;
const double gtol = 1.0e-6;
const double ftol = 1.0e-2;
const size_t n = fdf->n;
const size_t p = fdf->p;
gsl_multifit_nlinear_workspace *work =
gsl_multifit_nlinear_alloc(T, params, n, p);
gsl_vector * f = gsl_multifit_nlinear_residual(work);
gsl_vector * y = gsl_multifit_nlinear_position(work);
int info;
double chisq0, chisq, rcond;
/* initialize solver */
gsl_multifit_nlinear_init(x, fdf, work);
printf("ini done\n");
/* store initial cost */
gsl_blas_ddot(f, f, &chisq0);
printf("store done\n");
/* iterate until convergence */
gsl_multifit_nlinear_driver(max_iter, xtol, gtol, ftol,
callback, NULL, &info, work);
printf("convergenced,info = %d\n" ,info);
/* store final cost */
gsl_blas_ddot(f, f, &chisq);
/* store cond(J(x)) */
gsl_multifit_nlinear_rcond(&rcond, work);
gsl_vector_memcpy(x, y);
/* print summary */
fprintf(stderr, "NITER = %zu\n", gsl_multifit_nlinear_niter(work));
fprintf(stderr, "NFEV = %zu\n", fdf->nevalf);
fprintf(stderr, "NJEV = %zu\n", fdf->nevaldf);
fprintf(stderr, "NAEV = %zu\n", fdf->nevalfvv);
fprintf(stderr, "initial cost = %.12e\n", chisq0);
fprintf(stderr, "final cost = %.12e\n", chisq);
fprintf(stderr, "final x = (%.12e, %.12e, %12e)\n",
gsl_vector_get(x, 0), gsl_vector_get(x, 1), gsl_vector_get(x, 2));
fprintf(stderr, "final cond(J) = %.12e\n", 1.0 / rcond);
gsl_multifit_nlinear_free(work);
}
int
main (void)
{
const size_t n = 300; /* number of data points to fit */
const size_t p = 6; /* number of model parameters */
const double a = 1; //energy
const double b = 2.5; //sigema,distance
const double m = 12; /* m */
const double n2 = 6; /* n */
const double r0 = 0; /* cutoff */
const double B = 1; // Extra attractive
const gsl_rng_type * T = gsl_rng_default;
gsl_vector *f = gsl_vector_alloc(n);
gsl_vector *x = gsl_vector_alloc(p);
gsl_multifit_nlinear_fdf fdf;
gsl_multifit_nlinear_parameters fdf_params =
gsl_multifit_nlinear_default_parameters();
struct data fit_data;
gsl_rng * r;
size_t i;
gsl_rng_env_setup ();
r = gsl_rng_alloc (T);
fit_data.t = malloc(n * sizeof(double));
fit_data.y = malloc(n * sizeof(double));
fit_data.n = n;
/* generate synthetic data with noise */
for (i = 0; i < n; ++i)
{
double t = (double)i / (double) n+1.5;
double y0 = a*(pow(b/(t-r0),m)-B*pow(b/(t-r0),n2));
double dy = gsl_ran_gaussian (r, 0.1 * y0);
fit_data.t[i] = t;
fit_data.y[i] = y0 + dy*0;
}
/* define function to be minimized */
fdf.f = func_f;
#if 0
fdf.df = func_df;
fdf.fvv = func_fvv;
#else
fdf.df = NULL;
fdf.fvv = NULL;
#endif
fdf.n = n;
fdf.p = p;
fdf.params = &fit_data;
/* starting point */
gsl_vector_set(x, 0, 1.0);
gsl_vector_set(x, 1, 2.4);
gsl_vector_set(x, 2, 12.0);
gsl_vector_set(x, 3, 6.0);
gsl_vector_set(x, 4, 0.0);
gsl_vector_set(x, 5, 1.0);
fdf_params.trs = gsl_multifit_nlinear_trs_dogleg;
solve_system(x, &fdf, &fdf_params);
/* print data and model */
{
double A = gsl_vector_get(x, 0);
double B = gsl_vector_get(x, 1);
double M = gsl_vector_get(x, 2);
double N = gsl_vector_get(x, 3);
double R0 = gsl_vector_get(x, 4);
double BB = gsl_vector_get(x, 5);
for (i = 0; i < n; ++i)
{
double ti = fit_data.t[i];
double yi = fit_data.y[i];
double fi = A*(pow(B/(ti-R0),M)-BB*pow(B/(ti-R0),N));
printf("%f %f %f\n", ti, yi, fi);
}
fprintf(stderr, "a = %.4f, b = %.4f, m = %.4f ,n = %.4f, r0 = %.4f, B = %.4f\n",
A,//gsl_vector_get(x, 0),
B,//gsl_vector_get(x, 1),
M,//gsl_vector_get(x, 2),
N,//gsl_vector_get(x, 3),
R0,//gsl_vector_get(x, 4),
BB/*gsl_vector_get(x, 5)*/);
}
gsl_vector_free(f);
gsl_vector_free(x);
gsl_rng_free(r);
return 0;
}
