数値計算法
2011/5/25
林田 清
(大阪大学大学院理学研究科)
レポート課題1(締め切りは5/25)
平均値と標準偏差を求めるプログラム 入力:データの数、データ データは以下の10個(例えばある月の最高気温(℃)10日分) 34.3,25.0,32.2,34.6,22.9,27.7,30.6,25.8,23.0,31.3 出力:(標本)平均値、標準偏差 ソースプログラムと出力結果をメイルの本文にして khclass@ess.sci.osaka-u.ac.jp までメイルせよ。 実行形 式は添付しないこと。 メイルのタイトルは report1_学籍番 号とすること。 他の人のソースプログラムと結果をそのままコピーするのは ダメ。 計算結果が正しいことを別プログラムor別データで計算しようあてはめの良さ
(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
(直線モデルの場合
) は自由度 はパラメータの数、直線の場合 で2) の に従う。 期待値は 。 これがあてはめの良さ(仮定したモデル関数の妥当性、 パラメータ が適当であること、測定誤差が正しく評価 されていること)の基準になる。 を自 カイ自乗分布 由度 で割った をreduced chi-squareという。 ( ) i i i i y y x t は中心0,標準偏差1の正規分布に従う 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値(以 上の値)を得る 確率を表から調 べる。 • 確率があまりに も小さければ何 か間違っている。 (例えばモデル が適当でない) reduced-χ2の値の表(対応するχ2の 値を超える確率Pと自由度の関数と して表示されている)
• http://cluster.f7.ems.okayama-u.ac.jp/~yan/jscscd/table/chi.htmlに
も同様の表(但しreduced chi-squaredではなくchi-squaredの値)が掲 載されている。
パラメータの推定誤差
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
任意の関数形 をモデルに採用した場合でも を最小にするようパラメータを決定する。 パラメータの数を として は自由度 = の 分布に従うことが期待される。 パラメータの誤差の推定: を最小にするパラメータ値 に対して、 を1だけ増加させる ( ) の値、 、 を探す。 の誤差範囲(1パラメータ68%信頼水準)はa2 min aからa2 min a。カイ二乗フィットのパラメータ誤差推定
(パラメータの数による信頼区間の違い)
Numerical Recipes in C, 技術評論社より転載。 上の表で自由度とは(注 目する)パラメータの数。 パラメータa1,a2それぞれのの 68%信頼区間はΔχ2=1である が、(a1,a2)の組の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/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)
最小二乗(カイ二乗)フィットのまとめ
最尤法が根拠。 ただし、測定値yのモデル点からのば らつきが正規分布で近似できる場合に限定。 2を最小にするパラメータが最良推定値。 あてはめの良さ、モデルの妥当性は2の値が自由度 n-mに近いかどうかで評価できる。 パラメータの誤差(信頼区間)は 2から推定できる。最尤法の直接的な利用1
K0中間子の寿命の測定 K0中間子の生成点は生成 に伴う二次荷電粒子の飛跡 から、崩壊点と運動量は崩 壊後のパイ中間子の飛跡と 運動量の測定から決められ る 点線の領域内で崩壊が起 こった現象だけ取り扱う 観測した崩壊イベントの平 均が寿命の最良推定値にな るか? NoData Reduction and Error Analysis for the Physical Sciences, Bevington &
Robinson より 参考)
最尤法の直接的な利用2
/ / ( ; ) ( ; ) ~ i i i t i i i i i t i i i i i t P A p t A e A p t e t t dt A t 0 0 0 時間 だけ生き延びるK 中間子を観測する確率 ここで は、定められた領域内で崩壊が起こり検出できる効率、 K 中間子の生成点、崩壊点の位置や運動量、寿命 に依存する。 は寿命 の粒子が の間に崩壊する確率。 は や と独立ではないことに注意。 生成点と運動量が決まっているK 中間子に対して、点線領域に入るまでの 距 / / 1 1 , , 1 ( ) b b i a a i a b a b i t t t i i i i t t N N t i i i i d d t t A Pdt A e dt N L P A e
離を 出る(崩壊が起こらなかったとして)までの距離を とし、対応する 時間を とする。 は次のように規格化する。 個のイベントについて尤度は これを最大にするような が求めたい答え 参考)最尤法の直接的な利用3
/ 1 1 2 / 0 ( ) ( ) ln ( ) ln 0, 1 1/ ( ) ln ( ) 1 0 / ( ) 1/ i b i N N t i i i i i i a b i i i i a b t t i L P A e t M L A t t A M t N dM N t t N d t t A e のかわりに を最大にすることを考える 例1) のとき(粒子の寿命に対して測定領域が十分大きい場合) で より 例2)( 0), が共通の値である場合
/ / 1 1 ( ) ln ln ln ln 1 ( ) 0 b b i t t i i i dt e t t M L A N N e dM d をみたす が寿命の最尤推定値 参考)最尤法の直接利用と最小二乗法
最小二乗法を使えないとき=分布が正規分布でないとき ビンまとめし、ヒストグラムをつくると、1ビンあたりに含ま れるデータ数が十分大きい場合、正規分布で近似できる。 この場合最小二乗法が使えるようになる。 ただし、もともとのデータ数が小さい場合は適用付加。。。 最尤法の直接利用 複雑なモンテカルロ計算が必要になるような場合(例:K中 間子の寿命測定)も最尤法の直接利用が効果的 M=1/2より最尤法で決めたパラメータ誤差を推定できる しかし、最尤法の直接利用ではあてはめの良さを評価す る適当な指標(最小二乗法のχ2のような)がない。確率分布、検定、区間推定
いろいろな確率分布
二項分布、ポアッソン分布、正規(ガウス分布) t分布、χ2乗分布 統計的検定
仮説の当否を統計的に検証する 区間推定
真の値の範囲を統計的に推定する 相関係数
2個のパラメータ間の関連を調べる二項分布、ポアッソン分布
2 ! ( ; , ) (1 ) ( )! ! (1 ) (1 ) x n x B n P x n p p p n x x x pn np p p 二項分布 2 1 ( ; ) ! x p e P x x x ポアッソン分布 二項分布で の極限 0 0.1 0.2 0.3 0.4 0 5 10 15 Poisson Distribution 1 2 3 4 5 10 x 参考)視聴率の誤差について http://www.videor.co.jp/rating/wh/07.htmポアッソン分布の導出その1
2 1/ 0 0 0 ! 1 ! ( ; , ) (1 ) (1 ) (1 ) ( )! ! ! ( )! (1 ) (1 ) ! (for ) ( )! (1 ) 1 1 (1 ) (1 ) ( ;lim
lim
lim
x n x x x n B x x n p p p B p n n P x n p p p p p p n x x x n x pn np p p p n n x n n x p px p p e e P x 二項分布 において を一定に保ったまま、 1の極限を考える , ) ( ; ) ! x p n p P x e x ポアッソン分布
ポアッソン分布の例 放射線源の1秒あたりの崩壊数 放射線源の測定で1時間当たりの検出カウント数 1000人の集団の中で今日が誕生日の人の数 ポアッソン分布の統計誤差 平均値の平方根 (複数回の測定ができないとき)1回の測定値の平方根で置き換え るときもある ポアッソン分布と正規分布 平均値が大きいとき(例えば20以上)ではポアッソン分布は平 均値、分散2 の正規分布で近似できる。 0 0.1 0.2 0.3 0.4 0 5 10 15 Poisson Distribution 1 2 3 4 5 10 x 正規分布
Bevington &Robinsonより 2 2 1 ( ; , ) exp 2 2 x P x
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 、 はしかし自由度 の 分布 参考)統計的検定
(statistical test)
例)xの10回の測定平均値が0.45、標準偏差が0.05 仮説H:(例)母集団での平均値は0.5である 本当は対立仮説H':”母集団での平均値は0.5でない”を示したい ので、Hを帰無仮説という。 H':”母集団での平均値は0.5より小さい(大きい)”の場合も有り 得る。 両側検定、片側検定。 平均値0.5標準偏差0.05の母集団から10個の標本をサン プルした場合に平均値が0.45以下になる(あるいは0.45 以下、0.55以上になる)確率Pは? Pが定められた危険率(有意水準)aより 小さい:仮説は誤り。 正しい可能性を棄てる危険性aを伴って。 大きい:仮説は否定できない。 危険率(有意水準)=significance levelいろいろな検定
母平均の検定:正規分布
母集団の分散2が既知でない場合->t分布 母平均の差の検定->t分布 母分散の検定:
2分布
母分散の比の検定:F分布 相関の有無の検定:相関係数の表
区間推定
2 2 2 2 ) / / 1 ) / / / / x x s n n t x s n x s n x s n 2 2 N-1 N-1 N-1 N-1 例)n回の測定の平均値が と求まったとき 母平均の存在する範囲はどのように推定できるか? 母集団の分布は正規分布( , )と仮定すると、標本平均は 正規分布( , /n)に従う。 ( は自由度 の 分布に従う。 確率1- となる区間は -t ( /2) ( t ( /2) 変形して -t ( /2) t ( /2) が 2 2 100 (1 - ) / / x s n x s n %での母平均 の nが大きいときにはt分布のかわりに正規分布を使い -z( /2) z( /2) で近似 信頼係数 信頼 すると 区間 きもある f(t) t 1- /2 N-1 +t ( /2) N-1 -t ( /2)信頼区間の推定
正規分布の場合
-<x-<にくる確率68.3% -2<x-<2にくる確率95.5% -3<x-<3にくる確率99.7% -1.96<x-<1.96にくる確率 95% -2.58<x-<2.58にくる確率 99%相関係数
二つの測定量x,yの間に(線形)相関があるかどうか
1に近ければ正の相関、-1に近ければ負の相関、ゼロ なら相関なし
2 1/ 2
2 1/ 2 2 2 i i i i i i i i N x y x y r N x x N y y
r=0.89 r=-0.05 r=-0.95相関係数の検定
Data Reduction and Error Analysis for the
Physical Sciences, Bevington & Robinson より