円の最小二乗法のプログラム

今、研究でVBで円の最小二乗法のプログラムを使いたいのですが全然わかりません。誰か詳しい人お願いします。

No.4 回答者： usokoku 回答日時： 2008/10/16 22:53

一番です。

今帰ってきたと子なので大雑把に。
鞍部問題でいやってほど苦しんだ人間ですので、尋常な方法ではありません。

点(Xi, Yi), i=1,2,...,n を通る円の方程式を求めるとして、
円の中心点(X0, Y0)の座標は
Xiの最大値Xmaxと最小値Xminの間に X0が存在する
Yiの最大値Ymaxと最小値Yminの間に Y0が存在する
ことが「偏りがないサンプリング」で、十分大きな資料数がある場合にはは言えるでしょう(私が解いた問題では、偏りがあって苦労した)。
XmaxとXminの間を16分割して(Xmax-Xmin)/128の幅(Xhaba)で
YmaxとYminの間を16分割して(Ymax-Ymin)/128の幅(Yhaba)で
16*16 =256 個の点のどこかに中心点があるとして、すべての観測点との距離を計算して、距離の平均と、計算でもとのた距離との差の絶対値が最小となった点(Xk, Yk)をもとめます(私の場合には、Yk=Ymaxとなって本来求めたいYk=Ymin側の中心点が消えてしまった)。

点(Xk, Yk)の周り , Xk-Xhaba から Xk+Xhabaの間, Yk-Yhaba から Yk-Yhaba の間について、同様に、X方向・y方向ともに16分割して計算して行きます。
これを何回か繰り返してゆくと、残差の変化がほとんどない点が最適解としてえられます。

これで求められた点と、観測点の重心を中点とする反対側に鞍部の向こう側の解となる中心点(Xl, Yl)があると仮定して
点(Xl, Yl)の周り , Xl-Xhaba から Xl+Xhabaの間, Yl-Yhaba から Yl-Yhaba の間について、同様に、X方向・y方向ともに16分割して計算して行きます。
これを何回か繰り返してゆくと、残差の変化がほとんどない点が最適解としてえられます。
以上の方法て゛、2つの解がえられます。定義域などの情報から一方の解は捨てられるでしょう。

個の方法は、計算時間がやたらかかりますが、鞍部問題に引っかかりにくいので便利です。点(1.1,1.1), (1,-1), (-1,1), (-1,-1)の４点が与えられているときに、
点(0,0)付近を中心とした半径1.4位の円を想定するか、
点(10,10)(くらいかな)付近を中心とする半径10くらいの円を想定するか、
点(-10,-10)(くらいかな)付近を中心とする半径10くらいの円を想定するか、
です。

手作業で行うのであれば、任意の点３点を選んで、３点を通る円の式を解析で求めて、これを初期値(Xm, Ym)とする。
点(Xm, Ym)の周り , Xm-Xhaba から Xm+Xhabaの間, Ym-Yhaba から Ym-Yhaba の間について、同様に、X方向・y方向ともに16分割して計算して行きます。
という方法もあります。「任意の点」というよりも、点間距離が大きくなる点を選んだ(総当りで３てん感の距離の合計を求めて最大となった点を初期値とする。ただし正三角形に近い形の場合が収束が早い)ほうが、最適解になりやすいです。

方法はいくつかありますので、プログラミングの技術(含プログラム作成時間がどのくらいさけるか、最初の方法ですとFor-Nextの繰り返しだけだから３－４日で何とか形になる、を含む)と相談して選んでください。「研究する」ことが重要であって、プログラミングは研究するための道具でしかないのですから、多少計算時間はかかっても(夜間放置して解を得ることになっても)、研究のほうが重要ですので、プログラミングは短時間で終了させる必要があります。

No.3 回答者： ShowMeHow 回答日時： 2008/10/16 18:31

間違えました。

ＡＮｏ２です。
http://oshiete1.goo.ne.jp/qa3712186.html
の2番目の回答でした。
申し訳ありません。

No.2 回答者： ShowMeHow 回答日時： 2008/10/16 16:31

全く詳しくないド素人（プログラムも数学も）ですが、

プログラムの手順としては
1. Xj,Yjを読み込む。　ファイル形式で読み込むのか、数値入力するのか、ビジュアルに入れるのかその辺はきまっていますか？
2. x0,y0を計算で出す。
ヒントとしては、
S=Σ{(xj-x0)^2 +(yj-y0)^2}
=Σ(xj-x0)^2 +Σ(yj-y0)^2
=x1＾2-2x1x0+x0＾2+x2＾2-2x2x0+x0＾2+.....+Σ(yj-y0)^2
=Σxj＾2-Σ2xjx0+Σx0＾2+Σ(yj-y0)^2
∂S/∂x0=-Σ(2xj-2x0)
⇒Σxj=n*x0
x0 = Σxj/n
ということでちょっと気になっている、http://oshiete1.goo.ne.jp/qa4381786.htmlの2番目の回答と同じ。
これがどういうことかというと、読み込んだxjを全部足して読み込んだ個数で割るという意味。
3.pを求める。
4.結果を表示する。

どう見ても研究というより課題だけど、自分でやらないと意味がないと思う。

No.1 回答者： usokoku 回答日時： 2008/10/15 19:24

「円の」の意味が解釈不能。

VBプログラムを作るとして、求める方程式と、どの残差を最小にしたいか、

それとも、どなたかが作成したプログラムでしょうか、こちらですと、作成者に聞いたほうが早いのですが。

この回答へのお礼

ここでの円は座標上での円です。

残差を最小にして計算を行いたいです。

S=Σ{(xj-x0)^2 +(yj-y0)^2}
として、
∂S/∂x0=0
∂S/∂y0=0
より、(x0,y0)が求め、次に半径r0を
ρ=(1/N)Σ√{(xj-x0)^2 +(yj-y0)^2}
として、
r0＝pを求めるプログラムを作りたいのです。
ちなみにどなたかがこのプログラムを作成していないか探したところ見つかりませんでした。

内容不足ですみませんでした。

お礼日時：2008/10/16 11:54