.-- --

スポンサーサイト

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

K-means法理解

こんにちわ。Uedaです。

今回はクラスタリングの代表的なアルゴリズムであるk-means法を検証してみました。

参考サイト
http://tech.nitoyon.com/ja/blog/2009/04/09/kmeans-visualise/

ほうほう。では早速実装。


...
#ifdef _k_means_

std::vector dataArray;
const int DATA_COUNT = 100;
const int CLASS_COUNT = 5;
for (int i = 0 ; i < DATA_COUNT; i++)
{
double *data = new double[3];
data[0] = [self random] * 10.0;
data[1] = [self random] * 10.0;
double classVal = [self random];
for (int j = 0; j < CLASS_COUNT; j++)
{
if ( classVal < ((j+1) / (double)CLASS_COUNT) )
{
data[2] = j;
break;
}
}
dataArray.push_back(data);
NSLog(@"%f,%f,%d, %f", data[0], data[1], (int)data[2], classVal);
}

std::vectorcenterArray;
for (int i = 0; i < CLASS_COUNT; i++)
{
double *center = new double[2];
centerArray.push_back(center);
}

Boolean bFinish = false;
int calcCount = 0;
while (!bFinish)
{
for (int i = 0; i < CLASS_COUNT; i++)
{
double *center = centerArray.at(i);
double x = 0.0, y = 0.0;
int count = 0;
for (int j = 0; j < DATA_COUNT; j++)
{
double *data = dataArray.at(j);
if ((int)data[2] == i)
{
x += data[0];
y += data[1];
count++;
}
}
center[0] = x / (double)count;
center[1] = y / (double)count;
}

bFinish = true;
for (int i = 0; i < DATA_COUNT; i++)
{
double *data = dataArray.at(i);
double minDist = 0.0;
int min = -1;
int prev = (int) data[2];
for (int j = 0; j < CLASS_COUNT; j++)
{
double* center = centerArray.at(j);
double dx = data[0] - center[0];
double dy = data[1] - center[1];
double dist = sqrt(dx*dx + dy*dy);
if (min < 0 || dist < minDist)
{
min = j;
minDist = dist;
data[2] = j;
}
}
int cur = (int) data[2];
if (prev != cur) bFinish = false;
}
calcCount++;
}
NSLog(@"計算回数 %d", calcCount);

for (int i = 0 ; i < DATA_COUNT; i++)
{
double *data = dataArray.at(i);
NSLog(@"%f,%f,%d", data[0], data[1], (int)data[2]);
}
NSLog(@"重心");
for (int i = 0 ; i < CLASS_COUNT; i++)
{
double *center = centerArray.at(i);
NSLog(@"%f,%f,%d", center[0], center[1], i);
}

#endif

...

検証データ

Before

k-means_prev.png

After

k-means_cur.png

みごと、データを5つのクラスに分けることができました。

スポンサーサイト
.26 2013

Eigenでカルマンフィルタ理解

こんにちわ。Uedaです。

今日はEigenを使ってカルマンフィルタの検証を行います。

wiki
http://ja.wikipedia.org/wiki/カルマンフィルター

カルマンフィルタを軽く説明します。
wiki曰く、
「誤差のある観測値を用いて、ある動的システムの状態を推定あるいは制御するための、無限インパルス応答フィルターの一種である。」
ということらしいです。

例題から入った方がわかりやすいと思いますので、この度課題とさせて頂いたサイトを紹介
http://satomacoto.blogspot.jp/2011/06/python.html

以下で筆者が実装したプログラムをお見せしますが、このサイトのプログラムをコピペした訳ではないですよ!
なにせい、筆者はPython読めませんw

例題を読みましょう。
2次元座標において、あるロボットがt=0に原点を出発して、速度(2,2)で動くとする。ロボットの進路は風などの影響を受け(σx=σy=1)、毎秒ごとに観測できるGPSによる位置座標には計測誤差(σx=σy=2)があるとする。このとき、観測された軌跡から実際の軌跡を推定する。

なるほど、問題はロボットの軌跡を推定することです。与えられるものは、ロボットの制御値(例題でいうところの速度(2,2)で動く)と観測値(例題でいうところのGPSによる位置座標)です。これをもとにロボットの軌跡を推定していくわけです。

しかし!!! ロボットの制御にも、観測値にも誤差が含まれます。(どちらか一方が誤差なしなら、そもそも推定する必要がありませんよね

ようは、
制御値だけでは風の影響(誤差)があるので、精度のいい位置が計測できない。
観測値だけではGPSには誤差が含まるので、精度のいい位置が計測できない。
っということで、二つ合わせてもっと精度のいい推定をしよう!このセンサの統合にカルマンフィルタを使おう!ということです。

では、プログラムと検証データです。


#ifdef _kalmanfilter_

const int DATA_COUNT = 10;
// const double NOIZE1 = 2.0;
// const double NOIZE2 = 4.0;
const double NOIZE1 = 4.0;
const double NOIZE2 = 8.0;
std::vector correctList;
std::vector observeList;
std::vector controlList;


const double VEL = 2.0;
double correctX = 0.0, correctY = 0.0;
for (int t = 0; t < DATA_COUNT; t++)
{

correctX += VEL + [self random] * NOIZE1;
correctY += VEL + [self random]* NOIZE1;
double* correct = new double[2];
correct[0] = correctX;
correct[1] = correctY;
correctList.push_back(correct);

double* control = new double[2];
control[0] = VEL;
control[1] = VEL;
controlList.push_back(control);

double* observe = new double[2];
observe[0] = correctX + [self random] * NOIZE2;
observe[1] = correctY + [self random] * NOIZE2;
observeList.push_back(observe);
}

Vector2d mu;
mu(0) = 0.0; mu(1) = 0.0;
Matrix2d sigma = Matrix2d::Zero(2, 2);
Matrix2d Q = Matrix2d::Identity(2, 2);
Matrix2d R = Matrix2d::Zero(2, 2);
R(0, 0) = 2.0; R(1, 1) = 2.0;
for (int i = 0; i < DATA_COUNT; i++)
{
double *correct = correctList.at(i);
double *observe = observeList.at(i);
double *control = controlList.at(i);

// 予測
mu(0) += control[0];
mu(1) += control[1];
sigma += Q;

// 更新
MatrixXd y(2, 1);
y(0, 0) = observe[0] - mu(0);
y(1, 0) = observe[1] - mu(1);
MatrixXd S = sigma + R;
MatrixXd K = sigma * S.inverse();
MatrixXd _mu(2, 1);
_mu(0, 0) = mu(0);
_mu(1, 0) = mu(1);
MatrixXd ans = _mu + K * y;
mu(0) = ans(0, 0);
mu(1) = ans(1, 0);
sigma -= K * sigma;

NSLog(@"%f,%f,%f,%f,%f,%f", correct[0], correct[1], observe[0], observe[1], mu(0), mu(1));
}

#endif


検証データ
カルマンフィルタ

観測値だけの赤いグラフより、カルマンフィルタで推定した緑のグラフの方が、青い実際にロボットが動いた軌跡に近いことがわかります!
カルマンフィルタは、wikiで紹介されている式に当てはめることができれば、異なるデータをブレンドしてよりいい (コーヒーが) 結果を得る (飲める) ことができる。というフィルタです。(ただし、こういったフィルタは他にもいろいろあります〜
.25 2013

Eigenでクォータニオン理解

こんにちわ。Uedaです。

今回はEigenを使ってクォータニオンを作ってみました。
前回、回転行列を作ったので、いけるだろーと思っていました...
が、
符号が合わないやら、行列の要素逆に格納してしまうわで、てこずってしまいました。。

以下コードです。


#ifdef _quaternion_

double p = 0.5 -[self random] / 10.0;
double y = 0.5 - [self random] / 10.0;
double r = 0.5 - [self random] / 10.0;
Matrix3d R_pitch;
R_pitch(0, 0) = 1.0; R_pitch(0, 1) = 0.0; R_pitch(0, 2) = 0.0;
R_pitch(1, 0) = 0.0; R_pitch(1, 1) = cos(p); R_pitch(1, 2) = -sin(p);
R_pitch(2, 0) = 0.0; R_pitch(2, 1) = sin(p); R_pitch(2, 2) = cos(p);
Matrix3d R_yaw;
R_yaw(0, 0) = cos(y); R_yaw(0, 1) = 0.0; R_yaw(0, 2) = sin(y);
R_yaw(1, 0) = 0.0; R_yaw(1, 1) = 1.0; R_yaw(1, 2) = 0.0;
R_yaw(2, 0) = -sin(y); R_yaw(2, 1) = 0.0; R_yaw(2, 2) = cos(y);
Matrix3d R_roll;
R_roll(0, 0) = cos(r); R_roll(0, 1) = -sin(r); R_roll(0, 2) = 0.0;
R_roll(1, 0) = sin(r); R_roll(1, 1) = cos(r); R_roll(1, 2) = 0.0;
R_roll(2, 0) = 0.0; R_roll(2, 1) = 0.0; R_roll(2, 2) = 1.0;
MatrixXd R = R_roll * R_pitch * R_yaw;

double q0, q1, q2, q3;
double r11 = R(0, 0); double r12 = R(0, 1); double r13 = R(0, 2);
double r21 = R(1, 0); double r22 = R(1, 1); double r23 = R(1, 2);
double r31 = R(2, 0); double r32 = R(2, 1); double r33 = R(2, 2);

q0 = ( r11 + r22 + r33 + 1.0f) / 4.0f;
q1 = ( r11 - r22 - r33 + 1.0f) / 4.0f;
q2 = (-r11 + r22 - r33 + 1.0f) / 4.0f;
q3 = (-r11 - r22 + r33 + 1.0f) / 4.0f;
if(q0 < 0.0f) q0 = 0.0f;
if(q1 < 0.0f) q1 = 0.0f;
if(q2 < 0.0f) q2 = 0.0f;
if(q3 < 0.0f) q3 = 0.0f;
q0 = sqrt(q0);
q1 = sqrt(q1);
q2 = sqrt(q2);
q3 = sqrt(q3);
if(q0 >= q1 && q0 >= q2 && q0 >= q3) {
q0 *= +1.0f;
q1 *= 1.0 - 2.0 * signbit(r32 - r23);
q2 *= 1.0 - 2.0 * signbit(r13 - r31);
q3 *= 1.0 - 2.0 * signbit(r21 - r12);
} else if(q1 >= q0 && q1 >= q2 && q1 >= q3) {
q0 *= 1.0 - 2.0 * signbit(r32 - r23);
q1 *= +1.0f;
q2 *= 1.0 - 2.0 * signbit(r21 + r12);
q3 *= 1.0 - 2.0 * signbit(r13 + r31);
} else if(q2 >= q0 && q2 >= q1 && q2 >= q3) {
q0 *= 1.0 - 2.0 * signbit(r13 - r31);
q1 *= 1.0 - 2.0 * signbit(r21 + r12);
q2 *= +1.0f;
q3 *= 1.0 - 2.0 * signbit(r32 + r23);
} else if(q3 >= q0 && q3 >= q1 && q3 >= q2) {
q0 *= 1.0 - 2.0 * signbit(r21 - r12);
q1 *= 1.0 - 2.0 * signbit(r31 + r13);
q2 *= 1.0 - 2.0 * signbit(r32 + r23);
q3 *= +1.0f;
} else {
}
double normalize = sqrt(q0*q0 + q1*q1 + q2*q2 + q3*q3);
q0 /= normalize;
q1 /= normalize;
q2 /= normalize;
q3 /= normalize;

double sx = q0 * q0;
double sy = q1 * q1;
double sz = q2 * q2;
double cx = q1 * q2;
double cy = q0 * q2;
double cz = q0 * q1;
double wx = q3 * q0;
double wy = q3 * q1;
double wz = q3 * q2;

Matrix3d Rq(3, 3);

Rq(2, 2) = 1.0 - 2.0 * (sy + sz);
Rq(2, 1) = 2.0 * (cz + wz);
Rq(2, 0) = -(2.0 * (cy - wy));
Rq(1, 2) = -(2.0 * (cz - wz));
Rq(1, 1) = -(1.0 - 2.0 * (sx + sz));
Rq(1, 0) = 2.0 * (cx + wx);
Rq(0, 2) = 2.0 * (cy + wy);
Rq(0, 1) = 2.0 * (cx - wx);
Rq(0, 0) = -(1.0 - 2.0 * (sx + sy));

for (int i = 0; i < 3; i++)
{
NSLog(@"%f,%f,%f", R(i, 0), R(i, 1), R(i, 2));
}
NSLog(@"\n");
for (int i = 0; i < 3; i++)
{
NSLog(@"%f,%f,%f", Rq(i, 0), Rq(i, 1), Rq(i, 2));
}

#endif


出力結果
2013-04-25 17:25:01.295 TestEigen[5495:c07] 0.823784,-0.309281,0.475105
2013-04-25 17:25:01.297 TestEigen[5495:c07] 0.417921,0.897580,-0.140332
2013-04-25 17:25:01.297 TestEigen[5495:c07] -0.383043,0.314159,0.868667
2013-04-25 17:25:01.298 TestEigen[5495:c07]
2013-04-25 17:25:01.299 TestEigen[5495:c07] 0.823784,-0.309281,0.475105
2013-04-25 17:25:01.300 TestEigen[5495:c07] 0.417921,0.897580,-0.140332
2013-04-25 17:25:01.300 TestEigen[5495:c07] -0.383043,0.314159,0.868667

値は合ってますが、すみません、プログラムは自信ありません。。
.24 2013

Eigenで回転行列理解

やってしまいました...
本日寝坊してしまった、Uedaです。
目覚まし時計、もう一個買っておこう..

さて、今回はまたEigenにもどりまして3次元座標系での回転行列を作ってみました。

クォータニオンは次にとっておくとして、今回は各軸回転(yaw, pitch, roll)についてのプログラムを備忘録として残しておくことにします。

pitch → x軸中心の回転
yaw → y軸中心の回転
roll → z軸中心の回転

以下のコードは、
1. ランダムにpitch, yaw, rollの角度を算出
2. 算出した角度で回転行列を作る
3. 作った回転行列から、pitch, yaw, rollを抽出。元データとの確認。
4. 回転行列の行列としての性質として、ある座標を掛けて、掛けた座標を回転行列の転置行列を掛けると元通りになるというものがあります。その検証。

の順で、処理しています。


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

...
#ifdef _rotate_matrix_

NSLog(@"回転行列のテスト");

double p = 0.5 - [self random] / 10.0;
double y = 0.5 - [self random] / 10.0;
double r = 0.5 - [self random] / 10.0;

Matrix3d R_pitch;
R_pitch(0, 0) = 1.0; R_pitch(0, 1) = 0.0; R_pitch(0, 2) = 0.0;
R_pitch(1, 0) = 0.0; R_pitch(1, 1) = cos(p); R_pitch(1, 2) = -sin(p);
R_pitch(2, 0) = 0.0; R_pitch(2, 1) = sin(p); R_pitch(2, 2) = cos(p);
Matrix3d R_yaw;
R_yaw(0, 0) = cos(y); R_yaw(0, 1) = 0.0; R_yaw(0, 2) = sin(y);
R_yaw(1, 0) = 0.0; R_yaw(1, 1) = 1.0; R_yaw(1, 2) = 0.0;
R_yaw(2, 0) = -sin(y); R_yaw(2, 1) = 0.0; R_yaw(2, 2) = cos(y);
Matrix3d R_roll;
R_roll(0, 0) = cos(r); R_roll(0, 1) = -sin(r); R_roll(0, 2) = 0.0;
R_roll(1, 0) = sin(r); R_roll(1, 1) = cos(r); R_roll(1, 2) = 0.0;
R_roll(2, 0) = 0.0; R_roll(2, 1) = 0.0; R_roll(2, 2) = 1.0;

MatrixXd R = R_roll * R_pitch * R_yaw;

double r11 = R(0, 0); double r12 = R(0, 1); double r13 = R(0, 2);
double r21 = R(1, 0); double r22 = R(1, 1); double r23 = R(1, 2);
double r31 = R(2, 0); double r32 = R(2, 1); double r33 = R(2, 2);

double yaw, roll, pitch;
double cp = sqrt(r12*r12 + r22*r22);

if (cp == 0.0)
{
pitch = r32/M_PI_2;
pitch = r32*M_PI_2;
yaw = 0.0;
roll = atan2(r21, r11);
}
else
{
pitch = atan2(r32, cp);
yaw = atan2(-r31, r33);
roll = atan2(-r12, r22);
}

NSLog(@"正解 %f, %f, %f", p, r, y);
NSLog(@"計算 %f, %f, %f", pitch, roll, yaw);

MatrixXd data(3, 1);
data(0, 0) = 0.5 - [self random] / 10.0;
data(1, 0) = 0.5 - [self random] / 10.0;
data(2, 0) = 0.5 - [self random] / 10.0;
MatrixXd data2 = R * data;
MatrixXd data3 = R.transpose() * data2;

NSLog(@"写像前 %f, %f, %f", data(0, 0), data(1, 0), data(2, 0));
NSLog(@"写像後 %f, %f, %f", data3(0, 0), data3(1, 0), data3(2, 0));

#endif
...


出力
回転行列のテスト
正解 0.319571, 0.331831, 0.415307
計算 0.319571, 0.331831, 0.415307
写像前 0.328536, 0.304225, 0.457576
写像後 0.328536, 0.304225, 0.457576

すばらしい!ぴったりです♪

参考記事
http://www.japan-iss.co.jp/wp-content/uploads/2010/02/RollPitchYaw.pdf
.23 2013

objローダー探検記 ~iPhone

こんにちわ。Uedaです。

ARでどうしても通らないといけないのが3D描画。

iPhone上で3D描画したいなーと思いましたが、
慣れないObjective-CでOpenGLがりがり書くのは面倒で、
とりあえず、「簡単に」3D描画できないか模索・探索してみました。。

<3Dモデル描画は表現手段でアルゴリズム検証に使うので、とりあえず描画できればいい!って感じです。あと、無料なやつ。

まず、見つけたのがこれです。Wave Front OBJ
http://bill.dudney.net/roller/objc/entry/wave_front_obj_textures_working
↑ hereってのがリンク

obj読めたら、まー便利だろうと思い、
早速突っ込んでみましたが、アルゴリズムで推定した位置・姿勢をOpenGL側に適切に渡さなければいけないので、バグ混入がすごく不安...
アルゴリズムも研究最中のたいしたことないので、アルゴリズムのバグかもしれないですしw
まー、勉強すればいいじゃんってな感じですが、とりあえず保留しちゃいました。
(でも、結構便利だと思いましたよ!

次に見つけたのがassimp
http://assimp.sourceforge.net
調べてみるとassimpは読み込めるフォーマットが多いです!!
すごく賢い子なのですが、導入に手こずりそうだったので、とりあえず保留しちゃいました。

最後に見つけたのがNinevehGLという3D描画エンジン。
http://nineveh.gl
結論からいうと、こいつは半端ではない使いやすさです!!
3D描画が4、5行で実装できます!
beta版ですが、Tutorial動画があり取っ付きやすいです。こいつのせいで昨日は少し夜更かししてしまいましたw

サンプルコードいじってみて紹介したかったのですが、会社のコンピュータでインストールする訳にはいかないので、今日ひとりでいじって楽しみますw
お気に入りobjモデルとiPhoneがあれば一度試してみてはいかがでしょうか。
.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

この関数で表せるものっていうとー...
「いい塩梅な塩加減」っとかでしょうか (お後がよろしいようで
.19 2013

Eigenを使って最小2乗法を解いてみた

こんにちわ。Uedaです。

前回はEigenを使って固有値問題を解きました。今回は最小2乗法を解いてみます。

最小2乗法はフィッティング問題の常套方法だと思います。
問題へのアプローチが、誤差の2乗和を最小にするというもので直感的なので筆者は大好きです。

今回はストーリー仕立てでお送りします。

Aさんがどんなお味噌汁が好きか調査します。
お味噌汁の味付けは塩・味噌・鰹だしとし、Aさんはどの割合で調理したお味噌汁が好みか調べます。
それぞれの調味料の分量を変えたお味噌汁を飲んでもらっておいしかったかどうか聞きます。
例えば、
塩 1g, 味噌 0g, 鰹だし 0g. Aさんの評価→最悪 (そらそうだw
塩 1g, 味噌 1g, 鰹だし 0g. Aさんの評価→さっきよりはまし!
みたいなことを繰り返し行います。

塩がx, 味噌がy, 鰹だしがz, Aさんの評価をwとします。
ちょっとはしょりますが、
最小2乗法でAさんの好みを推定するプログラムは以下になります!
当然、Aさんの評価や調味料の分量には誤差が含まれるので、これらにはノイズを乗せてます。
それでは調査のためにAさんにはお味噌汁を100杯飲んでもらいましょう〜w



・・・・・

const double ANS_A = 1.0;
const double ANS_B = 2.0;
const double ANS_C = 3.0;
const double ANS_D = 4.0;
const int DATA_COUNT = 100;
const double NOISE = 1.0;

double* xArr = new double[DATA_COUNT];
double* yArr = new double[DATA_COUNT];
double* zArr = new double[DATA_COUNT];
double* wArr = new double[DATA_COUNT];

// Aさんに100杯お味噌汁を飲んでもらって、実験開始!! ------------------
for (int i = 0; i < DATA_COUNT; i++) {
double x = 0.5 - random() / 1000000000.0;
double y = 0.5 - random() / 1000000000.0;
double z = 0.5 - random() / 1000000000.0;
double noize = (0.5 - random() / 1000000000.0) * NOISE;
double w = x * ANS_A + y * ANS_B + z * ANS_C + ANS_D + noize;
wArr[i] = w; // Aさんの評価
xArr[i] = x; // 塩の分量
yArr[i] = y; // 味噌の分量
zArr[i] = z; // 鰹だしの分量
NSLog(@"%f, %f, %f, %f, %f", x, y, z, w, noize);
}
// Aさんご苦労様 ------------------

// Eigen使って最小2乗法 ------------
MatrixXd X(DATA_COUNT, 4);
MatrixXd Y(DATA_COUNT, 1);
for (int i = 0; i < DATA_COUNT; i++) {
X(i, 0) = xArr[i];
X(i, 1) = yArr[i];
X(i, 2) = zArr[i];
X(i, 3) = 1.0;
Y(i, 0) = wArr[i];
}
MatrixXd Xt = X.transpose();
MatrixXd ans = (Xt * X).inverse() * Xt * Y;

// 出力 ------------
NSLog(@"正解 %f, %f, %f, %f", ANS_A, ANS_B, ANS_C, ANS_D);
// → 正解 1.000000, 2.000000, 3.000000, 4.000000
NSLog(@"推定 %f, %f, %f, %f", ans(0, 0), ans(1, 0), ans(2, 0), ans(3, 0));
// → 推定 0.891576, 2.057038, 3.063395, 3.384292

・・・・・


結果は、
正解 1.000000, 2.000000, 3.000000, 4.000000
推定 0.891576, 2.057038, 3.063395, 3.384292
ですかー。
Aさんには少し塩っけの足りないみそ汁になっちゃいそうですが、おおむね推定できてると思います。

今回はEigenを使って最小2乗法を解くということで、
最小2乗法の計算方法の中身の紹介ができませんでした。
以下のサイトが詳しいです。
http://itmath.blog60.fc2.com/blog-entry-111.html

MatrixXd Xt = X.transpose();
MatrixXd ans = (Xt * X).inverse() * Xt * Y;

このサイト上の a = (Gt G)^-1 Gt Y
が味噌です! (お後がよろしいようで
.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)を使って、見事数学の恩恵を簡単に承ることができました!

固有値問題の計算について、これ参考にしました
.17 2013

行列計算ライブラリ「Eigen」の導入

はじめまして、Uedaです。
ただいまObjective-C+iPhone開発の勉強中です。

みなさん、行列計算は好きでしょうか?私はあまり好きではありません。
そんな私が、Objective-Cかじりたての私が、Objective-Cで「逆行列計算」「固有値計算」「特異値分解」...なんかを自前で実装するのは非常に危険です。
ベースにしたい計算にバグ混入の恐れや、速度を気にしたくありませんし...

そもそもXcodeにはAccelerateなるもで行列計算できるみたいですが、
乗算したいだけなのに関数呼ばないとならないらしく、大変そうです...
ちなみに固有値計算するだけでもこんな感じ...
http://sakura.math.kyushu-u.ac.jp/wiki/index.php?計算機設定%2FCLAPACK%2FXcode

そこで、C++で実装されたEigenを導入しました。
(ちなみにJava使っていたときはJamaを使っていました。)

Eigenはオペレータ演算子(+,-,*)で行列の足し算、引き算、積ができ、コードもすっきりします。
また、Eigenを使うと固有値計算はこんな感じになります。


int main()
{
Matrix2f A;
A << 1, 2, 2, 3; // 行列の初期化 (Matlabっぽく書くと [1 2; 2 3])
SelfAdjointEigenSolver eigensolver(A); // ←固有値計算
eigensolver.eigenvalues(); // 固有値
eigensolver.eigenvectors(); // 固有ベクトル
}


行列計算がとてもすっきりします!

Eigen
Jama
iPhoneへの導入方法
 HOME 

ブログ内検索

関連リンク

製品情報

最新記事

カテゴリ

プロフィール

neoxneo



NEXT-SYSTEM iOS Developers Blog


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


  • Ehara:
    ...


  • Hayate:
    ...


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


  • Ueda:
    ...



リンク

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