http://mathforum.org/library/drmath/view/69239.html


ここで説明されている解を得るまでのステップについて、実数を使って実際の流れをみせてもらえないでしょうか。

(英語はわかりますが、計算のほうのオペレーションがわかりません。)

回答の条件
  • 1人1回まで
  • 登録:2009/09/12 10:39:39
  • 終了:2009/09/15 16:07:35

ベストアンサー

id:ita No.2

ita回答回数204ベストアンサー獲得回数482009/09/13 19:59:14

ポイント46pt

ーーー計算のバックグラウンドーーー

i番目のデータの実測値と推定値のズレ Eiは

Ei = a*(xi-x0) + b*(yi-y0) + z0 -zi

これの二乗和 E=ΣEi^2 を最小にするaとbをもとめます。それにはEをaで微分して0かつEをbで微分して0

という条件を満たすa、bを求めます。

aで微分すると (d/da) ΣEi^2 = 2ΣEi (d/da)Ei = 2Σ (a*(xi-x0)+b*(yi-y0) + z0 - zi)*(xi-x0)

=2aΣ(xi-x0)(xi-x0) +2bΣ(yi-y0)(xi-x0)+2Σ(z0-zi)(xi-x0)

Σ(xi-x0)(yi-y0)=Vxyなどとおくと

=2a Vxx + 2b Vxy -2Vzx =0 これが第一の式

bで微分すると同様な計算で

2a Vxy + 2b Vyy - 2Vyz = 0

行列で書くと

 |Vxx  Vxy | |a|   |Vzx|
 |Vxy  Vyy | |b| = |Vzy|
行列の逆行列を|Vzx Vzy|にかけるとa,bが出る。逆行列は
   1              |  Vyy  -Vxy |
 ---------------- |            |
 Vxx Vyy - Vxy^2  |  -Vxy  Vxx |

したがって最終的に

a= (Vyy*Vzx - Vxy*Vzy)/(Vxx*Vyy-Vxy*Vxy)

b= (-Vxy*Vzx +Vxx*Vzy)/(Vxx*Vyy-Vxy*Vxy)

Vxxなどは以下のように計算

double x0=0, y0=0, z0=0
for(i=0;i<N;i++){
x0 += x[i];
y0 += y[i];
z0 += z[i];
}
x0/=N;
y0/=N;
z0/=N;

double Vxx=0, Vyy=0, Vzz=0, Vxy=0, Vyz=0, Vzx=0;
for(i=0;i<N;i++)
{
  double dx = x[i]-x0;
  double dy = y[i]-y0;
  double dz = z[i]-z0;

  Vxx += dx*dx;
  Vyy += dy*dy;
  Vzz += dz*dz;
  Vxy += dx*dy;
  Vyz += dy*dz;
  Vzx += dz*dx;
}

その他の回答(1件)

id:pyopyopyo No.1

pyopyopyo回答回数348ベストアンサー獲得回数862009/09/12 11:30:08

ポイント24pt

要は、最小2乗法のオペレーションが判らない、ということですね。

以下のページに、具体的な説明と、数字を使った例があります。

http://edu.isc.chubu.ac.jp/hsuzuki/iip/2007-katsuyou/w5/lms2.htm...

この例では、2次元の点群から、その点群にフィットする直線を求めています。

spin6536 さんがご覧になっている英語の例では、これと同じオペレーションで、3次元の点群から、その点群にフィットする平面を求めることになります。

つまり

Now take the partial derivative of g with respect to a and set it

equal to zero. Do the same for b and c. This will give you three

linear equations to solve for a, b and c.

ということです。

id:spin6536

教えていただいたサイトの説明は非常に細かくてよかったです。でも、自分の数学力不足のため、三次元に応用して、実際にどのようにやったらよいのかがわからないのです。あと、細かな基本的と思われるところも。偏微分(かな?)の式をどのように立てて、それをどのように解いて、といった話です(さすがに式があるていど単純ならそれを微分することぐらいはできますが)。そこで、デモンストレーションをここでお願いしたのです。

2009/09/12 11:51:40
id:ita No.2

ita回答回数204ベストアンサー獲得回数482009/09/13 19:59:14ここでベストアンサー

ポイント46pt

ーーー計算のバックグラウンドーーー

i番目のデータの実測値と推定値のズレ Eiは

Ei = a*(xi-x0) + b*(yi-y0) + z0 -zi

これの二乗和 E=ΣEi^2 を最小にするaとbをもとめます。それにはEをaで微分して0かつEをbで微分して0

という条件を満たすa、bを求めます。

aで微分すると (d/da) ΣEi^2 = 2ΣEi (d/da)Ei = 2Σ (a*(xi-x0)+b*(yi-y0) + z0 - zi)*(xi-x0)

=2aΣ(xi-x0)(xi-x0) +2bΣ(yi-y0)(xi-x0)+2Σ(z0-zi)(xi-x0)

Σ(xi-x0)(yi-y0)=Vxyなどとおくと

=2a Vxx + 2b Vxy -2Vzx =0 これが第一の式

bで微分すると同様な計算で

2a Vxy + 2b Vyy - 2Vyz = 0

行列で書くと

 |Vxx  Vxy | |a|   |Vzx|
 |Vxy  Vyy | |b| = |Vzy|
行列の逆行列を|Vzx Vzy|にかけるとa,bが出る。逆行列は
   1              |  Vyy  -Vxy |
 ---------------- |            |
 Vxx Vyy - Vxy^2  |  -Vxy  Vxx |

したがって最終的に

a= (Vyy*Vzx - Vxy*Vzy)/(Vxx*Vyy-Vxy*Vxy)

b= (-Vxy*Vzx +Vxx*Vzy)/(Vxx*Vyy-Vxy*Vxy)

Vxxなどは以下のように計算

double x0=0, y0=0, z0=0
for(i=0;i<N;i++){
x0 += x[i];
y0 += y[i];
z0 += z[i];
}
x0/=N;
y0/=N;
z0/=N;

double Vxx=0, Vyy=0, Vzz=0, Vxy=0, Vyz=0, Vzx=0;
for(i=0;i<N;i++)
{
  double dx = x[i]-x0;
  double dy = y[i]-y0;
  double dz = z[i]-z0;

  Vxx += dx*dx;
  Vyy += dy*dy;
  Vzz += dz*dz;
  Vxy += dx*dy;
  Vyz += dy*dz;
  Vzx += dz*dx;
}

コメントはまだありません

この質問への反応(ブックマークコメント)

「あの人に答えてほしい」「この質問はあの人が答えられそう」というときに、回答リクエストを送ってみてましょう。

これ以上回答リクエストを送信することはできません。制限について

絞り込み :
はてなココの「ともだち」を表示します。
回答リクエストを送信したユーザーはいません