y=a/(x-b)+cの最小二乗法

y=a/(x-b)+cの最小二乗法

y=a/(x-b)+c
という、反比例の式をｘ方向に＋ｂ、ｙ方向に＋ｃ平行移動したような曲線の係数a,b,cを求めるための最小二乗法の方法を教えていただけないでしょうか。

工夫してみたのですが、なかなかうまくいきませんでした。
すみませんが、力を貸してください。

No.2 ベストアンサー 回答者： stomachman 回答日時： 2009/12/30 07:08

　データ(x[k], y[k]) (k=1,2,...,N)

が与えられた時に、残差(residue)
r[k] = a/(x[k]-b)+c-y[k]
の二乗和(square sum of residue)
E(a,b,c) = Σ(r[k]^2) （Σはk=1,...,Nについての総和）
を最小にする係数(a,b,c)を求む、というのが「y=a/(x-b)+cの最小二乗法」の意味するところですね。

　Eが極値を取る係数(a,b,c)を知るために「Eをa,b,cそれぞれで偏微分したものが0になる」という条件を表す連立方程式
∂E/∂a = 0
∂E/∂b = 0
∂E/∂c = 0
を作ってみますと、左辺を計算すると未知数a,b,cについて一次方程式にはならず、これじゃ解けない。

　連立一次方程式にならない場合を「非線形最小二乗法」と（そして、連立一次方程式になってくれる場合を「線形最小二乗法」と）呼びます。非線形最小二乗法では、適当な出発値(a[0], b[0], c[0])から始めて、Eが小さくなるようにちょっとずつ(a,b,c)を改良して行く繰り返し計算が必要です。
　改良のための計算方法はいろいろ工夫（Gauss-Newton法、最急降下法、Marquardt法など）されていますが、最も単純でオバカなやり方は、
(1) b,cを定数、aだけを変数と考えてEを最小化する。
(2) a,cを定数、bだけを変数と考えてEを最小化する。
(3) a,bを定数、cだけを変数と考えてEを最小化する。
を順繰りに何度も何度も繰り返すことです。（教科書は「最小二乗法による実験データの解析」（東京大学出版会）をお薦めします。）

　案外重要なのが出発値(a[0], b[0], c[0])の選び方。ここで言う「適当な」とは、テキトーじゃ駄目なんで、（本来の意味での）適当な、つまり真の解に十分近い出発値を選ばないと改良の計算が収束しないことがあります。
　出発値を選ぶひとつの方法は、「もし残差が0になるようなデータが与えられたならばEに一致して、しかも簡単に計算できる」ような値Fを定義して、Fを最小化する(a,b,c)を(a[0], b[0], c[0])として利用する、というやりかたです。具体的には、
r[k] = a[0]/(x[k]-b[0])+c[0]-y[k]
を
r[k](x[k]-b[0]) = a[0]+(c[0]-y[k])(x[k]-b[0])
と変形しておいて右辺を展開すると、
r[k](x[k]-b[0]) = (a[0]-b[0]c[0])+b[0]y[k]+c[0]x[k]-y[k]x[k]
ここで
A = a[0]-b[0]c[0]
B = b[0]
C = c[0]
とおくと、
r[k](x[k]-b[0]) = A+By[k]+Cx[k]-y[k]x[k]
と書けますから、この右辺を
s[k] = A+By[k]+Cx[k]-y[k]x[k]
と定義して
F(A,B,C) = Σ(s[k]^2) （Σはk=1,...,Nについての総和）
を最小化するA,B,Cを考える。
すると、s[k]はA,B,Cの一次式になっているので、「普通の（線形）」最小二乗法でA,B,Cが決定できます。そしてA,B,Cからa[0],b[0],c[0]が決まる。（もしr[k]が全て0であれば、E(a[0],b[0],c[0])=0 になるのは自明でしょう。）
　こうして得たa[0],b[0],c[0]はEを最小化しません（∵r[k]が全て0にならない限り、E≠F）が、もしも「Fが十分小さくできるような旨いa,b,cが存在する」のであれば、そのa,b,cはEも小さくしますから、出発値に利用できる訳です。
　言い換えると、もしEが十分小さくなるようなa,b,cが存在する（つまり、モデル y=a/(x-b)+cがデータに良く合う）のであれば、Fを最小化するa,b,cと、Eを最小化するa,b,cとは非常に近い値になる。
　なので実用上は、(a[0],b[0],c[0])だけで十分な精度が得られる（改良のための繰り返し計算をしないで済んでしまう）ことも多々あります。

この回答へのお礼

ありがとうございます！
助かりました。
やはりニュートン法などの方法をとるのがベストですよね。
初期値の与え方、すごく勉強になりました。
細かい説明と式展開、非常にわかりやすく、すぐに実践できそうです。

本当にありがとうございました。

お礼日時：2010/01/01 00:15

No.1 回答者： sanori 回答日時： 2009/12/29 21:37

こんばんは。

反比例の最小二乗の経験はないですが、やってみます。


あるデータ（ｘk，ｙk）に対して、誤差εk を
εk　＝　ａ／（ｘk － ｂ）　＋　ｃ　－　ｙk
と定義します。

二乗誤差の合計Ｓは、
Ｓ　＝　Σ[k=1→n] εk^2
　＝　Σ[k=1→n]｛ａ／（ｘk － ｂ）　＋　ｃ　－　ｙk ｝^2
ですので、これを極小にすることを方針とすればよいわけです。

ここで、
εk^2　＝　｛ａ／（ｘk － ｂ）　＋　ｃ　－　ｙk ｝^2
　＝　ａ^2／（ｘk － ｂ）^2　＋　ｃ^2　－　ｙk^2　＋　２ａｃ／（ｘk － ｂ）　－　２ａｙk／（ｘk － ｂ）　－　２ｃｙk
であり、これをａ、ｂ、ｃのそれぞれで偏微分すれば、

∂εk^2／∂ａ　＝　２ａ／（ｘk － ｂ）^2　＋　２ｃ／（ｘk － ｂ）　－　２ｙk／（ｘk － ｂ）　　・・・（あ）

∂εk^2／∂ｂ　＝　－２ａ^2／（ｘk － ｂ）^3　－　２ａｃ／（ｘk － ｂ）^2　＋　２ａｙk／（ｘk － ｂ）^2　　・・・（い）

∂εk^2／∂ｃ　＝　２ｃ　＋　２ａ／（ｘk － ｂ）　－　２ｙk　　・・・（う）

よって、

∂Ｓ／∂ａ　＝　Σ[k=1→n]｛２ａ／（ｘk － ｂ）^2　＋　２ｃ／（ｘk － ｂ）　－　２ｙk／（ｘk － ｂ）｝　　・・・（あ’）

∂Ｓ／∂ｂ　＝　Σ[k=1→n]｛－２ａ^2／（ｘk － ｂ）^3　－　２ａｃ／（ｘk － ｂ）^2　＋　２ａｙk／（ｘk － ｂ）^2｝　　・・・（い’）

∂Ｓ／∂ｃ　＝　Σ[k=1→n]｛２ｃ　＋　２ａ／（ｘk － ｂ）　－　２ｙk｝　　・・・（う’）

最小二乗ということは、Ｓが極小ということなので、
∂Ｓ／∂ａ＝０　かつ　∂Ｓ／∂ｂ＝０　かつ　∂Ｓ／∂ｃ＝０

つまり、
０　＝　Σ[k=1→n]｛２ａ／（ｘk － ｂ）^2　＋　２ｃ／（ｘk － ｂ）　－　２ｙk／（ｘk － ｂ）｝　　・・・（あ’’）

０　＝　Σ[k=1→n]｛－２ａ^2／（ｘk － ｂ）^3　－　２ａｃ／（ｘk － ｂ）^2　＋　２ａｙk／（ｘk － ｂ）^2｝　　・・・（い’’）

０　＝　Σ[k=1→n]｛２ｃ　＋　２ａ／（ｘk － ｂ）　－　２ｙk｝　　・・・（う’’）


式を簡単にすると、
０　＝　Σ[k=1→n]｛ａ　＋　ｃ（ｘk － ｂ）　－　ｙk（ｘk － ｂ）｝　　・・・（あ’’’）

０　＝　Σ[k=1→n]｛ａ　＋　ｃ（ｘk － ｂ）　－　ｙk（ｘk － ｂ）｝　　・・・（い’’’）

０　＝　Σ[k=1→n]｛ａ　＋　ｃ（ｘk － ｂ）　－　ｙk（ｘk － ｂ）｝　　・・・（う’’’）

ここまで来て気づきましたが、３つの式は同じですね。
１本の式で３つの未知数を求める連立方程式はありませんから、
普通の最小二乗法はできないということです。