Levenberg-Marquardt法の結果がEigenとGnuplotで異なるのはなぜでしょうか
C++の数値計算ライブラリEigenに含まれる、非線形最小二乗法(Levenberg-Marquardt法)で関数フィッティングを実行してみました
開発環境はVisual Studio 2015です
フィッティング対象の関数は、y = a×cos(x)^b で、ノイズを含む観測データx,yからパラメータa,bを求めます
#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();
}
実行結果は以下のようになりました
TooManyFunctionEvaluation
0.6, 30
収束しておらずパラメータも初期値のままなので、うまくいっていない様子です
しかし、Gnuplot(version 4.6)のfitコマンドを用いて同じデータ(保存した"output1.txt")でフィッティングを行ったところ、収束をし、パラメータも真値に近づきました
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
調べたところ、GnuplotのfitコマンドはLevenberg-Marquardt法を用いているようなので、EigenとGnuplotでフィッティングがうまくいく場合といかない場合があるのはなぜでしょうか？