.-- --

スポンサーサイト

上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。
スポンサー広告 comment(-) trackback(-)
.20 2013

Eigenを使って非線形最小2乗法 ~Levenberg-Marquardt法 補足

こんにちわ。Uedaです。

前回の補足します。

<問題> z = a / x + b / y + c の非線形の式のa, b, cを、ノイズを乗せたM個のデータxi, yi, zi(0 < i < M)から推定する。
・正解a, b, cはそれぞれ
a = 1.0 ,b = 2.0, c = 3.0
・初期データ_a, _b, _cはそれぞれ
_a = a + noise2, _b = b + noise2, _b = b + noise
・-0.5 < xi < 0.5, -0.5 < yi < 0.5, zi = z + noise1
・-2.5 < noise1 < 2.5
・-50.0 < noise2 < 50.0

をEigenで実装されているLM法を使って推定します。

まず、データを作っていきます。これらは直接計算に関係ありません。計算の準備です。

const int DATA_COUNT = 100;
const double ANS_A = 1.0;
const double ANS_B = 2.0;
const double ANS_C = 3.0;
const double NOIZE1 = 5.0;
const double NOIZE2 = 100.0;

double *xa = new double[DATA_COUNT];
double *ya = new double[DATA_COUNT];
double *za = new double[DATA_COUNT];
for (int i = 0; i < DATA_COUNT; i++)
{
double x = 0.5 - [self random];
double y = 0.5 - [self random];
double noize = (0.5 - [self random]) * NOIZE1;
double z = ANS_A/x + ANS_B/y + ANS_C + noize;
xa[i] = x;
ya[i] = y;
za[i] = z;
}


初期値を適当に決定します。今回は正解データにnoise2を足します。(ただし, -50.0 < noise2 < 50.0

VectorXd p(3); // beta1とbeta2の初期値(適当)
double a = ANS_A + (0.5 - [self random]) * NOIZE2;
double b = ANS_B + (0.5 - [self random]) * NOIZE2;
double c = ANS_C + (0.5 - [self random]) * NOIZE2;
p(0) = a;
p(1) = b;
p(2) = c;


関数を定義します。
int operator()(const VectorXd& b, VectorXd& fvec) constの第一引数 bは求めたい係数です。今回ではa, b, cにあたります。初期計算では初期値がこれになります。fvec[i]はあるデータi(xi, yi, zi)に対して, 推定中の関数だとどれだけの誤差があるかを入れるものです。

struct misra1a_functor : Functor
{
・・・
int operator()(const VectorXd& b, VectorXd& fvec) const
{
for (int i = 0; i < values_; ++i) {
fvec[i] = z[i] - b[0]/x[i] - b[1]/y[i] - b[2]; // fvec[i] = zi - a / xi - b / yi - c
}
return 0;
}
・・・
}


実際に計算します。

misra1a_functor functor(3, DATA_COUNT, xa, ya, za);

NumericalDiff<misra1a_functor> numDiff(functor);
LevenbergMarquardt<NumericalDiff<misra1a_functor>> lm(numDiff); // 注意 "<"が全角
int info = lm.minimize(p);


1行目: 関数を生成します。
2行目: 生成した関数の(遍)微分した形の関数を用意させています。手で関数の微分をしてもいいのですが, 数値微分させて代用しています。
3, 4行目: Eigenのレーベンバーグ・マーカート法を使って計算させる。結果を配列pに出力しています。

以上です。
関連記事
スポンサーサイト

Comment

Post comment

  • comment
  • secret
  • 管理者にだけ表示を許可する

Trackback

trackbackURL:http://appteam.blog114.fc2.com/tb.php/264-ddaa0961

ブログ内検索

関連リンク

製品情報

最新記事

カテゴリ

プロフィール

neoxneo



NEXT-SYSTEM iOS Developers Blog


  • UTO:
    カナダ版iPhone4Sは、マナーモードでシャッター音がならない…


  • Ehara:
    ...


  • Hayate:
    ...


  • Tasaki:
    Developer登録完了...したのはいいけど


  • Ueda:
    ...



リンク

上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。