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

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

我已经成功地获得了一个工作代码，该代码使用C++包装器从拟合分析中估计出正确的参数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和c。但我想找到一个解决方案，这样我就可以向固定参数a输入任何值，而不需要每次修改高斯函数。

好的。以下是基于中链接的代码的答案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语言的原始代码相比，这段代码非常简单。初始化此处存储为向量xs和ys的x-y值对，并执行一个包含4个易于理解的参数的单个函数：要拟合到数据的函数、函数所依赖的拟合参数的初始值、函数必须拟合到的数据的x数值和相应的

1y你的问题是如何保持这个高级接口，但使用它来拟合只有一些参数的函数；"免费"；，也就是说，可以在拟合过程中更改，而其他值必须固定。例如，使用函数可以访问的全局变量可以很容易地实现这一点，但我们讨厌全局变量，并且在没有真正原因的情况下永远不会使用它们。

我建议使用一个众所周知的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);
这几乎与您用于原始包装器的代码相同，但有两个区别：
现在使用对象名(此处：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>
(顺便说一句：IMHO，这种右手边有已知类型的->语法不是最好的语法，但顺其自然)。其思想是迫使编译器在编译时从array的大小中读取拟合参数的数量。
然后你需要在curve_fit_impl的参数列表中进行类似的调整，这几乎就是了
在这里，我花了很多时间试图弄清楚为什么这个代码不起作用。事实证明，它一直有效，秘密是，如果你将一个函数与一些数据相匹配，你最好提供与解决方案相当接近的初始值。这就是为什么使用这个初始化器std::array{0.444, 0.11}而不是原始的{0, 1}，因为后者不会收敛到任何接近正确答案的地方。
我们真的需要使用显式函数对象吗？也许兰达斯可以？是的，他们会的——这会按照预期编译和运行：

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

相关内容

最新更新

热门标签：