C++里实现scipy里的curve_fit

在C++中实现Scipy中的curve_fit函数可以使用非线性最小二乘回归方法来拟合非线性函数。以下是一个简单的实现：

首先，需要定义一个函数去拟合，这个函数可以是任何非线性函数。比如，定义一个指数函数：

double func(double x, double a, double b, double c) {
    return a*exp(-b*x) + c;
}

然后，创建一个vector去存储x和y的值：

std::vector<double> x_data;
std::vector<double> y_data;

将x和y的值分别加入到vector中：

x_data.push_back(x1);
y_data.push_back(y1);
x_data.push_back(x2);
y_data.push_back(y2);
//...

接下来，使用非线性最小二乘回归方法来拟合函数：

#include <gsl/gsl_multifit_nlin.h>

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

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

    double a = gsl_vector_get(x, 0);
    double b = gsl_vector_get(x, 1);
    double c = gsl_vector_get(x, 2);

    for (size_t i = 0; i < n; i++) {
        double t = x_data[i];
        double yi = a * exp(-b * t) + c;
        gsl_vector_set(f, i, yi - y_data[i]);
    }

    return GSL_SUCCESS;
}

int main() {
    const gsl_multifit_fdfsolver_type *T;
    gsl_multifit_fdfsolver *s;
    int status;
    size_t i, iter = 0;

    const size_t n = x_data.size();
    const size_t p = 3;

    gsl_matrix *covar = gsl_matrix_alloc(p, p);
    double x_init[3] = {1.0, 1.0, 1.0};
    gsl_vector_view x = gsl_vector_view_array(x_init, p);

    gsl_multifit_function_fdf f;
    struct data d = {n, x_data.data(), y_data.data()};

    f.f = &func_f;
    f.df = NULL;
    f.fdf = NULL;
    f.n = n;
    f.p = p;
    f.params = &d;

    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);

        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);

    gsl_matrix_fprintf(stdout, covar, "%g");

    gsl_multifit_fdfsolver_free(s);
    gsl_matrix_free(covar);

    return status;
}

这里使用了GSL库中的非线性最小二乘回归方法。首先，需要定义一个结构体去存储x和y的值：

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

然后，定义一个函数去计算残差，即拟合函数的值与实际值之间的差距：

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

    double a = gsl_vector_get(x, 0);
    double b = gsl_vector_get(x, 1);
    double c = gsl_vector_get(x, 2);

    for (size_t i = 0; i < n; i++) {
        double t = x_data[i];
        double yi = a * exp(-b * t) + c;
        gsl_vector_set(f, i, yi - y_data[i]);
    }

    return GSL_SUCCESS;
}

然后，初始化拟合函数的参数，使用gsl_multifit_fdfsolver_type中的lmsder方法进行拟合：

    const gsl_multifit_fdfsolver_type *T;
    gsl_multifit_fdfsolver *s;
    int status;
    size_t i, iter = 0;

    const size_t n = x_data.size();
    const size_t p = 3;

    gsl_matrix *covar = gsl_matrix_alloc(p, p);
    double x_init[3] = {1.0, 1.0, 1.0};
    gsl_vector_view x = gsl_vector_view_array(x_init, p);

    gsl_multifit_function_fdf f;
    struct data d = {n, x_data.data(), y_data.data()};

    f.f = &func_f;
    f.df = NULL;
    f.fdf = NULL;
    f.n = n;
    f.p = p;
    f.params = &d;

    T = gsl_multifit_fdfsolver_lmsder;
    s = gsl_multifit_fdfsolver_alloc(T, n, p);
    gsl_multifit_fdfsolver_set(s, &f, &x.vector);

最后，使用gsl_multifit_test_delta函数判断迭代是否结束，并输出协方差矩阵：

    do {
        iter++;
        status = gsl_multifit_fdfsolver_iterate(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);

    gsl_matrix_fprintf(stdout, covar, "%g");

    gsl_multifit_fdfsolver_free(s);
    gsl_matrix_free(covar);

    return status;

这个函数的返回值即为拟合的状态，若为0则表示拟合成功。

以上就是一个简单的C++实现Scipy中的curve_fit函数的方法。