非线性最小二乘 GSL 中拟合函数的固定参数

Question

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

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

现在，我想修复函数的一些参数以使其适合。而且我想以这样一种方式修改函数，即我已经可以输入要固定的参数值。

知道怎么做吗？我在这里展示了完整的代码。

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

#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 和 c。但我想找到一个解决方案，让我可以在固定参数a中输入任何值，而不需要每次都修改高斯函数。

Answer 1

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

所以，基本上，主要问题是 GSL 是一个用纯 C 编写的库，而您使用的是用 C++ 编写的高级包装器，已在上述 link 中发布。虽然包装器在现代 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 和 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); }
};

本struct/class引入了一种新型的函数对象。在构造函数中，参数 a 被传递并存储在一个对象中。然后，有一个只接受 3 个参数而不是 4 个参数的函数调用运算符，从其存储值中替换 a。这个对象可以伪装成具有 a 固定常数和只有其他 3 个参数 x、b 和 c 的高斯分布。

我们想这样使用它：

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

这与您用于原始包装器的代码几乎相同，但有 2 处不同：

您现在使用对象名称（此处：g）而不是函数名称
你必须将参数的数量作为编译时常量传递给curve_fit，因为它随后在内部使用它来调用参数化的模板这个号码。我通过使用 std::array 来实现它，编译器可以根据需要在编译时推断出它的大小。另一种方法是使用讨厌的模板语法，curve_fit<2>(...

为此，您需要更改 curve_fit 的界面，从

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>

至

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

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

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

非线性最小二乘 GSL 中拟合函数的固定参数

Fixing parameters of a fitting function in Nonlinear Least-Square GSL

c

curve-fitting

gsl