

C++の数値計算ライブラリEigenに含まれる、非線形最小二乗法(Levenberg-Marquardt法)で関数フィッティングを実行してみました  
開発環境はVisual Studio 2015です  
フィッティング対象の関数は、y = a×cos(x)^b で、ノイズを含む観測データx,yからパラメータa,bを求めます

[code]

    #include <iostream>
    #include <vector>
    #include <fstream>
    #include <string>
    #include <random>
    #define _USE_MATH_DEFINES
    #include <cmath>
    
    #include <Eigen/Dense>
    #include <unsupported/Eigen/NonLinearOptimization>
    #include <unsupported/Eigen/NumericalDiff>
    
    using namespace std;
    using namespace Eigen;
    
    template<typename _Scalar, int NX = Dynamic, int NY = Dynamic>
    struct Functor
    {
        typedef _Scalar Scalar;
        enum {
            InputsAtCompileTime = NX,
            ValuesAtCompileTime = NY
        };
        typedef Matrix<Scalar, InputsAtCompileTime, 1> InputType;
        typedef Matrix<Scalar, ValuesAtCompileTime, 1> ValueType;
        typedef Matrix<Scalar, ValuesAtCompileTime, InputsAtCompileTime> JacobianType;
    };
    
    struct myFunctor : Functor<double> {
        myFunctor(int inputs, int values, double*x, double*y)
            : inputs_(inputs), values_(values), x(x), y(y) {}
    
        double *x, *y;
    
        int operator()(const VectorXd &b, VectorXd &fvec) const {
            // 目的関数
            for (int i = 0; i < values_; i++) {
                fvec[i] = (b[0] * pow(cos(x[i]), b[1])) - y[i];
            }
            return 0;
        }
    
        const int inputs_;
        const int values_;
        int inputs() const { return inputs_; }
        int values() const { return values_; }
    };
    
    void showLMStatus(int);
    void writeDatas(vector<double> &x, vector<double> &y);
    
    int main()
    {
        const int seed = 20;    // 乱数生成シード値
        const int dataSize = 60;    // データ数
    
        vector<double> x(dataSize), y(dataSize);
    
        const double ans[] = { 0.2, 100 };  // 真のパラメータ
    
        mt19937 mt(seed);
        // xデータをランダムに生成
        uniform_real_distribution<double> x_rand(-M_PI / 2.0, M_PI / 2.0);
        // x,yそれぞれの(正規分布に従う)ノイズを生成
        normal_distribution<double> x_noise(0.0, M_PI / 90.0);
        normal_distribution<double> y_noise(0.0, 0.01);
    
        for (int i = 0; i < dataSize; i++) {
            // 真値
            x[i] = x_rand(mt);
            y[i] = ans[0] * pow(cos(x[i]), ans[1]);
    
            // ノイズを加える
            x[i] += x_noise(mt);
            y[i] += y_noise(mt);
        }
    
        // データの保存
        writeDatas(x, y);
    
        VectorXd p(2);
        p(0) = 0.6;
        p(1) = 30;
    
        myFunctor functor(2, dataSize, &x[0], &y[0]);
    
        NumericalDiff<myFunctor> numDiff(functor);
        LevenbergMarquardt<NumericalDiff<myFunctor> > lm(numDiff);
        lm.parameters.maxfev = (dataSize + 1) * 100;    // max function evaluation
    
        int info = lm.minimize(p);
        showLMStatus(info);
    
        cout << p[0] << ", " << p[1] << endl;
    }
    
    void showLMStatus(int info)
    {
        string status;
        switch (info) {
        case 0:
            status = "ImproperInputParameters";
            break;
        case 1:
            status = "RelativeReductionTooSmall";
            break;
        case 2:
            status = "RelativeErrorTooSmall";
            break;
        case 3:
            status = "RelativeErrorAndReductionTooSmall";
            break;
        case 4:
            status = "CosinusTooSmall";
            break;
        case 5:
            status = "TooManyFunctionEvaluation";
            break;
        case 6:
            status = "FtolTooSmall";
            break;
        case 7:
            status = "XtolTooSmall";
        case 8:
            status = "GtolTooSmall";
            break;
        case 9:
            status = "UserAsked";
            break;
        default:
            status = "Unknown";
            break;
        }
    
        cout << status << endl;
    }
    
    void writeDatas(vector<double> &x, vector<double> &y) {
        ofstream outputFile("./output1.txt");
    
        for (int i = 0; i < x.size(); i++) {
            outputFile << x[i] << "\t" << y[i] << endl;
        }
    
        outputFile.close();
    }
    
[/code]

実行結果は以下のようになりました

[code]

    TooManyFunctionEvaluation
    0.6, 30
    
[/code]

収束しておらずパラメータも初期値のままなので、うまくいっていない様子です  
しかし、Gnuplot(version
4.6)のfitコマンドを用いて同じデータ(保存した"output1.txt")でフィッティングを行ったところ、収束をし、パラメータも真値に近づきました

[code]

    gnuplot> f(x) = a * cos(x)**b
    gnuplot> a = 0.6
    gnuplot> b = 30
    gnuplot> fit f(x) "output1.txt" via a, b
    
    
    Iteration 0
    WSSR        : 1.32176           delta(WSSR) / WSSR : 0
    delta(WSSR) : 0                 limit for stopping : 1e-005
    lambda : 0.217494
    
    initial set of free parameter values
    
    a = 0.6
    b = 30
    /
    
    Iteration 1
    WSSR        : 0.0219191         delta(WSSR) / WSSR : -59.302
    delta(WSSR) : -1.29984          limit for stopping : 1e-005
    lambda : 0.0217494
    
    resultant parameter values
    
    a = 0.125571
    b = 30.0348
    /
    
    Iteration 2
    WSSR        : 0.0214163         delta(WSSR) / WSSR : -0.0234765
    delta(WSSR) : -0.000502779      limit for stopping : 1e-005
    lambda : 0.00217494
    
    resultant parameter values
    
    a = 0.122193
    b = 30.7021
    /
    
    Iteration 3
    WSSR        : 0.0124169         delta(WSSR) / WSSR : -0.724764
    delta(WSSR) : -0.00899934       limit for stopping : 1e-005
    lambda : 0.000217494
    
    resultant parameter values
    
    a = 0.146067
    b = 54.2227
    /
    
    Iteration 4
    WSSR        : 0.00728884        delta(WSSR) / WSSR : -0.703554
    delta(WSSR) : -0.00512809       limit for stopping : 1e-005
    lambda : 2.17494e-005
    
    resultant parameter values
    
    a = 0.185555
    b = 95.4605
    /
    
    Iteration 5
    WSSR        : 0.00708739        delta(WSSR) / WSSR : -0.028424
    delta(WSSR) : -0.000201452      limit for stopping : 1e-005
    lambda : 2.17494e-006
    
    resultant parameter values
    
    a = 0.196069
    b = 106.342
    /
    
    Iteration 6
    WSSR        : 0.00708611        delta(WSSR) / WSSR : -0.000180134
    delta(WSSR) : -1.27645e-006     limit for stopping : 1e-005
    lambda : 2.17494e-007
    
    resultant parameter values
    
    a = 0.196878
    b = 107.308
    /
    
    Iteration 7
    WSSR        : 0.00708611        delta(WSSR) / WSSR : -6.06245e-007
    delta(WSSR) : -4.29592e-009     limit for stopping : 1e-005
    lambda : 2.17494e-008
    
    resultant parameter values
    
    a = 0.196921
    b = 107.368
    
    After 7 iterations the fit converged.
    final sum of squares of residuals : 0.00708611
    rel.change during last iteration : -6.06245e-007
    
    degrees of freedom(FIT_NDF) : 58
    rms of residuals(FIT_STDFIT) = sqrt(WSSR / ndf) : 0.0110532
    variance of residuals(reduced chisquare) = WSSR / ndf : 0.000122174
    
    Final set of parameters            Asymptotic Standard Error
    == == == == == == == == == == == = == == == == == == == == == == == == ==
    
    a = 0.196921 + / -0.009823     (4.988%)
    b = 107.368 + / -9.945        (9.263%)
    
    
    correlation matrix of the fit parameters :
    
                    a      b
    a               1.000
    b               0.708  1.000 
    
[/code]

調べたところ、GnuplotのfitコマンドはLevenberg-
Marquardt法を用いているようなので、EigenとGnuplotでフィッティングがうまくいく場合といかない場合があるのはなぜでしょうか？


CC BY-SA 3.0
2018-01-05T06:12:01.883
40738
2018-01-05T09:17:58.450
2018-01-05T09:17:58.450
4236
26888
question
c++ アルゴリズム
Levenberg-Marquardt法の結果がEigenとGnuplotで異なるのはなぜでしょうか