最小二乗法 楕円放物面

数千のｘｙｚ座標データを、最小二乗法を用いて、ax^2+by^2+2cz=1の楕円放物面に近似したいのですが、どのようにしたらいいですか？
ご存知の方教えてください。
宜しくお願いします。

No.5 回答者： mocapapa89 回答日時： 2008/11/10 11:49

No.3です。

No.1で例示されたURLには線形最小二乗法とかかれており、モデル関数が線形の場合は連立方程式を解けば求まりますが、本問は非線形なので数値的に解くしかないと思います。私がNo.3でa, b, cについて線形と申し上げたのは誤りなので訂正させてください。lp_solveは使えません。ちなみに経験者も一般人の誤りでした(寝ぼけてたのかな)。

さて、非線形のカーブフィッティングはoctaveでもできるかとは思いますが、私もここで質問して(答えは得られませんでしたが)調べたところwnlibというAPIで解くことができました。私の場合はlogistic curveへのフィッティングでしたが、ご参考まで。

参考URL：http://www.willnaylor.com/mantext/wnnlp.txt

No.4 回答者： echoes_x86 回答日時： 2008/11/07 09:09

おはようございます．

まだ閉じていなかったのですね．

さて，k番目のデータの真値(xk,yk,zk)が(xk,yk,zk+dk)と観測されれば，
残差は次のように書けます．
dk=zk-(axk^2+byk^2+1)/2c ...(1)
これはzのみに誤差がのっている場合に対応します．
とりあえず，これの平方和を最小化する方法を，と言うことであれば，
cの逆数を定義するか何かすれば，
非負制約つき二次計画(QP)辺りを利用できるかもしれません．
線型計画(LP)でないのは「二乗誤差」の最小化が目的だからです．
裸のoctaveでは結構難しいと思います．
この意味ではExcelソルバのGRG2は最も良い選択肢の一つでしょう．

しかしながら，z以外にも誤差がある場合は事情が異なります．
これは(xk+dxk, yk+dyk, zk+dzk)を元の式に代入すればすぐに分かりますが，
式(1)のように誤差関数を簡単には書けません．
したがって，パラメターの現在値で書ける曲面上にあって，
観測データに最も近い点を(場合によっては数値的に)
計算する必要が生じるわけです(直交射影)．
誤差の分布によっては「評価関数の計算」からして骨が折れますが，
これが計算できなければどんな方法でもパラメターは更新できません．
Excelソルバも同様でしょう．
だから必要以上に問題が難しくならないように「補足要求」をしたのです．
参考までに述べると「最も一般的に」考えると，
修士以上の学位論文が執筆可能な内容になると思われます．

No.3 回答者： mocapapa89 回答日時： 2008/11/06 18:10

普通なら(線形だろうが非線形だろうが)エクセルソルバで最小二乗法は解くのですが、数千となるとだめなのでしょうか。

となるとパラメータa, b, cについて線形なのでlp_solveはだめでしょうか？

http://lpsolve.sourceforge.net/5.5/

あるいはoctaveを使うとか。

http://ja.wikipedia.org/wiki/GNU_Octave

No.2 回答者： echoes_x86 回答日時： 2008/10/17 15:49

こんにちは．

データの分布と誤差の乗り方によってはかなり難しくなりますよ．
というのも，ax^2+by^2+2cz=1は楕円の「標準型」で，
適切な直交変換をかけてxy平面の断面が楕円に，長軸と短軸がx軸とy軸になるようにしたものだからです．
入力する観測データは上記のことが期待できるのでしょうか？
また，誤差の分布はどのようになっていますか？
観測データの真値(x^,y^,z^)が楕円放物面に拘束されているとして，
観測データ(x,y,z)が(x^+u,y^+v,z^+w)という誤差を含んで観測されるとなると，
パラメターの現在値(a,b,c)で書ける楕円放物面への直交射影を計算せねばなりません．
さらにパラメターはa>0，b>0なる非負制約のもとで更新する必要があります
(単純な最小二乗法ではa>0，b>0となる保証がなく，「楕円放物面」以外の二次曲面が当てはまる恐れがあります)．
これらの情報が不足しています．問題設定をもう少しきちんとした方が回答もつきやすいと思います．

No.1 回答者： amue 回答日時： 2008/10/17 15:25

最小二乗法とは残差二乗和を各係数で偏微分した式を=0とした連立方程式を立てて解くことですので、

残差二乗和をS=Σ(ax^2+by^2+2cz-1)^2　として、
各係数a,b,c毎に偏微分した式で連立方程式を立てると

(Σx^4)a+(Σx^2 y^2)b+2(Σx^2 z)c = Σx^2
(Σx^2 y^2)a+(Σy^4)b+2(Σy^2 z)c = Σy^2
(Σx^2 z)a+(Σy^2 z)b+2(Σz^2)c = Σz
となります。
この連立方程式を解けば係数a,b,cが求まります。

参考URL：http://www.eli.hokkai-s-u.ac.jp/~kikuchi/ma2/cha …

この回答への補足

回答ありがとうございます。

例えば（x1,y1,z1)における誤差dは
d=z1-(aXk^2+bYk^2+1)になるのではないのですか？

ようするに
残和二乗和Σd=Σ{Zk-(aXk^2+bYk^2+1)}になると思うのですが、
どうでしょうか？
ご意見お願いします。

補足日時：2008/10/17 16:03