過去に、ROOTとPythonで正規分布のフィッティングをした。
だが、ROOTはWindowsでx64と共存しないので、x64と共存するライブラリだけを使ったフィッティングを試みた。まずは、
ガウス分布のフィッティング用のサンプルデータ Sample data
を、gaussian_sample.txt としてダウンロードしておく。
正規分布へのフィッティングは Eigenライブラリの NonLinearOptimization を使った。
パラメータの誤差を求めるときは、ほぼコピペの二分法を使った。通常、パラメータの誤差はカイ二乗が1.0増えるときであるが、 Constant と Sigma は 1.5増える時がROOTのパラメータと近かったので、そちらを採用した。なぜ1.0ではなく1.5なのか知ってる人はコメントください。
ROOTのパラメータ
Constant 402.531 +/- 4.95816 Mean -0.00397889 +/- 0.00989405 Sigma 0.983115 +/- 0.00704207 Chi2 83.6418
下記のコードで求めたパラメータ
Constant:402.531 +/- 4.9502 Mean :-0.00397891 +/- 0.0102539 Sigma :0.983115 +/- 0.00708008 Chi2 :83.6418
ソースコード
#include <Eigen/Dense> #include <unsupported/Eigen/NonLinearOptimization> #include <unsupported/Eigen/NumericalDiff> #include <iostream> #include <sstream> #include <fstream> namespace { // Generic functor template<typename _Scalar, int NX = Eigen::Dynamic, int NY = Eigen::Dynamic> struct Functor { typedef _Scalar Scalar; enum { InputsAtCompileTime = NX, ValuesAtCompileTime = NY }; typedef Eigen::Matrix<Scalar, InputsAtCompileTime, 1> InputType; typedef Eigen::Matrix<Scalar, ValuesAtCompileTime, 1> ValueType; typedef Eigen::Matrix<Scalar, ValuesAtCompileTime, InputsAtCompileTime> JacobianType; }; double gaussign_function(const Eigen::VectorXd& b, const double& x) { return b[0] * exp(-(x - b[1]) * (x - b[1]) / (2 * b[2] * b[2])); } double gaussign_function(const double& constant, const double& mean, const double& sigma, const double& x) { return constant * exp(-(x - mean) * (x - mean) / (2 * sigma * sigma)); } struct gaussian_functor2 : Functor<double> { gaussian_functor2(int inputs, size_t values, double* x, double* y) : inputs_(inputs), values_(values), x(x), y(y) {} double* x; double* y; int operator()(const Eigen::VectorXd& b, Eigen::VectorXd& fvec) const { for (int i = 0; i < values_; ++i) { if (y[i] == 0) { fvec[i] = 0; } else { double dy = std::abs(gaussign_function(b, x[i]) - y[i]); fvec[i] = dy / std::sqrt(y[i]); } } return 0; } const int inputs_; const int values_; int inputs() const { return inputs_; } int values() const { return values_; } }; double chi2_gaussign_function(const std::vector<double>& x, const std::vector<double>& y, const Eigen::VectorXd& q, double d0 = 0, double d1 = 0, double d2 = 0) { double chi2 = 0; for (int i = 0; i < x.size(); i++) { if (y[i] == 0)continue; double y_exp = gaussign_function(q[0] + d0, q[1] + d1, q[2] + d2, x[i]); chi2 += (y[i] - y_exp) * (y[i] - y_exp) / y[i]; } return chi2; } /// <summary> /// 2分法 /// </summary> /// <param name="func"></param> /// <param name="xr"></param> /// <param name="max_iter"></param> /// <param name="x_found"></param> /// <param name="error"></param> void bisection(std::function<double(const double&)> func, double xr, const int& max_iter, double& x_found, double& error) { double xl = 0; double f = 0; double fmid = 0; for (int k = 0; k < max_iter; k++) { f = func(xl); fmid = func(xr); if (f * fmid < 0.0) break; xr *= 2; } if (f * fmid >= 0.0) { throw std::exception("bisection q1_error"); } double dx = fabs(xr - xl), xmid; for (int k = 0; k < max_iter; ++k) { xmid = 0.5 * (xl + xr); // 中点 dx *= 0.5; // 中点での関数値を求める fmid = func(xmid); // 収束判定 if (dx < error || fmid == 0.0) { x_found = xmid; error = dx; return; } // 新しい区間 if (f * fmid < 0) { xr = xmid; } else { xl = xmid; f = fmid; } } throw std::exception("bisection q2_error"); } Eigen::VectorXd error_gaussign_function(const std::vector<double>& x, const std::vector<double>& y, const Eigen::VectorXd& q, double chi2) { auto func0 = [&x, &y, &q, &chi2](const double& diff) { return chi2_gaussign_function(x, y, q, diff, 0, 0) + chi2_gaussign_function(x, y, q, -diff, 0, 0) - chi2 * 2 - 1.5*2; }; auto func1 = [&x, &y, &q, &chi2](const double& diff) { return chi2_gaussign_function(x, y, q, 0, diff, 0) + chi2_gaussign_function(x, y, q, 0, -diff, 0) - chi2 * 2 - 1.0*2; }; auto func2 = [&x, &y, &q, &chi2](const double& diff) { return chi2_gaussign_function(x, y, q, 0, 0, diff) + chi2_gaussign_function(x, y, q, 0, 0, -diff) - chi2 * 2 - 1.5*2; }; double q0_error, q1_error, q2_error; double bisection_error = 0.001; bisection(func0, 1, 100, q0_error, bisection_error); bisection(func1, 1, 100, q1_error, bisection_error); bisection(func2, 1, 100, q2_error, bisection_error); Eigen::VectorXd q_error(3); q_error << q0_error, q1_error, q2_error; return q_error; } } void main() { std::ifstream ifs("gaussian_sample.txt"); const int n = 3; Eigen::VectorXd q(n); q << 400, 0, 1; std::vector<double> xbuf, ybuf; double x, y; while (ifs >> x >> y) { xbuf.emplace_back(x); ybuf.emplace_back(y); } gaussian_functor2 functor(n, xbuf.size(), xbuf.data(), ybuf.data()); Eigen::NumericalDiff<gaussian_functor2> numDiff(functor); Eigen::LevenbergMarquardt<Eigen::NumericalDiff<gaussian_functor2> > lm(numDiff); int info = lm.minimize(q); double chi2 = chi2_gaussign_function(xbuf, ybuf, q); Eigen::VectorXd error = error_gaussign_function(xbuf, ybuf, q, chi2); std::cout << "Constant:" << q[0] << " +/- " << error[0] << std::endl; std::cout << "Mean :" << q[1] << " +/- " << error[1] << std::endl; std::cout << "Sigma :" << q[2] << " +/- " << error[2] << std::endl; std::cout << "Chi2 :" << chi2 << std::endl; }
参考文献
