logの入った最小二乗法について

質問

Y=ln(x-b)-ln(x-c)
といったように、ｘが二つに分かれてしまっていて、さらに式にあるようにｂやｃが入っているときの最小二乗法のやり方を教えてください。お願いします。

回答No.3

stomachman

No.1の補足を見ると、どうやらstomachman流「超平面法」の出番のようですね。

(1)　lnY(lnA) = aln(x-b) - bln(x-c)
左辺は(ln(Y))×(ln(A))の意味でしょうか、あるいはln(Y×(ln(A)))でしょうか。どちらにせよ、Aは既知らしいですから、まとめて
(2)　ln(y) = lnY(lnA)
と書くことにしましょう。
それより、bが２箇所に現れる。分かりにくくなっちゃうので、とりあえず、
(3)　ln(y) =aln(x-b) + dln(x-c)、ただしd=-b
ということにしましょう。

こいつの両辺をxで微分します。ここで
(d/dx)ln(y)=(dy/dx)/y
(d/dx)aln(x-b)=a/(x-b)
ですから、微分方程式
(4)　(x-b)(x-c)(dy/dx) =((x-c)a + (x-b)d)y
が得られる。x,yの項別に整理すると
(5)　(x^2)(dy/dx) -(b+c)x(dy/dx) + bc(dy/dx) = (a+d)xy - (ac+bd)y
この両辺をxで積分するんです。部分積分を使って
∫(x^n)(dy/dx)dx = (x^n)y - n∫(x^(n-1))y dx
ですから、未知の定数項tを入れて
(5)　(x^2)y- 2∫xy dx -(b+c)(xy -∫y dx) + (bc)y = (a+d)∫xy dx - (ac+bd)∫y dx + t
となる。これを整理すると
(6)　(x^2)y = (b+c)(xy -∫y dx) + (a+d+2)∫xy dx - (ac+bd+b+c)∫y dx - (bc)y + t
です。ここで
　p=(b+c)
　q=(a+d+2)
　r=-(ac+bd+b+c)
　s=-bc
と定義すれば、p,q,r,s,tの５つのパラメータを持つ線形モデルになります。
積分はデータx[k]が得られている範囲について、数値積分でやるんです。つまり、データ(x[k],y[k]) (k=0,1,2,...,K)を、x[k]が小さい順に並ぶようにsortしておきます。
たとえば台形則で積分をやるものとすれば、
　P[k] = x[k]y[k] -R[k]
　Q[k] = Σ(y[j]x[j]+y[j-1]x[j-1])(x[j]-x[j-1])/2
　R[k] = Σ(y[j]+y[j-1])(x[j]-x[j-1])/2
　S[k] = y[k]
　T[k] = (x[k]^2)y[k]
として（Σはいずれもj=1～kについて取ります。）、モデルは
(7)　T[k] = pP[k]+qQ[k]+rR[k]+sS[k]+t
です。p,q,r,s,tに関する線形モデルですね。

さて、このモデルをデータ(x[k],y[k])にフィッティングするとp,q,r,s,tの５個のパラメータが出てきてしまう。もともとa,b,cの３個しかパラメータがないはずなのに、これじゃ過剰です。過剰ですが、データがモデルと良く合っている（残差がごく小さくなるa,b,cが存在する）ならば、p,q,r,sのどの組み合わせを使って求めたa,b,cもほぼ同じになる筈です。
　なんでこんな過剰が出たかというと、tは積分定数ですから上記の処方で出てきた余計なモノで、単に無視して良し。そしてp,q,r,sはa,b,c,dの４パラメータを含むモデルなら丁度良い訳です。

　以上のようにして、a,b,cを決めたら、さらに非線形最小二乗法で改良することもできます。a,b,cが既に「正解」にごく近い近似値になっているので、ガウス・ニュートン法など、線形近似法を反復させて簡単に収束させられます。

　なお、stomachmanは計算間違い、書き間違いの常習犯です。チェック宜しく。