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

このQ&Aに関連する最新のQ&A

A 回答 (4件)

No.3の補足説明です。


lnY(lnA) = aln(x-b) - dln(x-c)
において
(2) ln(y) = lnY(lnA)
とする代わりに、
(2') y=lnY(lnA)
とした方が易しいですね。こうすれば微分方程式は
(4') (x-b)(x-c)(dy/dx) =((x-c)a + (x-b)d)
となり、
(5') (x^2)(dy/dx) -(b+c)x(dy/dx) + bc(dy/dx) = (a+d)x - (ac+bd)
積分すると
 (x^2)y- 2∫xy dx -(b+c)(xy -∫y dx) + (bc)y = ((a+d)/2)(x^2)- (ac+bd)x + t
より
(6') (x^2)y- 2∫xy dx = (b+c)(xy -∫y dx) +((a+d)/2)(x^2)- (ac+bd)x - (bc)y + t
 p=b+c
 q=(a+d)/2
 r=-(ac+bd)
 s=-(bc)
とおいて、
 P[k] = x[k]y[k] -Σ(y[j]+y[j-1])(x[j]-x[j-1])/2
 Q[k] = x[k]^2
 R[k] = x[k]
 S[k] = y[k]
 T[k] = (x[k]^2)y[k] - Σ(y[j]x[j]+y[j-1]x[j-1])(x[j]-x[j-1])
として
(7) T[k] = pP[k]+qQ[k]+rR[k]+sS[k]+t
の方が多少すっきりしてますね。いや、やってることは同じですけど。
    • good
    • 0

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が2箇所に現れる。分かりにくくなっちゃうので、とりあえず、
(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の5つのパラメータを持つ線形モデルになります。
積分はデータ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の5個のパラメータが出てきてしまう。もともとa,b,cの3個しかパラメータがないはずなのに、これじゃ過剰です。過剰ですが、データがモデルと良く合っている(残差がごく小さくなるa,b,cが存在する)ならば、p,q,r,sのどの組み合わせを使って求めたa,b,cもほぼ同じになる筈です。
 なんでこんな過剰が出たかというと、tは積分定数ですから上記の処方で出てきた余計なモノで、単に無視して良し。そしてp,q,r,sはa,b,c,dの4パラメータを含むモデルなら丁度良い訳です。

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

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

●正攻法だと非線形最小二乗法を使います。

過去の質問から「最小二乗法」を検索すれば、参考資料等いろいろ分かるはず。

●最小二乗法の一般的注意事項として、何を最小にするかをはっきりさせる必要があります。モデル
y=ln(x-b)-ln(x-c)

exp(y) = (x-b)/(x-c)

(x-c) exp(y)=x-b
と同じ意味ですけれど、
残差
ε[k] = y[k] - (ln(x[k]-b)-ln(x[k]-c))
の二乗和を最小にするb, cは、残差
δ[k] = exp(y[k]) - (x[k]-b)/(x[k]-c)
の二乗和を最小にするb, cや
γ[k] = (x[k]-c) exp(y[k])-(x[k]-b)
の二乗和を最小にするb, cとは違いますからね。
 しかし、データ(x[k],y[k]) (k=1,2,.....,K) がモデルと良く合っていて残差がごく少なくできる(そのようなb,cが存在する)場合には、どれでやってもほとんど同じb,cが得られます。そういう時は、実用上どれでやっても大差ない訳です。

●以上をふまえた上で、ご質問の例の場合、モデルは
x[k] (exp(y[k])-1)=c exp(y[k])-b
と表すことができます。だから、c,bについて見ると線形最小二乗法になってます。どういうことかというと、
X[k]=exp(y[k])
Y[k]=x[k] (exp(y[k])-1)
とおけば、X[k],Y[k]は共にデータ(x[k],y[k])が与えられれば値が確定する。そしてモデル
Y[k]=c X[k] - b
のb,cを決めろ、という問題です。これなら簡単でしょ?

このやり方は、飽くまでも残差
β[k] = x[k] (exp(y[k])-1) - (c exp(y[k])-b)
の二乗和を最小にするb, cを求めるのであって、ε[k]の二乗和を最小にするb, cとは完全には一致しない。これは頭に入れて置いてくださいね。

と回答しようとしたら、siegmund先生が既に......なはは。せっかく書いたからupしちゃおうっと。
    • good
    • 0

(1)  Y = ln{(x-b)/(x-c)}


として,
(2)  e^Y = (x-b)/(x-c)  ⇒  (x-c)e^Y = x-b
   ⇒  cz - b - xz + x = 0
と変形すれば(z = e^Y) ,定めるべき係数 b,c について線型になりますから
普通のよく本に載っている方法がそのまま適用できます.

細かいことを言えば,Y について 測定値がガウス分布することと,
e^Y についてガウス分布することとは違うのですが,
通常余り気にしません.

この回答への補足

 丁寧に答えていただき、ありがとうございます。
教えていただいた方法で解いてみようと思ったのですが、
問題のタイプが少し違うようでした。
 私がいま直面している問題は
    lnY(lnA) = aln(x-b) - bln(x-c) (Aは定数です。)
という形だたのですが、私が勝手に簡略化してしまっていました。
すいません。
これをどうしてもよい形にもっていくことができません。
もしよい方法をご存知でしたら、お教え頂きたいと思います。
よろしくお願いします。

      

補足日時:2001/06/29 21:18
    • good
    • 0

このQ&Aに関連する人気のQ&A

お探しのQ&Aが見つからない時は、教えて!gooで質問しましょう!

このQ&Aを見た人が検索しているワード

このQ&Aと関連する良く見られている質問

Q最小二乗法での指数関数の計算

最小二乗法での指数関数の計算
最小二乗法での指数関数の計算方法が良く分からないのですが公式などありますか?
y=ae^bxでしたらやり方があるのですがy=ae^bx+cの方法がわかりません・・・・・・

Aベストアンサー

No1の方が書いておられるように、これは厳密には非線形の最小二乗法になります。
ただし、単にそう書いても、質問者の方にはわからないと思いますので、具体的かつ実用的なやり方を書きましょう。
と言っても、簡単です。
まず、cとして、適当な値c1を仮定します。
すると、
y-c1=a1e^b1x
の最小二乗法に帰着され、a1,b1が求まります。
この時の残差を、z1とします。
z1=Σ(yi-c1-a1e^b1xi)^2

次に、cの増分値をΔcとして、
c2←c1+Δc
と置いて、また最小二乗法を適用し、残差z2を求めます。
z2=Σ(yi-c2-a2e^b2xi)^2

さらに、
c3←c1+2Δc
と置いて、また最小二乗法を適用し、残差z3を求めます。

この段階で、cに対して残差zをプロットし、2次関数近似してみましょう。
うまくいけば、極小値がこの範囲内にみつかります。
この2次関数から直接、または、初期値c1と増分値Δcをもっと適切な値となるように設定するなりして極小値の存在範囲を狭めてから2次関数近似して、この極小値cを求めましょう。
最後にもういちど、cに対する最小二乗法を適用して、係数a,bを求めれば、それが答です。

もし、cに対するzが、単調増加、あるいは単調減少になってしまったら?
この場合、Δcが大きすぎると、極小値が範囲内にあるにも関わらず、単調増加あるいは単調減少になってしまっている可能性も高いので、c1,c2,c3は変えず、Δcを1/2にして、c1とc2、c2とc3の間の値を求め、全5点をプロットしてみます。
本質的に単調増加あるいは単調減少なら、傾向に変化はないはずです。この場合には、zが小さくなる方向にcの初期値と増分値Δcを選びなおして、再計算してみましょう。
区間内に極小値がある場合には、傾向が”それらしく"変化しますから、増分値をもっと細かくして、極小値を押さえれば良いのです。

慣れてくれば、2次関数近似といわず、最初から絨毯爆撃的にcを規則的に変化させながらzの動きを把握していくなどの小ざかしい方法を覚えるようになりますが、最後は極小値の存在する区間の3個の値から2次関数近似することには変わりがありません。(2次関数近似は扱いが簡単で便利ですから。)

No1の方が書いておられるように、これは厳密には非線形の最小二乗法になります。
ただし、単にそう書いても、質問者の方にはわからないと思いますので、具体的かつ実用的なやり方を書きましょう。
と言っても、簡単です。
まず、cとして、適当な値c1を仮定します。
すると、
y-c1=a1e^b1x
の最小二乗法に帰着され、a1,b1が求まります。
この時の残差を、z1とします。
z1=Σ(yi-c1-a1e^b1xi)^2

次に、cの増分値をΔcとして、
c2←c1+Δc
と置いて、また最小二乗法を適用し、残差z2を求めます。
z2=Σ(yi-c2-a2e^b2xi)...続きを読む

Q誤差の二乗を最小にする理由

収集したデータをある関数でフィッティングする際、収集したデータと関数の差を二乗した合計が最小になるよう、関数を求める方法がありますが、なぜ二乗なのでしょうか。
統計的な根拠があるという話を聞いたのですが、WEBで検索しても手法の説明や実際の計算の仕方ばかり検索され、根拠がなかなかみつかりません。
なぜ、絶対値の合計や3乗、4乗、平方根ではなく、二乗の和を使用するのでしょうか。

Aベストアンサー

[1] 厳密な話ならばANo.1の通りです。

 最小二乗法が厳密な意味でモデルのパラメータの最尤推定法になるためには、以下の前提が必要です。
(1) 誤差、すなわち測定値y[i]と真値y0[i]の差
  ε[i] = (y[i]-y0[i])
は平均0、分散(σ[i])^2の正規分布 N(0, (σ[i])^2)に従う。ただし、
  σ[i] = s[i] σ0
で、σ0は未知でも良いが、s[i]は既知である。
(2) 各測定は互いに独立で、それぞれの誤差の共分散は0である。
(3) パラメータのベクトルxを含む既知のモデルM(x,i)があって、真値y0[i]について、あるパラメータx0が存在して、
  M(x0,i)=y0[i]
である。

 つまり、正解のx0を使って計算した残差y[i]-M(x0,i)は、正規分布 N(0, (σ[i])^2)に従う独立なランダム変数n個それぞれから取ったサンプルになっている、ということです。

 さて、真値の推定値としてモデルが与える値 M(x,i) を使ったとき、その尤度(すなわち、測定値がy[i]であるときに真値がM(x,i)である確率)は、上記(1)(2)の仮定から、
  L(M(x,i) | y[i]) = {1/(σ[i]√(2π))} exp( -((y[i]-M(x,i))^2) / (2(σ[i])^2) )
である。xの尤度L(x | y[i](i=1,...,n))(すなわち、測定値がy[i](i=1,...,n)であるときにパラメータの正解がxである確率)は、i=1...nについてのL(M(x,i) | y[i]) の積
  L(x | y[i](i=1,...,n)) = Π L(M(x,i) | y[i])
   = [1/{((2π)^(n/2))Πσ[i]} ] exp[-E(x)/(2(σ0^2))]
である。ただし
E(x) = Σ{((y[i]-M(x,i))/s[i])^2}
としました。
 確率の意味で最も尤もらしいx0の推定値は、尤度L(x | y[i](i=1,...,n)) を最大にするようなxである(最尤推定)。そして、L(x | y[i](i=1,...,n)) の式から明らかなように
  L(x | y[i](i=1,...,n))を最大にするようなx ⇔ E(x)を最小にするx
です。

[2] しかしながら、実務においては「必ずしも最尤推定だと保証できなくてもいい」という場合が多々ある。また、「データy[i]が含む測定誤差が正規分布に従い、その分散の相対値s[i]が分かっている」という条件や「真値は(パラメータさえ正しければ)モデルで誤差なしに説明できる」という条件を満たせない場合も多い。なので、実務のセンスで言いますと、ANo.2の説明もまた適切であろうと思います。

[1] 厳密な話ならばANo.1の通りです。

 最小二乗法が厳密な意味でモデルのパラメータの最尤推定法になるためには、以下の前提が必要です。
(1) 誤差、すなわち測定値y[i]と真値y0[i]の差
  ε[i] = (y[i]-y0[i])
は平均0、分散(σ[i])^2の正規分布 N(0, (σ[i])^2)に従う。ただし、
  σ[i] = s[i] σ0
で、σ0は未知でも良いが、s[i]は既知である。
(2) 各測定は互いに独立で、それぞれの誤差の共分散は0である。
(3) パラメータのベクトルxを含む既知のモデルM(x,i)があって、真値y0[i]について、あるパラメー...続きを読む


人気Q&Aランキング