最小自乗法

質問

実験結果をグラフにプロットしたところ、二次関数のような形になりました。
ここで
Y=aXX+bX=c
とおいたとき、最小自乗法から a,b,cの値を求めるにはどうしたら
いいのでしょうか?

回答No.6

noname#259682 · ベストアンサー

【基本的な考え方】
・二つの変数(x , y)の観測値がｎ組み得られてて、
その関係が y = ax^2 + bx +c と言う式で関係
付けられる時、ｎ組の値を x と y を関係付ける
式の代入して、観測値yと x と y を関係付ける式
　で得られる値との差を記すと

e1 = y1 - ( ax1^2 + bx1 + c )
e2 = y2 - ( ax2^2 + bx2 + c )
: : : : : : :
en = yn - ( axn^2 + bxn + c )

・これらの、二乗和を求めると

E = e1^2 + e2^2 + + en^2

・この E を最小にする、(a , b , c )を求めるのが
ご質問の回答ですね。

【具体的な手順】
・上に記した方向付けで、手順を示します。
(Σの範囲は、1～n です。)
E=Σ(yi-(axi^2 + bxi + c))^2

・Eを最小にする a,b,c を求めるには、上の式を
a,b,c で偏微分して極値を求めます。

∂E/∂a=-2Σxi^2(yi-(axi^2 + bxi + c))
∂E/∂b=-2Σxi(yi-(axi^2 + bxi + c))
∂E/∂a=-2Σ(yi-(axi^2 + bxi + c))

・上の三つの式から次の連立方程式が得られます。

aΣxi^4 + bΣxi^3 + cΣxi^2 = Σxi^2yi
aΣxi^3 + bΣxi^2 + cΣxi = Σxiyi
aΣxi^2 + bΣxi + cΣ = Σyi

・この解が、求める a,b,c です。
さて、参考になったでしょうか？

お礼

はい、とても参考になりました。一次式の傾きと切片を求める公式は解ったのですが、二次以上となるとさっぱりわからなかったので、今回非常に勉強になりました。ありがとうございます。

回答No.5

noname#21649

プログラムが既に出ているので．手計算の方法を
３点を使って方程式を立てて．あ、ｂ、ｃを求めます。
これを初期値として．他の点すべてで計算上のＹと実際のＹの差を求めま．全体の２情話をけいさんします。
あ、ｂ、ｃの値を適当に変化させて．全体の２情話が少なくなる点を求めます。

まあ．答え一発カシオミニを使って計算すると．キーボートの修理時間を除いて．１ヶ月くらい．連日連夜キーボードをたたきつづけると．適当な答えが求まります。「適当に変化」させる方法にはいろいろ有ります。非線型最小二乗法の本を読んでください。東京大学出版界の本（ＵＰおうようせんしょ）が結構わかりやすいです。

お礼

「あ」は「a」,「2情話」は「自乗和」のことですよね?丁寧な解答ありがとうございました。

回答No.4

iykm

念のため、７行目と８行目の間にも、
______Y1 = 0.0D0
を入れてください。

お礼

今回は質問が曖昧だったためか、少しもとめていたものとはちがっていますが、
お忙しい中解答して下さって真にありがとうございます。

回答No.3

iykm

______IMPLICIT REAL*8(A-H,O-Z)
______X4 = 0.0D0
______X3 = 0.0D0
______X2 = 0.0D0
______X1 = 0.0D0
______X2Y = 0.0D0
______XY = 0.0D0
______N = 0
__100 WRITE (*,*) 'X,Y? ( End = Ctrl+D )'
______READ (*,*,END=200) X,Y
______X4 = X4 + X**4
______X3 = X3 + X**3
______X2 = X2 + X**2
______X1 = X1 + X
______X2Y = X2Y + X*X*Y
______XY = XY + X*Y
______Y1 = Y1 + Y
______N = N + 1
______GOTO 100
__200 A = (X3*(X1*Y1-N*XY)-X2**2*Y1+X2*(X1*XY+N*X2Y)-X1**2*X2Y)
_____& /(X4*(X2*N-X1**2)-X3**2*N+X2*(2*X3*X1-X2**2))
______B = -(X4*(X1*Y1-N*XY)+X3*(N*X2Y-X2*Y1)+X2*(X2*XY-X1*X2Y))
_____& /(X4*(X2*N-X1**2)-X3**2*N+X2*(2*X3*X1-X2**2))
______C = (X4*(X2*Y1-X1*XY)-X3**2*Y1+X3*(X2*XY+X1*X2Y)-X2**2*X2Y)
_____& /(X4*(X2*N-X1**2)-X3**2*N+X2*(2*X3*X1-X2**2))
______WRITE (*,*) 'A =',A
______WRITE (*,*) 'B =',B
______WRITE (*,*) 'C =',C
______END

飲み会帰りの酔い覚まし中ですが作ってみました。Fortran です。
先頭の下線＿は空白に置き換えてください。

回答No.2

siegmund

最小自乗法の質問は過去たくさんありますので，
質問検索で検索してみてください．

http://oshiete1.goo.ne.jp/kotaeru.php3?q=92718
http://oshiete1.goo.ne.jp/kotaeru.php3?q=24627
あたりは重要です．

質問の
Y = aX^2 + bX + c　　(ですよね)
は最も基本的な形です．
大学初年級の物理実験書などには大抵載っています．
たぶん，学生実験ですよね．
実験をやった後のデータ処理で，最小自乗法は基本中の基本です．
これからいくらでも出てくる事項ですので，
この際適当な本で基本からマスターされるようにおすすめします．

Excel なら一発，というのは pm-770c さんおっしゃるとおりですが．

お礼

一日でも早く最小自乗法を使いこなせるように頑張ります。具体的な指導に
とても感謝しています。ありがとうございました。

回答No.1

pm-770c

もし、Excelがあるなら簡単に出来ます。

お礼

おっしゃる通りですね。的確な助言ありがとうございます。