GSL非线性最小二乘曲线拟合中固定待拟合函数参数的实现方案问询

核心思路 #

• 将原代码中硬编码的拟合参数数量改为模板参数，适配不同数量的待拟合参数
• 利用C++ lambda表达式对原拟合函数的指定参数做值绑定，生成参数更少的新拟合函数，无需修改原函数定义

具体修改步骤 #

第一步：修改curve_fit函数，支持动态指定待拟合参数数量 #

把原来的curve_fit函数替换为以下代码，移除硬编码的参数数量限制：

template<size_t n_params, 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>
{
    assert(initial_params.size() == n_params);

    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_params>, nullptr, nullptr, params,  fd);
}

第二步：使用lambda绑定固定参数 #

无需修改原高斯函数定义，只需要在调用拟合前用lambda绑定要固定的参数值即可，以下是不同固定场景的示例：

1. 固定a=5.0，拟合b和c（共2个待估参数）

// 绑定第一个参数a为固定值5.0，新函数仅接收b、c两个待拟合参数
auto fixed_a_func = [](double x, double b, double c) {
    return gaussian(x, 5.0, b, c);
};
// 调用时指定模板参数为2，初始参数只传b、c的初始值
auto r = curve_fit<2>(fixed_a_func, {0.0, 1.0}, xs, ys);
// 输出结果r[0]为拟合得到的b，r[1]为拟合得到的c，a固定为5.0

2. 固定b=0.4，拟合a和c

auto fixed_b_func = [](double x, double a, double c) {
    return gaussian(x, a, 0.4, c);
};
auto r = curve_fit<2>(fixed_b_func, {1.0, 1.0}, xs, ys);

3. 固定c=0.15，拟合a和b

auto fixed_c_func = [](double x, double a, double b) {
    return gaussian(x, a, b, 0.15);
};
auto r = curve_fit<2>(fixed_c_func, {1.0, 0.0}, xs, ys);

如果需要固定两个参数、仅拟合1个参数也可以直接扩展，比如固定a和b、只拟合c：

auto fixed_ab_func = [](double x, double c) {
    return gaussian(x, 5.0, 0.4, c);
};
auto r = curve_fit<1>(fixed_ab_func, {1.0}, xs, ys);

完整测试示例 #

修改后的main函数测试代码：

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

    // 测试固定a=5.0，拟合b、c
    auto fixed_a_func = [](double x, double b, double c) {
        return gaussian(x, 5.0, b, c);
    };
    auto r = curve_fit<2>(fixed_a_func, {0.0, 1.0}, xs, ys);

    std::cout << "固定a=5.0，拟合结果：b = " << r[0] << "，c = " << r[1] << '\n';
    std::cout << "误差：b误差 = " << r[0] - b << "，c误差 = " << r[1] - c << '\n';

    return 0;
}

内容的提问来源于stack exchange，提问作者CafféSospeso