過去に、ROOTとPythonで正規分布のフィッティングをした。
phst.hateblo.jp
だが、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 {
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;
}
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;
}
参考文献
非線形方程式の根 - PukiWiki for PBCG Lab