修复非线性最小二乘 GSL 中拟合函数的参数

我正在编写一些代码，这些代码使用 [GNU 科学库 (GSL)][1] 的非线性最小二乘算法进行曲线拟合。

我已经成功地获得了一个工作代码，该代码使用来自的 C++ 包装器从拟合分析中估计了正确的参数https://github.com/Eleobert/gsl-curve-fit/blob/master/example.cpp https://github.com/Eleobert/gsl-curve-fit/blob/master/example.cpp.

现在，我想修复函数的一些参数以使其适合。我想修改该函数，以便我已经可以输入要固定的参数值。

关于如何做有任何想法吗？ 我在这里展示完整的代码。

这是执行非线性最小二乘拟合的代码：

#include <gsl/gsl_vector.h>
#include <gsl/gsl_multifit_nlinear.h>
#include <iostream>
#include <random>
#include <vector>
#include <cassert>
#include <functional>

template <typename F, size_t... Is>
auto gen_tuple_impl(F func, std::index_sequence<Is...> )
{
    return std::make_tuple(func(Is)...);
}

template <size_t N, typename F>
auto gen_tuple(F func)
{
    return gen_tuple_impl(func, std::make_index_sequence<N>{} );
}

auto internal_solve_system(gsl_vector* initial_params, gsl_multifit_nlinear_fdf *fdf,
                           gsl_multifit_nlinear_parameters *params) -> std::vector<double>
{
    // This specifies a trust region method
    const gsl_multifit_nlinear_type *T = gsl_multifit_nlinear_trust;
    const size_t max_iter = 200;
    const double xtol = 1.0e-8;
    const double gtol = 1.0e-8;
    const double ftol = 1.0e-8;

    auto *work = gsl_multifit_nlinear_alloc(T, params, fdf->n, fdf->p);
    int info;

    // initialize solver
    gsl_multifit_nlinear_init(initial_params, fdf, work);
    //iterate until convergence
    gsl_multifit_nlinear_driver(max_iter, xtol, gtol, ftol, nullptr, nullptr, &info, work);

    // result will be stored here
    gsl_vector * y    = gsl_multifit_nlinear_position(work);
    auto result = std::vector<double>(initial_params->size);

    for(int i = 0; i < result.size(); i++)
    {
        result[i] = gsl_vector_get(y, i);
    }

    auto niter = gsl_multifit_nlinear_niter(work);
    auto nfev  = fdf->nevalf;
    auto njev  = fdf->nevaldf;
    auto naev  = fdf->nevalfvv;

    // nfev - number of function evaluations
    // njev - number of Jacobian evaluations
    // naev - number of f_vv evaluations
    //logger::debug("curve fitted after ", niter, " iterations {nfev = ", nfev, "} {njev = ", njev, "} {naev = ", naev, "}");

    gsl_multifit_nlinear_free(work);
    gsl_vector_free(initial_params);
    return result;
}

auto internal_make_gsl_vector_ptr(const std::vector<double>& vec) -> gsl_vector*
{
    auto* result = gsl_vector_alloc(vec.size());
    int i = 0;
    for(const auto e: vec)
    {
        gsl_vector_set(result, i, e);
        i++;
    }
    return result;
}

template<typename C1>
struct fit_data
{
    const std::vector<double>& t;
    const std::vector<double>& y;
    // the actual function to be fitted
    C1 f;
};


template<typename FitData, int n_params>
int internal_f(const gsl_vector* x, void* params, gsl_vector *f)
{
    auto* d  = static_cast<FitData*>(params);
    // Convert the parameter values from gsl_vector (in x) into std::tuple
    auto init_args = [x](int index)
    {
        return gsl_vector_get(x, index);
    };
    auto parameters = gen_tuple<n_params>(init_args);

    // Calculate the error for each...
    for (size_t i = 0; i < d->t.size(); ++i)
    {
        double ti = d->t[i];
        double yi = d->y[i];
        auto func = [ti, &d](auto ...xs)
        {
            // call the actual function to be fitted
            return d->f(ti, xs...);
        };
        auto y = std::apply(func, parameters);
        gsl_vector_set(f, i, yi - y);
    }
    return GSL_SUCCESS;
}

using func_f_type   = int (*) (const gsl_vector*, void*, gsl_vector*);
using func_df_type  = int (*) (const gsl_vector*, void*, gsl_matrix*);
using func_fvv_type = int (*) (const gsl_vector*, const gsl_vector *, void *, gsl_vector *);


auto internal_make_gsl_vector_ptr(const std::vector<double>& vec) -> gsl_vector*;


auto internal_solve_system(gsl_vector* initial_params, gsl_multifit_nlinear_fdf *fdf,
                           gsl_multifit_nlinear_parameters *params) -> std::vector<double>;

template<typename C1>
auto curve_fit_impl(func_f_type f, func_df_type df, func_fvv_type fvv, gsl_vector* initial_params, fit_data<C1>& fd) -> std::vector<double>
{
    assert(fd.t.size() == fd.y.size());

    auto fdf = gsl_multifit_nlinear_fdf();
    auto fdf_params = gsl_multifit_nlinear_default_parameters();

    fdf.f   = f;
    fdf.df  = df;
    fdf.fvv = fvv;
    fdf.n   = fd.t.size();
    fdf.p   = initial_params->size;
    fdf.params = &fd;

    // "This selects the Levenberg-Marquardt algorithm with geodesic acceleration."
    fdf_params.trs = gsl_multifit_nlinear_trs_lmaccel;
    return internal_solve_system(initial_params, &fdf, &fdf_params);
}

template<typename Callable>
auto curve_fit(Callable f, const std::vector<double>& initial_params, const std::vector<double>& x, const std::vector<double>& y) -> std::vector<double>
{
    // We can't pass lambdas without convert to std::function.
    constexpr auto n = 3;
    assert(initial_params.size() == n);

    auto params = internal_make_gsl_vector_ptr(initial_params);
    auto fd = fit_data<Callable>{x, y, f};
    return curve_fit_impl(internal_f<decltype(fd), n>, nullptr, nullptr, params,  fd);
}

// linspace from https://github.com/Eleobert/meth/blob/master/interpolators.hpp
template <typename Container>
auto linspace(typename Container::value_type a, typename Container::value_type b, size_t n)
{
    assert(b > a);
    assert(n > 1);

    Container res(n);
    const auto step = (b - a) / (n - 1);
    auto val = a;
    for(auto& e: res)
    {
        e = val;
        val += step;
    }
    return res;
}

这是我用于拟合的函数：

double gaussian(double x, double a, double b, double c)
{
    const double z = (x - b) / c;
    return a * std::exp(-0.5 * z * z);
}

最后几行创建了观察数据的假数据集（带有一些正态分布的噪声）并测试拟合曲线函数。

int main()
{
    auto device = std::random_device();
    auto gen    = std::mt19937(device());

    auto xs = linspace<std::vector<double>>(0.0, 1.0, 300);
    auto ys = std::vector<double>(xs.size());

    double a = 5.0, b = 0.4, c = 0.15;

    for(size_t i = 0; i < xs.size(); i++)
    {
        auto y =  gaussian(xs[i], a, b, c);
        auto dist  = std::normal_distribution(0.0, 0.1 * y);
        ys[i] = y + dist(gen);
    }

    auto r = curve_fit(gaussian, {1.0, 0.0, 1.0}, xs, ys);

    std::cout << "result: " << r[0] << ' ' << r[1] << ' ' << r[2] << '\n';
    std::cout << "error : " << r[0] - a << ' ' << r[1] - b << ' ' << r[2] - c << '\n';

}

在这种情况下，我想修复其中之一a, b, c参数并估计其余两个。例如，修复a并估计b and c。但我想找到一个解决方案，以便我可以向固定参数输入任何值a，无需每次都修改高斯函数。

好的。这是基于 linkedin 代码的答案http://github.com/Eleobert/gsl-curve-fit/blob/master/example.cpp http://github.com/Eleobert/gsl-curve-fit/blob/master/example.cpp。但是，这不是问题中发布的代码：您应该相应地更新您的问题，以便其他人可以从问题和答案中受益。

因此，基本上，主要问题是 GSL 是一个用纯 C 编写的库，而您使用的是用 C++ 编写的高级包装器，在上述链接中发布。虽然包装器用现代 C++ 写得很好，但它有一个基本问题：它是“僵化的”——它只能用于它设计的问题的子类，而这个子类是所提供功能的相当狭窄的子集由原始的C代码。

让我们尝试对其进行一些改进，并从包装器的使用方式开始：

double gaussian(double x, double a, double b, double c)
{
    const double z = (x - b) / c;
    return a * std::exp(-0.5 * z * z);
}

int main()
{
    auto device = std::random_device();
    auto gen = std::mt19937(device());

    auto xs = linspace<std::vector<double>>(0.0, 1.0, 300);
    auto ys = std::vector<double>(xs.size());

    double a = 5.0, b = 0.4, c = 0.15;

    for (size_t i = 0; i < xs.size(); i++)
    {
        auto y = gaussian(xs[i], a, b, c);
        auto dist = std::normal_distribution(0.0, 0.1 * y);
        ys[i] = y + dist(gen);
    }

    auto result = curve_fit(gaussian, {1.0, 0.0, 1.0}, xs, ys);
   // use result
}

与原始的 C 语言代码相比，这段代码简单得惊人。一是初始化x-y值对，此处存储为向量xs and ys并执行一个带有 4 个易于理解的参数的函数：要拟合数据的函数、函数所依赖的拟合参数的初始值、x值和相应的y函数必须拟合的数据值。

您的问题是如何保留这个高级接口，但将其用于拟合函数，其中只有一些参数是“自由”的，即可以在拟合过程中更改，而其他参数的值必须固定。这可以很容易地实现，例如使用函数可以访问的全局变量，但我们讨厌全局变量，并且在没有真正原因的情况下从不使用它们。

我建议使用众所周知的 C++ 替代方案：函子。看：

struct gaussian_fixed_a
{
    double a;
    gaussian_fixed_a(double a) : a{a} {}
    double operator()(double x, double b, double c) const { return gaussian(x, a, b, c); }
};

This struct/class引入了一种新类型的函数对象。在构造函数中，参数a被传递并存储在一个对象中。然后，有一个函数调用运算符只需要 3 个参数而不是 4 个，替换为a从它的存储值。这个对象可以假装是高斯分布a固定常量和只有其他 3 个参数，x, b, and c，这可能会有所不同。

我们想这样使用它：

    gaussian_fixed_a g(a);
    auto r2 = curve_fit(g, std::array{0.444, 0.11}, xs, ys);

这与您用于原始包装器的代码几乎相同，但有两个差异：

1. 您现在使用对象名称（此处：g) 而不是函数名
2. 您必须将参数数量传递给curve_fit as a 编译时常量，因为它随后在内部使用它来调用由该数字参数化的模板。我通过使用实现它std::array，编译器可以根据需要在编译时推断其大小。另一种选择是使用令人讨厌的模板语法，curve_fit<2>(...

为此，您需要更改接口curve_fit, from

template <typename Callable>
auto curve_fit(Callable f, const std::vector<double>& initial_params, const std::vector<double>& x,
               const std::vector<double>& y) -> std::vector<double>

to

template <typename Callable, auto n>
auto curve_fit(Callable f, const std::array<double, n>& initial_params, const std::vector<double>& x,
               const std::vector<double>& y) -> std::vector<double>

（顺便说一句：这个->右侧带有众所周知类型的语法不是最好的语法，恕我直言，但就这样吧）。这个想法是强制编译器在编译时从参数的大小中读取拟合参数的数量。array.

然后你需要在参数列表中进行类似的调整curve_fit_impl——差不多就是这样了。

在这里，我花了很多时间试图找出为什么这段代码不起作用。事实证明它一直有效，秘密在于，如果您将函数拟合到某些数据，则最好提供相当接近解决方案的初始值。这就是为什么使用这个初始化器std::array{0.444, 0.11}而不是原来的{0, 1}，因为后者不会收敛到任何接近正确答案的结果。

我们真的需要使用显式函数对象吗？也许 lambda 可以吗？是的，他们会的 - 这会按预期编译和运行：

auto r3 = curve_fit([a](double x, double b, double c) { return gaussian(x, a, b, c); }, std::array{0.444, 0.11}, xs, ys);

这是完整的diff原始代码和修改后的代码之间（没有 lambda）：

7a8
> #include <array>
72c73,74
< auto internal_make_gsl_vector_ptr(const std::vector<double>& vec) -> gsl_vector*
---
> template<auto n>
> auto internal_make_gsl_vector_ptr(const std::array<double, n>& vec) -> gsl_vector*
158,159c160,161
< template <typename Callable>
< auto curve_fit(Callable f, const std::vector<double>& initial_params, const std::vector<double>& x,
---
> template <typename Callable, auto n>
> auto curve_fit(Callable f, const std::array<double, n>& initial_params, const std::vector<double>& x,
163,164c165,166
<     constexpr auto n = 3;
<     assert(initial_params.size() == n);
---
>     //    constexpr auto n = 2;
>     //    assert(initial_params.size() == n);
194a197,204
> 
> struct gaussian_fixed_a
> {
>     double a;
>     gaussian_fixed_a(double a) : a{a} {}
>     double operator()(double x, double b, double c) const { return gaussian(x, a, b, c); }
> };
> 
212c222,224
<     auto r = curve_fit(gaussian, {1.0, 0.0, 1.0}, xs, ys);
---
>     auto r = curve_fit(gaussian, std::array{1.0, 0.0, 1.0}, xs, ys);
>     gaussian_fixed_a g(a);
>     auto r2 = curve_fit(g, std::array{0.444, 0.11}, xs, ys);
215a228,230
>     std::cout << "\n";
>     std::cout << "result: " << r2[0] << ' ' << r2[1] << '\n';
>     std::cout << "error : " << r2[0] - b << ' ' << r2[1] - c << '\n';