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

Transcription

1 数値計算法 009 5/7 林田清 ( 大阪大学大学院理学研究科 )

2 最尤法 (Maxmum Lkelhood Method) 回の ( 独立な ) 測定 xで, x,..., x 1 母集団が平均値 μgauss) 標準偏差 の正規 ( 分布の場合 1 回の測定で xから( xの間の値を観測する確率は + dx) dq = Pdx 1 1 x µ P exp π µ は不可知 推定値をとする µ ' 尤度が最大になるはどう決められるか µ '

3 最尤法 平均値 µ ' 標準偏差 ' = の正規分布を仮定すると xを観測する確率は 1 1 x µ ' P ( µ ') = exp π 回の測定でx, x,..., xを観測する確率 ( 尤度 ) は = 1 1 P( µ ') = P( µ ') 1 1 x µ ' = exp π = 1 P( µ ') を最大にするµ ' が最も確からしいµ の推定値 最尤法 ( 正規分布の場合の例 ) 考え方 : 最も確率の高い標本分布 ( 測定値の組 ) が実現されているはず

4 最尤法 3 最も確からしい母集団平均 (mea) の推定値は加算平均 (average) ( ') ' 1 ' 0 ' 1 ' P X x X x dx d x x µ µ µ µ µ = = = = = = = = を最大にすることは次のを最小にするのと同じ

5 誤差が異なるデータの場合 ( 重みつき平均 ) 各測定値 xにつく誤差が異なる の場合 1 1 x µ ' P( µ ') = exp = 1 π µ ' の最尤推定値は = 1 x µ ' x µ ' x = 0 µ ' = より = = 1 = 1 d 1 ( / ) dµ ' (1/ ) 1 また推定値 µ ' に関する誤差は (1/ ) µ ' =

6 どうやって誤差を評価するか? 例えば使用説明書に書いてある測定器の精度を使用するのは一般には不十分 Coservatve な測定値の範囲を示すには有効 一般に系統誤差を評価するのは困難 全く独立な実験を行い結果を比較する 測定を同じ条件で複数回繰り替えすことができる場合は測定値の ( 標本 ) 標準偏差が ( 偶然 ) 誤差の推定値を与える 最尤法を使い誤差を推定することもできる 統計誤差の場合 理論的に推定できることがある 例 ) 放射線源を決まった時間だけ計測する際の計数 x はポアソン分布に従う この場合統計誤差は x

7 レポート課題 1( 締め切りは 6/3) 平均値と標準偏差を求めるプログラム 入力 : データの数 データ データは以下の10 個 ( 例えばある月の最高気温 ( )10 日分 ) 35.5,5.0,34.,34.6,.8,7.7,30.6,6.8,3.0,39.3 出力 :( 標本 ) 平均値 標準偏差 ソースプログラムと出力結果をメイルの本文にして khclass@ess.sc.osaka-u.ac.jp までメイルせよ 実行形式は添付しないこと メイルのタイトルは report1_ 学籍番号とすること 他の人のソースプログラムと結果をそのままコピーするのはダメ

8 データ入力 C の場合 t ; double a; scaf("%d %lf",&,&a); Fortra の場合 c3456 teger real*8 a read(*,*),a 1. 可能であれば データをファイルに書き込み ファイルから読み出すようにプログラムしてみよう. fope, fscaf, fclose を利用するのが一つの方法 授業ではリダイレクト (< あるいは >) を利用する方法を紹介した 3. ファイルに含まれるデータの数が任意の場合にも対応できるようにできるか?

9 繰り返しループ C の場合 t,; double a,sum; scaf("%d",&); sum=0; for(=0; <; ++) { scaf("%lf",&a); sum = sum+a; } Fortra の場合 c3456 teger, real*8 a,sum read(*,*) sum=0 do 10 =1, read(*,*) a sum = sum+a 10 cotue

10 配列の利用 C の場合 t,; double a,amat[100]; scaf("%d",&); for(=0; <; ++ ) { scaf("%lf",&a); amat[] = a; } Fortra の場合 c3456 teger, real*8 a,amat(100) read(*,*) sum=0 do 10 =1, read(*,*) a amat()=a 10 cotue

11 データを配列に読み込み 和をとるプログラムの一例 ; 宿題の参考 #clude <stdo.h> #clude <math.h> t ma(){ t, ; double a,amat[100], sum; scaf("%d",&); for(=0;<;++) { scaf("%lf",&a); amat[]=a; } sum = 0.0; for(=0;<;++) { sum=sum+amat[];

12 データのモデル化 あてはめ (Ft) 回帰 ばらつきのある測定値に適当なモデル ( 直線や曲線 ) であてはめること モデル 直線の場合 線形回帰 多項式の場合 一般の関数の場合 データの誤差 各点共通の場合 各点で重みが異なる場合 モデル点のまわりのばらつき 正規分布の場合 それ以外の場合 X X

13 最小二乗フィット ( 例 : 直線モデル ) 1 測定値の組 ( x, y ) があり 独立変数 xと従属変数 yの間の関係を y( x) = ax + b で近似するとき ab, に関する最も確からしい推定値はどうやって決められるか? 母集団における係数をa, bとし 真 の関係式を y( x) = ax+ b さらに測定値 yは平均値 y ( x ) 標準偏差 の 0 正規分布に従うと仮定する 正規分布に従う母集団から標本を 1 個採ってくるのが測定 X

14 最小二乗フィット ( 例 : 直線モデル ) yを観測する確率 ( 密度 ) Pは 1 1 y y0( x) P = exp π 個の観測値 yの組を得る確率 ( 密度 ) は 1 1 y y0( x) Pa ( 0, b0) = P = exp = 1 = 1 π = 1 同様に任意の係数推定値 ab, に従うときに観測値 yの組を得る確率 ( 密度 ) は 1 1 y yx ( ) Pab (, ) = exp = 1 π = 1 観測は母集団 Pa (, 採取する操作 b) から 0 0 Pab (, ) の最大値を与えるような( ab, ) が( a, b定値 ) の最尤推 0 0 最尤法の考え方

15 最小二乗フィット ( 例 : 直線モデル ) 3 χ y( x) y ax b y = = 1 = 1 Pab (, ) を最大にする=χ を最小にする χ = 0, χ = 0 a b からχ を最小にするab, として 1 1 xy x y a = 1 x y x x y b = 1 x x ただし = 二乗の和を最小にするので最小二乗フィットと呼ぶ χ フィットともいう 各点の誤差が同一のとき 1 χ = y ax b を最小にする ( a, b) を = 1 ( ) ( x, y ) 求めることは 各測定点 とモデル点 ( x, ax + b) の距離のニ乗和を最小にする( a, b) を求めることと等価 1 a= ( xy x y ) 1 b= x y x xy ただし ( ) = x ( x)

16 あてはめの良さ (Goodess of Ft) t y yx ( ) = は中心 0, 標準偏差 1の正規分布に従う ( ) y y x y ax b χ ( 直線モデルの場合 ) = 1 = 1 は自由度 mm ( はパラメータの数 直線の場合 ab, で) のカイ自乗分布に従う 期待値は - m これがあてはめの良さ ( 仮定したモデル関数の妥当性 パラメータab, が適当であること 測定誤差が正しく評価されていること ) の基準になる χ を自由度 νで割ったχ ( χ / ν) をreduced ch-squareという ν guplot の ft では自由度は degrees of freedom (df) : として reduced χ は varace of resduals (reduced chsquare) = WSSR/df : として表示されている

17 χ 分布 自由度の ( カ χ 平均値 0, 標準偏差 1の正規分布に従う変数 xの自乗和 χ = x =1 分布を自由度の分布と呼ぶ 一般に自由度の分布は f χ = χ e Γ ν ν / 1 χ / ν / ν ( ) {( ) }/ ( /) 期待値 E イ二乗 ) 分布 χ ν χ χ = ν 分散 V χ = ν ( ) ( ) ( x µ ) 平均値 µ, 標準偏差 の正規分布に従う も自由度 のχ 分布 =1 ( x x) =1 はしかし自由度 1の分布 χ の従う

18 カイ二乗分布の確率分布の積分あてはめの良さの検定 reduced-χ の値の表 ( 対応する χ の値を超える確率 P と自由度 ν の関数として表示されている ) 最小二乗フィットによりモデルパラメータを最適化した際の χ 値を求める 上記の χ 値 ( 以上の値 ) を得る確率を表から調べる Data Reducto ad Error Aalyss for the Physcal Sceces, Bevgto & Robso より 確率があまりにも小さければ何か間違っている ( 例えばモデル

19 パラメータの推定誤差 参考 最適化したパラメータはあくまでもパラメータの真の値の推定値 必ず推定誤差がある 直線モデルの場合 誤差伝播側より計算できる a 1 1 = = a = 1 y = = b 1 x b = 1 y

20 任意関数の最小二乗 ( カイ二乗 ) フィット 任意の関数形 yx ( ) をモデルに採用した場合でも y yx ( ) χ = 1 を最小にするようパラメータを決定する パラメータの数をとしては自由度 = の分布に従うことが期待される m χ ν m χ パラメータの誤差の推定 : χ を最小にするパラメータ値 a に対して χ を1だけ増加させる χ m ( ) の値 を探す χ = 1 a aχm + a+ aχm a aの誤差範囲 (1パラメータ68% 信頼水準 ) はa m a から a m + a χ χ +

21 カイ二乗フィットのパラメータ誤差推定 ( パラメータの数による信頼区間の違い ) 参考 パラメータ a 1,a それぞれのの 68% 信頼区間は Δχ =1 であるが (a 1,a ) の組の 68% 信頼区間は Δχ =.3 の楕円で囲まれ Numercal Recpes C, 技術評論社より転載 上の表で自由度とは ( 注目する ) パラメータの数

22 グラフの書き方練習 guplot 端末 で guplotとうつと起動する簡単関数やファイルに書き込んだデータをプロットできる使い方は help というコマンドで参照できる インターネットで参照できる日本語のマニュアルもあり StarSute8 のスプレッドシート (Calc) Wdows, MacのExcelに近い使用方法 誤差を考慮した複雑な解析には使いにくい?

23 guplot の練習 データファイル ( 例えば ) xye.dat を用意する 端末で guplot とうつと起動する 続いて以下の操作を試してみる plot "xye.dat" usg 1::3 wth yerrorbars set xrage[0.0:7.0] set yrage[0.0:7.0] f(x)=a*x+b ft f(x) "xye.dat" usg 1::3 va a,b ( ここで表示されるフィット結果を理解せよ ) replot f(x) 次に 各点の誤差を無視した ( 重みづけなしの ) フィットを比較のためやってみると g(x)=c*x+d ft g(x) "xye.dat" usg 1: va c,d replot f(x),g(x)