.-- --

スポンサーサイト

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

Eigenを利用して固有値問題を解いてみた

こんにちわ。Uedaです。

前回Eigenの導入について記載させていただきました。
せっかくなので、固有値問題を解いてみました。
(Eigenっていう名前ですし、丁度いいかなと)

単純に解いても面白くないので、
データx, yからax+by+c=0の式の係数a, b, cの推定を行うプログラムを実装し、グラフにおこして検証します。
(ただし、a^2+b^2=1
観測データから係数を推定したいことって、よくありますもんね〜

ソースコードは突貫工事ですが、掲載しておきます。



#include "Eigen/Core"
#include "Eigen/Dense"

using namespace Eigen;

-(void)viewDidAppear:(BOOL)animated // ←突貫工事の名残
{
[super viewDidAppear:animated]; // ←突貫工事の名残

const int DATA_COUNT = 200;
const double NOISE = 0.1;
double ANS_A = 0.1;
double ANS_B = 0.7;
double ANS_C = 0.7;
NSLog(@"正解 %f, %f, %f", ANS_A, ANS_B, ANS_C);
double* observeX = new double[DATA_COUNT];
double* observeY = new double[DATA_COUNT];
double mx = 0.0, my = 0.0;
for (int i = 0; i < DATA_COUNT; i++)
{
double x = 0.5 - random() / 1000000000.0;
double y = (-ANS_C - ANS_A * x) / ANS_B + NOISE * (0.5 - random() / 1000000000.0);
observeX[i] = x;
observeY[i] = y;
mx += x;
my += y;
NSLog(@"%f, %f", x, y);
}
// ↑ ここまでで観測データx, yを生成

mx /= (double) DATA_COUNT;
my /= (double) DATA_COUNT;
double b11 = 0.0;
double b12 = 0.0;
double b22 = 0.0;
for (int i = 0; i < DATA_COUNT; i++)
{
double dx = observeX[i] - mx;
double dy = observeY[i] - my;
b11 += dx * dx;
b12 += dx * dy;
b22 += dy * dy;
}
Matrix2d B;
B(0, 0) = b11; B(0, 1) = b12;
B(1, 0) = b12; B(1, 1) = b22;
SelfAdjointEigenSolver eig(B); // ←固有値計算
double a = -eig.eigenvectors()(0, 0); // ←今回は固有ベクトルが推定値になります
double b = -eig.eigenvectors()(1, 0);
double c = -a*mx - b*my;
double n = sqrt(a*a + b*b + c*c);
NSLog(@"推定 %f, %f, %f", a/n, b/n, c/n);
}


検証グラフ
固有値問題を解いてみた


観測データのy成分が変に小さい気が・・・(観測データ生成してるとこにバグがありそう
でも、固有値問題を解くことで観測データから正解の式を推定できていることがわかりますね!

行列計算ライブラリ(Eigen)を使って、見事数学の恩恵を簡単に承ることができました!

固有値問題の計算について、これ参考にしました
関連記事
スポンサーサイト

Comment

Post comment

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

Trackback

trackbackURL:http://appteam.blog114.fc2.com/tb.php/246-5aabb9be

ブログ内検索

関連リンク

製品情報

最新記事

カテゴリ

プロフィール

neoxneo



NEXT-SYSTEM iOS Developers Blog


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


  • Ehara:
    ...


  • Hayate:
    ...


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


  • Ueda:
    ...



リンク

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