.-- --

スポンサーサイト

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

Eigenで最小2乗法を解いてみた 〜非線形編

こんにちわ。Uedaです。
今日は電話対応と名刺交換の研修をしてもらいましたー
まだ、自分の名刺を見るとにやけてしまいますし、
電話対応もままならないです...


さて、前回最小2乗法について書きました。
あれは線形でしたが、世の中はあんなに奇麗じゃないですよねー。

お味噌汁の味付けの数学モデルでさえ、あんなに単純に作れないかもしれません。
<味噌の濃いお味噌汁には、鰹だしはマイナス効果...みたいなことです。

ということで、今回は非線形の最小2乗法をやってみたいと思います。
もちろんEigenを使って!!
用いた式は y = a x^2 + b x + cです。
<この数式で表される場合も奇麗な場合だと思いますが...

ちなみに筆者はこの関数を扱うのは初めてですが、
行列ライブラリEigenを使うとものの30分で検証完了しました♪
先人のお知恵は骨も残さず活用すべきです!!!


...

const int DATA_COUNT = 200;
const double NOISE = 0.5;
double ANS_A = 2.0;
double ANS_B = 3.0;
double ANS_C = 4.0;
NSLog(@"正解 %f, %f, %f", ANS_A, ANS_B, ANS_C);
double* observeX = new double[DATA_COUNT];
double* observeY = new double[DATA_COUNT];
double X1 = 0.0, X2 = 0.0, X3 = 0.0, X4 = 0.0;
for (int i = 0; i < DATA_COUNT; i++)
{
// データ生成 ↓
double x = 0.5 - random() / 1000000000.0;
double y = ANS_A * x * x + ANS_B * x + ANS_C+ NOISE * (0.5 - random() / 1000000000.0);
observeX[i] = x;
observeY[i] = y;
NSLog(@"%f, %f", x, y);
   // データ生成 ↑
X1 += x;
X2 += x*x;
X3 += x*x*x;
X4 += x*x*x*x;
}
MatrixXd A(3, 3);
double n = (double)DATA_COUNT;
A(0, 0) = n; A(0, 1) = X1; A(0, 2) = X2;
A(1, 0) = X1; A(1, 1) = X2; A(1, 2) = X3;
A(2, 0) = X2; A(2, 1) = X3; A(2, 2) = X4;
MatrixXd B(3, DATA_COUNT);
MatrixXd C(DATA_COUNT, 1);
for (int i = 0; i < DATA_COUNT; i++)
{
double x = observeX[i];
double y = observeY[i];
B(0, i) = 1.0;
B(1, i) = x;
B(2, i) = x*x;
C(i, 0) = y;
}
MatrixXd ans = A.inverse() * B * C; // ← 最小2乗法
NSLog(@"%f,%f,%f", ans(0, 0), ans(1, 0), ans(2, 0));

...


検証グラフ
最小2乗法(非線形)

正解 2.000000, 3.000000, 4.000000
推定 1.979831, 2.960393, 3.695208

この関数で表せるものっていうとー...
「いい塩梅な塩加減」っとかでしょうか (お後がよろしいようで
関連記事
スポンサーサイト

Comment

Post comment

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

Trackback

trackbackURL:http://appteam.blog114.fc2.com/tb.php/248-8df72a44

ブログ内検索

関連リンク

製品情報

最新記事

カテゴリ

プロフィール

neoxneo



NEXT-SYSTEM iOS Developers Blog


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


  • Ehara:
    ...


  • Hayate:
    ...


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


  • Ueda:
    ...



リンク

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