最小二乗フィット、カイ二乗フィット、gnuplot

数値計算法

2009 5/27

林田 清

（大阪大学大学院理学研究科）

最尤法 (Maximum Likelihood

Method)

1 2 2

n

,

,....,

μGauss)

(

)

1

1

exp

2

2

'

'

n i i i i i i

x x

x

x

x

dx

dQ

Pdx

x

P

σ

µ

σ

σ

π

µ

µ

µ

+

=



_

₋

_



≡

_

−

_

_

_













回の（独立な）測定

で、

母集団が平均値 　　標準偏差 の正規（

分布の場

合

１回の測定で　 から

　の間の値を観測する確率は

　は不可知、推定値をとする。

尤度が最大になるはどう決められるか

最尤法２

最尤法（正規分布の場合の例） 2 1 2 1 2 1 ' ' ' 1 1 ( ') exp 2 2 , ,..., ( ') ( ') ' 1 1 exp 2 2 ( ') ' i i i n n i i n n i i x x P n x x x P P x P µ σ σ µ µ σ σ π µ µ µ σ σ π µ µ µ = = =  _{ − }  = − _{ }        =  ₋      = _{  }− _{  }     _ _

∏

∑

平均値 、標準偏差 の正規分布を仮定すると を観測する確率は 回の測定で を観測する確率（尤度）は を最大にする が最も確からしい の推定値 考え方： 最も確率の高い標本分布（測定 値の組）が実現されているはず

最尤法３



最も確からしい母集団平均(mean)の推定値は加

算平均(average)

2 1 2 1 1 ( ') ' 1 2 ' 0 ' 1 ' n i i n i i n i i P X x X x dX d x x n

µ

µ

σ

µ

µ

σ

µ

= = = −   = _{ }   −   = − _{ } =   = =

∑

∑

∑

を最大にすることは次の を最小にするのと同じ

誤差が異なるデータの場合

（重みつき平均）

2 1 1 2 ₂ 2 2 1 1 2 ' ₂ ' 1 1 ( ') exp 2 2 ' ( / ) ' ' 1 0 ' ' 2 (1/ ) 1 ' (1/ ) i i n n i i i _i i n n i i i i i _i i _{i i} i x x P x x x d d µ

σ

µ

µ

σ

σ

π

µ

σ

µ

µ

_µ

µ

σ

σ

σ

µ

σ

σ

= = = =     _{ − }   = _{ } − _{ }  _{ }    _{ }    −  − = − _{ } = =       =

∑

∏

∑

∑

∑

_∑

∑

各測定値 につく誤差が異なる の場合 の最尤推定値は より また推定値 に関する誤差は

どうやって誤差を評価するか？



例えば使用説明書に書いてある測定器の精度を使

用するのは一般には不十分

 Conservativeな測定値の範囲を示すには有効 

一般に系統誤差を評価するのは困難

 全く独立な実験を行い結果を比較する 

測定を

同じ条件で

複数回繰り替えすことができる場

合は測定値の（標本）標準偏差が（偶然）誤差の推

定値を与える

 最尤法を使い誤差を推定することもできる 

統計誤差の場合、理論的に推定できることがある

 例）放射線源を決まった時間だけ計測する際の計数ｘは ポアソン分布に従う。 この場合統計誤差は√x。

レポート課題１（締め切りは6/3）

 平均値と標準偏差を求めるプログラム  入力：データの数、データ  データは以下の１０個（例えばある月の最高気温（℃）１０日分）  35.5,25.0,34.2,34.6,22.8,27.7,30.6,26.8,23.0,39.3  出力：（標本）平均値、標準偏差  ソースプログラムと出力結果をメイルの本文にして [email protected] までメイルせよ。 実行形 式は添付しないこと。 メイルのタイトルは report1_学籍番 号とすること。  他の人のソースプログラムと結果をそのままコピーするのは ダメ。

データ入力



Cの場合

int n; double a; scanf("%d %lf",&n,&a); 

Fortranの場合

c23456 integer n real*8 a read(*,*) n,a 1. 可能であれば、データをファイルに書き込み、ファイルから読 み出すようにプログラムしてみよう。

2. fopen, fscanf, fcloseを利用するのが一つの方法。 授業では リダイレクト(<あるいは>）を利用する方法を紹介した。

3. ファイルに含まれるデータの数が任意の場合にも対応できる ようにできるか？

繰り返しループ

 Cの場合 int i,n; double a,sum; scanf("%d",&n); sum=0;

for(i=0; i<n; i++) {

scanf("%lf",&a); sum = sum+a； }  Fortranの場合 c23456 integer i,n real*8 a,sum read(*,*) n sum=0 do 10 i=1,n read(*,*) a sum = sum+a 10 continue

配列の利用

 Cの場合

int i,n;

double a,amat[100]; scanf("%d",&n);

for(i=0; i<n; i++ ) { scanf("%lf",&a); amat[i] = a； }  Fortranの場合 c23456 integer i,n real*8 a,amat(100) read(*,*) n sum=0 do 10 i=1,n read(*,*) a amat(i)=a 10 continue

#include <stdio.h> #include <math.h> int main(){

int i, n;

double a,amat[100], sum; scanf("%d",&n); for(i=0;i<n;i++) { scanf("%lf",&a); amat[i]=a; } sum = 0.0; for(i=0;i<n;i++) { sum=sum+amat[i]; printf("%d %lf¥n",i,amat[i]); データを配列に読み込み、和をとるプログラムの一例；宿題の参考

データのモデル化、あてはめ

(Fit)、回帰

 ばらつきのある測定値に適当なモ デル（直線や曲線）であてはめるこ と  モデル  直線の場合。。。線形回帰  多項式の場合  一般の関数の場合  データの誤差  各点共通の場合  各点で重みが異なる場合  モデル点のまわりのばらつき  正規分布の場合  それ以外の場合 0 5 10 15 0 2 4 6 8 10 X -1 0 1 2 3 4 5 0 2 4 6 8 10 X

最小二乗フィット

(例：直線モデル） １

0 0 0 0 0 0 ( ) , , ( ) ( ) i i i i i x y x y y x ax b a b a b y x a x b y y x σ = + = + 測定値の組( , )があり、独立変数 と従属変数 の間の関係を 　 で近似するとき　、 に関する最も確からしい推定値は どうやって決められるか？ 母集団における係数を とし、”真”の関係式を さらに測定値 は平均値 、標準偏差 の 正規分布に従うと仮定する。 0 5 10 15 0 2 4 6 8 10 X 正規分布に従う 母集団から標本 を1個採ってくる のが測定

最小二乗フィット

(例：直線モデル） ２

2 0 2 0 0 0 1 1 1 ( ) 1 1 exp 2 2 ( ) 1 1 ( , ) exp 2 2 , ( , i i i i i i i i n n n i i i i i i _{i i} i y P y y x P n y y y x P a b P a b y P a b σ σ π σ σ π = = =  _{ − }    = _{−  }  _{ }          −    _{ } = = _{  } − _{ }  _{ }       _

∑

_

∏

∏

を観測する確率（密度） は 個の観測値 の組を得る確率（密度）は 同様に任意の係数推定値 に従うときに観測値 の組を得る確率（密度）は 2 1 1 0 0 0 0 ( ) 1 1 ) exp 2 2 ( , ) ( , ) ( , ) ( , ) n n i i i i _i i y y x P a b P a b a b a b σ σ π = =    _{ } _{ − }   _{ } = _{  } − _{ }  _{ }     

∏

 _

∑

_ 観測は母集団 採取する操作。　から の最大値を与えるような が 定値。の最尤推

最尤法の考え方

最小二乗フィット

(例：直線モデル） ３

2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 2 2 0, 0 , 1 1 1 1 ( ) ( , ) i i i i i i i i n n i i i i i i i i i i i i i i i i i i i _i i _i a b x y x y a x y x x y y y x y ax b P a b b x x χ χ χ σ σ σ σ σ σ σ χ σ σ σ σ σ σ χ = = ∂ ₌ ∂ ₌ ∂ ∂   = _ − _ ∆     = _ − _ ∆     ∆ = −  −   − −  ≡ =          

∑ ∑

∑ ∑

∑ ∑

∑

∑

∑

∑

∑

a b から を最小 を最大にする にす ＝ を最小にす とし ただ る る て し 2  

∑ ∑

二乗の和を最小にするので 最小二乗フィットと呼ぶ。 χ2_{フィットともいう。} ( ) ( ) ( ) ( ) ( ) 2 2 2 1 2 2 1 , ) , , ) ( ) 1 , 1 n i i i i i i i i i i i i i i i i i i y ax b a b x y x a n x y x y b x y x x y n x b x ax a b χ σ = = − = − ∆ = − ∆ = − + ∆ − ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ 2 各点の誤差が同一のとき を最小にする( を 求めることは、各測定点 とモデル点 ( の距離のニ乗和を最小にする を求める ただ 価 し ことと等

あてはめの良さ

(Goodness of Fit)

2 2 2 1 1 2 2 2 ( ) ( , -, ( / ) n n i i i i i _i i _i y y x y ax b n m m a b n m a b ν

χ

σ

σ

χ

ν

χ

χ ν

= =  −   − −  ≡ _{   }     − ≡

∑

（直線モデルの場合

∑

） は自由度 はパラメータの数、直線の場合 で２） の に従う。　期待値は 。 これがあてはめの良さ（仮定したモデル関数の妥当性、 パラメータ が適当であること、測定誤差が正しく評価 されていること）の基準になる。 を自 カイ自乗分布 由度 で割った を_{reduced chi-square}という。 ( ) i i i i y y x t

σ

− = は中心0,標準偏差１の正規分布に従う gnuplotのfitでは自由度はdegrees of freedom (ndf) : として

χ

2

_分布

2 2 2 2 2 2 2 / 2 1 / 2 / 2 2 2 2 2 2 2 0 1 ( ) {( ) }/ 2 ( / 2 ) ( ) ( ) 2 ( ) i i x x n e E V x n n ν χ ν ν χ χ ν χ χ χ ν χ ν χ ν µ µ σ χ χ σ − − = = Γ = = −

∑

∑

n i=1 n i=1 平均値 ,標準偏差 の正規分布に従う変数 の自乗和 の従う 分布を自由度 の 分布と呼ぶ。　 一般に自由度 の 分布は f 期待値　 分散　　 平 自由度 の （カ 均値 ,標準偏差 の正規分布に従う 　も自由度 の イ二乗）分布 分布 2 2 2 ( ) 1 i x x n χ σ − −

∑

n i=1 、　 はしかし自由度 の 分布

カイ二乗分布の確率分布の積分

あてはめの良さの検定

Data Reduction and Error Analysis for the Physical Sciences, Bevington & Robinson より

• 最小二乗フィット によりモデルパ ラメータを最適 化した際のχ2値 を求める • 上記のχ2値（以 上の値）を得る 確率を表から調 べる。 • 確率があまりに も小さければ何 か間違っている。 （例えばモデル が適当でない） ｒeduced-χ2の値の表（対応するχ2の 値を超える確率Pと自由度νの関数と して表示されている）

パラメータの推定誤差

2 2 2 1 2 2 2 2 1 1 1 1 n a i i _{i i} n i b i i _{i i} a y x b y σ σ σ σ σ σ = =    _∂  = _{  =  } ∂ ∆        ∂  = _{  =  } ∂ ∆    

∑

∑

∑

∑



最適化したパラメータはあくまでもパラメータの

真の値の推定値。 必ず推定誤差がある。



直線モデルの場合、誤差伝播側より計算できる

参考

任意関数の最小二乗（カイ二乗）フィット

2 2 1 2 2 2 2 2 min 2 2 min 2 min ( ) ( ) 1 n i i i i y x y y x m n m a a a a a a a χ χ χ χ σ χ ν χ χ χ χ = + −  −  ≡ _{ }   − ∆ = + ∆ − ∆

∑

任意の関数形 をモデルに採用した場合でも を最小にするようパラメータを決定する。 パラメータの数を として は自由度 = の 分布に従うことが期待される。 パラメータの誤差の推定: を最小にするパラメータ値 に対して、 を１だけ増加させる （ ） の値、 、 を探す。 の誤差範囲（１パラメータ68%信頼水準）はa_{χ2 min} − ∆a₋からa_{χ2 min} + ∆a₊。

カイ二乗フィットのパラメータ誤差推定

（パラメータの数による信頼区間の違い）

Numerical Recipes in C, 技術評論社より転載。 上の表で自由度とは（注 目する）パラメータの数。 パラメータa₁,a₂それぞれのの 68%信頼区間はΔχ2_=1である が、(a₁,a₂)の組の68%信頼区 間はΔχ2_{=2.3の楕円で囲まれ} 参考

グラフの書き方練習



gnuplot



“端末”で gnuplotとうつと起動する



簡単関数やファイルに書き込んだデータをプロットできる



使い方は help というコマンドで参照できる



インターネットで参照できる日本語のマニュアルもあり

 http://lagendra.s.kanazawa-u.ac.jp/ogurisu/manuals/gnuplot-intro/  http://lagendra.s.kanazawa-u.ac.jp/ogurisu/manuals/gnuplot/ 

StarSuite8 のスプレッドシート(Calc)



Windows, MacのExcelに近い使用方法



誤差を考慮した複雑な解析には使いにくい？

gnuplotの練習



データファイル

（例えば）

xye.datを用意

する。

1.0 2.2 0.2 2.0 2.9 0.1 3.0 4.3 0.3 4.0 5.0 0.1 5.0 5.8 0.3 

端末でgnuplotとうつと起動する。続

いて以下の操作を試してみる。

 plot "xye.dat" using 1:2:3 with yerrorbars  set xrange[0.0:7.0]

 set yrange[0.0:7.0]  f(x)=a*x+b

 fit f(x) "xye.dat" using 1:2:3 via a,b

 （ここで表示されるフィット結果を理解せよ）  replot f(x)

 次に、各点の誤差を無視した（重みづけなしの）フィットを 比較のためやってみると

 g(x)=c*x+d

 fit g(x) "xye.dat" using 1:2 via c,d  replot f(x),g(x)