c++ - 使用 GSL 的非线性拟合

Question

所以我正在尝试修改我发现的一些代码 here以适应不同的功能，但我稍微修改后的版本无法收敛，我不明白为什么。

我试图找到最小二乘拟合的函数是“A + lambdalog(t) + blog(t)^2。这是代码

主要.cpp

#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"

#define N 40 // All "N"s are 40

void print_state (size_t iter, gsl_multifit_fdfsolver * s);//The prototype of some function

int main (void){

    const gsl_multifit_fdfsolver_type *T; //pointer to a newly allocated instance of a solver of type T
    gsl_multifit_fdfsolver *s;
    int status;
    unsigned int i, iter = 0;
    const size_t n = N; //N number of observations
    const size_t p = 3;//3 parameters

    gsl_matrix *covar = gsl_matrix_alloc (p, p); // creates a pxp gsl_matrix
    double y[N], sigma[N]; //declares vector variables that will hold the noise data. They are "N" long
    struct data d = { n, y, sigma}; //Populates struct d with variables n, y and sigma. Struct data is defined in expfit.c
    gsl_multifit_function_fdf f;
    double x_init[3] = { 0, -.1, -.1 }; //initial x values !These are initial guesses to the solution!
    gsl_vector_view x = gsl_vector_view_array (x_init, p);//view arrays allow one to litterally view elements of a certain array without modifying or created a copy of the original array. Essentially a pointer to the original data.
    const gsl_rng_type * type; // pointer to a new random number generator type. RNG type will be assigned later
    gsl_rng * r; //Pointer to a new RNG

    gsl_rng_env_setup();

    type = gsl_rng_default; //Assigns random number generator type
    r = gsl_rng_alloc (type); // Allocates memory for new RNG of type "type"

    f.f = &logb_f;
    f.df = &logb_df;
    f.fdf = &logb_fdf;
    f.n = n;
    f.p = p;
    f.params = &d;



    for (i = 0; i < n; i++){// This is where the data is being generated

        double t = i; // t is being redclared at each iteration for some reason wtf
        if(t==0){//since log(0) is undefined, I said they equal 0 at t=0
            y[i] = 2.0 -.5 * 0 - 0 + gsl_ran_gaussian (r, 0.1);
        }else{
            y[i] = 2.0 -.5 * log (t) - pow(log(t),2) + gsl_ran_gaussian (r, 0.1); //This is the noised up data
        }
        sigma[i] = .1; //not sure what this sigma does
        printf ("data: %u %g %g\n", i, y[i], sigma[i]);//Printing out the data at each iteration
    };

    T = gsl_multifit_fdfsolver_lmsder; // Not sure what this is doing
    s = gsl_multifit_fdfsolver_alloc (T, n, p);
    gsl_multifit_fdfsolver_set (s, &f, &x.vector);

    print_state (iter, s);

    do{
        iter++;
        status = gsl_multifit_fdfsolver_iterate (s);

        printf ("status = %s\n", gsl_strerror (status));

        print_state (iter, s);

        if (status)
            break;

        status = gsl_multifit_test_delta (s->dx, s->x, 1e-4, 1e-4);

    }while (status == GSL_CONTINUE && iter < 500);

    gsl_multifit_covar (s->J, 0.0, covar);

#define FIT(i) gsl_vector_get(s->x, i)
#define ERR(i) sqrt(gsl_matrix_get(covar,i,i))

    {
        double chi = gsl_blas_dnrm2(s->f);
        double dof = n - p;
        double c = GSL_MAX_DBL(1, chi / sqrt(dof));

        printf("chisq/dof = %g\n",  pow(chi, 2.0) / dof);

        printf ("A      = %.5f +/- %.5f\n", FIT(0), c*ERR(0));
        printf ("lambda = %.5f +/- %.5f\n", FIT(1), c*ERR(1));
        printf ("b      = %.5f +/- %.5f\n", FIT(2), c*ERR(2));
    }

    printf ("status = %s\n", gsl_strerror (status));

    gsl_multifit_fdfsolver_free (s);
    gsl_matrix_free (covar);
    gsl_rng_free (r);
    return 0;
}

void print_state (size_t iter, gsl_multifit_fdfsolver * s){
    printf ("iter: %3zu x = % 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_blas_dnrm2 (s->f));
}

expfit.c

//
//  expfit.c
//  test
//
//  Created by [] on 4/11/15.
//  Copyright (c) 2015 []. All rights reserved.
//
#include <stdlib.h>
#include <stdio.h>
#include <math.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.h"

/* expfit.c -- model functions for exponential + background */

struct data {
    size_t n;
    double * y;
    double * sigma;
};

int logb_f (const gsl_vector * x, void *data, gsl_vector * f){
    size_t n = ((struct data *)data)->n;
    double *y = ((struct data *)data)->y;
    double *sigma = ((struct data *) data)->sigma;

    double A = gsl_vector_get (x, 0);
    double lambda = gsl_vector_get (x, 1);
    double b = gsl_vector_get (x, 2);
    double Yi;//will hold the value of the function to be stored into the vector set
    double t;//time variable.
    size_t i;//iterative variable

    for (i = 0; i < n; i++){
        /* Model Yi = A + lambda*log(i) + b*lambda*log(i)^2 */
        t = i;

        if(t==0){ //need if statement to bypass log(0) when t==0 since the value of log is undefined there
            Yi = A + lambda * log(t) + b * pow(log(t),2);//function for t==0
        }else{
            Yi = A + lambda * log(t) + b * pow(log(t),2); // This is the function we used
        }
        gsl_vector_set (f, i, (Yi - y[i])/sigma[i]);
    }

    return GSL_SUCCESS;
}

int logb_df (const gsl_vector * x, void *data, gsl_matrix * J){
    size_t n = ((struct data *)data)->n;
    double *sigma = ((struct data *) data)->sigma;

    //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++){
        /* Jacobian matrix J(i,j) = dfi / dxj, */
        /* where fi = (Yi - yi)/sigma[i],      */
        /*       Yi = A + lambda*log(i) + b*log(i)^2  */
        /* and the xj are the parameters (A,lambda,b) */
        double t = i;
        double s = sigma[i];

        gsl_matrix_set (J, i, 0, 1/s);
        gsl_matrix_set (J, i, 1, log(t)/s);
        gsl_matrix_set (J, i, 2, pow(log(t),2)/s);
    }
    return GSL_SUCCESS;
}

int logb_fdf (const gsl_vector * x, void *data, gsl_vector * f, gsl_matrix * J){
    logb_f (x, data, f);
    logb_df (x, data, J);

    return GSL_SUCCESS;
}

这是头文件，以备不时之需

//
//  expfit.h
//  test
//
//  Created by [] on 4/11/15.
//  Copyright (c) 2015 []. All rights reserved.
//

#ifndef __test__expfit__
#define __test__expfit__

#include <stdio.h>

#endif /* defined(__test__expfit__) */

Answer 1

计算拟合函数时，考虑特殊情况 t=0 以避免 log(0)，但函数值没有差异:

if(t==0){ 
            Yi = A + lambda * log(t) + b * pow(log(t),2);//function for t==0
        }else{
            Yi = A + lambda * log(t) + b * pow(log(t),2); // This is the function we used
        }

此外，您在计算导数时不会考虑这种情况。因此，我将函数和导数更改如下:

int logb_f (const gsl_vector * x, void *data, gsl_vector * f){
    size_t n = ((struct data *)data)->n;
    double *y = ((struct data *)data)->y;
    double *sigma = ((struct data *) data)->sigma;

    double A = gsl_vector_get (x, 0);
    double lambda = gsl_vector_get (x, 1);
    double b = gsl_vector_get (x, 2);
    double Yi;//will hold the value of the function to be stored into the vector set
    double t;//time variable.
    size_t i;//iterative variable


    for (i = 0; i < n; i++){
        /* Model Yi = A + lambda*log(i) + b*lambda*log(i)^2 */
        t = i;
    if(t==0){ //need if statement to bypass log(0) when t==0 since the value of log is undefined there
        Yi = A ;
        }else{
            Yi = A + lambda * log(t) + b * pow(log(t),2); // This is the function we used
        }

    //Yi = A + lambda * log(t) + b * pow(log(t),2); // This is the function we used
        gsl_vector_set (f, i, (Yi - y[i])/sigma[i]);
    }
    return GSL_SUCCESS;
}

int logb_df (const gsl_vector * x, void *data, gsl_matrix * J){
    size_t n = ((struct data *)data)->n;
    double *sigma = ((struct data *) data)->sigma;

    //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++){
        /* Jacobian matrix J(i,j) = dfi / dxj, */
        /* where fi = (Yi - yi)/sigma[i],      */
        /*       Yi = A + lambda*log(i) + b*log(i)^2  */
        /* and the xj are the parameters (A,lambda,b) */
      // d 
        double t = i;
        double s = sigma[i];
        if(i == 0) 
          {
            gsl_matrix_set (J, i, 0, 0);
            gsl_matrix_set (J, i, 1, 0);
            gsl_matrix_set (J, i, 2, 0);
          }
        else 
          {
            gsl_matrix_set (J, i, 0, 1/s);
            gsl_matrix_set (J, i, 1, log(t)/s);
            gsl_matrix_set (J, i, 2, pow(log(t),2)/s);
          }
        //Yi = A + lambda * log(t) + b * pow(log(t),2); // This is the function we used
    }
    return GSL_SUCCESS;
}

在这种情况下迭代收敛:

iter:   0 x =      0.00000000     -0.10000000     -0.10000000 |f(x)| = 470.77
status = success
iter:   1 x =      2.08763815     -0.60282892     -0.97819822 |f(x)| = 5.5047
status = success
iter:   2 x =      2.08763815     -0.60282892     -0.97819822 |f(x)| = 5.5047
chisq/dof = 0.818964
A      = 2.08764 +/- 0.08245
lambda = -0.60283 +/- 0.07702
b      = -0.97820 +/- 0.01722
status = success

但是，我建议验证拟合函数对于 t -> 0 的可微性。如果有疑问，您还可以通过仅考虑 t > 0 的值来限制上述函数中的拟合范围: for (i = 1; i < n; i++) ...而不是 for (i = 0; i < n; i++) ...