.-- --

スポンサーサイト

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

Eigenで最小2乗法(2) 〜ロバスト推定

こんにちわ。Uedaです。

最小2乗法について、最初の方で触れました。前回、たくさんの観測データ(ノイズを含む)から、その観測データの傾向を抽出するという最適化問題として、利用しました。

今回は観測データに大きめのノイズがのったデータを含ませます。このノイズだらけのデータから、ノイズの影響をあまり受けないように推定する方法を検証します。

ちなみに、ノイズの影響をあまり受けずに推定する、ノイズに強い推定のことを、ロバストな推定、ロバスト推定といいます。

扱うのはTukeyのBiweight推定法です。
TukeyのBiweight推定法について
http://imagingsolution.blog107.fc2.com/blog-entry-32.html

では、検証実験内容です。
w=ax+by+cz+dの式の係数a, b, c, dをノイズの乗ったデータw, x, y, zから推定する。(a, b, c, dは1.0, 2.0, 3.0, 4.0とした) データw, x, y, zはそれぞれ130個あり、
100個は-0.5〜0.5のノイズをのせる。
30個は-50〜50のノイズをのせる。

TukeyのBiweight推定法は重み付け最小2乗法を繰り返し行って、解に近づけていきます。推定した係数をもとに、あるデータがその係数からなる線からどれだけ遠いか(どれだけの誤差があるか)を調べ、近いデータ(誤差が小さい)ほど重みを軽くし、次の最小2乗法での計算で重視します。(筆者の理解が正しければ...
今回繰り返し回数は3回にし、重みを決定するための値を40.0としました。

以下コードと出力です。


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_A = 100;
const int DATA_COUNT_B = 30;
const int DATA_COUNT = DATA_COUNT_A + DATA_COUNT_B;
const double NOISE_A = 1.0;
const double NOISE_B = 100.0;
const int CALC_ITERATE_COUNT = 3;
const double W = 40.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];

double* weightArr = new double[DATA_COUNT];

for (int i = 0; i < DATA_COUNT_A; i++) {
double x = 0.5 - [self random];
double y = 0.5 - [self random];
double z = 0.5 - [self random];
double noize = (0.5 - [self random]) * NOISE_A;
double w = x * ANS_A + y * ANS_B + z * ANS_C + ANS_D + noize;
wArr[i] = w;
xArr[i] = x;
yArr[i] = y;
zArr[i] = z;
weightArr[i] = 1.0;
NSLog(@"%f, %f, %f, %f, %f", x, y, z, w, noize);
}

for (int i = DATA_COUNT_A; i < DATA_COUNT; i++) {
double x = 0.5 - [self random];
double y = 0.5 - [self random];
double z = 0.5 - [self random];
double noize = (0.5 - [self random]) * NOISE_B;
double w = x * ANS_A + y * ANS_B + z * ANS_C + ANS_D + noize;
wArr[i] = w;
xArr[i] = x;
yArr[i] = y;
zArr[i] = z;
weightArr[i] = 1.0;
NSLog(@"%f, %f, %f, %f, %f", x, y, z, w, noize);
}


MatrixXd X(DATA_COUNT, 4);
MatrixXd Y(DATA_COUNT, 1);
MatrixXd ans;

for (int j = 0; j < CALC_ITERATE_COUNT; j++)
{
for (int i = 0; i < DATA_COUNT; i++)
{
X(i, 0) = xArr[i] * weightArr[i];
X(i, 1) = yArr[i] * weightArr[i];
X(i, 2) = zArr[i] * weightArr[i];
X(i, 3) = 1.0 * weightArr[i];
Y(i, 0) = wArr[i] * weightArr[i];
}

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

for (int i = 0; i < DATA_COUNT; i++)
{
double ax = xArr[i] * ans(0, 0);
double by = yArr[i] * ans(1, 0);
double cz = zArr[i] * ans(2, 0);
double d = ans(3, 0);
double w = wArr[i];
double e = w - ax - by - cz - d;
// NSLog(@"%f", e);
if (fabs(e) > W) weightArr[i] = 0.0;
else
{
e /= W;
weightArr[i] = (1.0 - e*e) * (1.0 - e*e);
}
}
NSLog(@"推定 %f, %f, %f, %f", ans(0, 0), ans(1, 0), ans(2, 0), ans(3, 0));
}

NSLog(@"正解 %f, %f, %f, %f", ANS_A, ANS_B, ANS_C, ANS_D);


出力
Msuitei2.png


一番目の推定結果はふつうに最小2乗法を解いた場合です。繰り返し計算することにより解に近づいていることがわかります。


関連記事
スポンサーサイト

Comment

Post comment

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

Trackback

trackbackURL:http://appteam.blog114.fc2.com/tb.php/261-7f105486
-
管理人の承認後に表示されます
2013.10.17 17:59 

ブログ内検索

関連リンク

製品情報

最新記事

カテゴリ

プロフィール

neoxneo



NEXT-SYSTEM iOS Developers Blog


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


  • Ehara:
    ...


  • Hayate:
    ...


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


  • Ueda:
    ...



リンク

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