物理の駅 by onsanai

Physics Station → PhSt 質問・疑問・間違いの指摘は、コメントに書くか、直接伝えるときっと良いことがあります。主にWindows or Ubuntu用の記事です

C++ Eigenを使った正規分布のフィッティング手法

過去に、ROOTとPython正規分布のフィッティングをした。

phst.hateblo.jp

だが、ROOTはWindowsでx64と共存しないので、x64と共存するライブラリだけを使ったフィッティングを試みた。まずは、

ガウス分布のフィッティング用のサンプルデータ Sample data

を、gaussian_sample.txt としてダウンロードしておく。

正規分布へのフィッティングは Eigenライブラリの NonLinearOptimization を使った。 パラメータの誤差を求めるときは、ほぼコピペの二分法を使った。通常、パラメータの誤差はカイ二乗が1.0増えるときであるが、 ConstantSigma は 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;
}

参考文献

非線形方程式の根 - PukiWiki for PBCG Lab