数値計算法

Transcription

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

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

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

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

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

6 にも同様の表 ( 但し reduced ch-squared ではなく ch-squared の値 ) が掲載されている Excel なら CHIDIST,CHIINV

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

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

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

10 グラフの書き方練習 gnuplot 端末 で gnuplotとうつと起動する簡単関数やファイルに書き込んだデータをプロットできる使い方は help というコマンドで参照できる インターネットで参照できる日本語のマニュアルもあり

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

12 最小二乗 ( カイ二乗 ) フィットのまとめ 最尤法が根拠 ただし 測定値 y のモデル点からのばらつきが正規分布で近似できる場合に限定 を最小にするパラメータが最良推定値 あてはめの良さ モデルの妥当性は の値が自由度 n-m に近いかどうかで評価できる パラメータの誤差 ( 信頼区間 ) は から推定できる

13 参考 ) 最尤法の直接的な利用 1 K 0 中間子の寿命の測定 K 0 中間子の生成点は生成に伴う二次荷電粒子の飛跡から 崩壊点と運動量は崩壊後のパイ中間子の飛跡と運動量の測定から決められる点線の領域内で崩壊が起こった現象だけ取り扱う観測した崩壊イベントの平均が寿命の最良推定値になるか? No Data Reducton and Error Analyss for the Physcal Scences, Bevngton & Robnson より

14 参考 ) 最尤法の直接的な利用 0 時間 t だけ生き延びるK 中間子を観測する確率 P A p( t ; ) A e t / t / ( ; ) は寿命 の粒子が ~ ここでAは 定められた領域内で崩壊が起こり検出できる効率 0 K 中間子の生成点 崩壊点の位置や運動量 寿命 に依存する p t e t t dtの間に崩壊する確率 Aはtやと独立ではないことに注意 0 生成点と運動量が決まっているK 中間子に対して 点線領域に入るまでの 距離を d, 出る ( 崩壊が起こらなかったとして ) までの距離をdとし 対応する 時間をt, t とする Aは次のように規格化する t t a b a a b tb t / t Pdt A e dt N L( ) P A e 1 1 a N 1 N個のイベントについて尤度は t / これを最大にするようなが求めたい答え b

15 参考 ) 最尤法の直接的な利用 3 L( ) P A e 1 1 のかわりに t M ( ) ln L( ) ln A を最大にすることを考える 例 1) t 0, t のとき ( 粒子の寿命に対して測定領域が十分大きい場合 ) 1 A 1/ でM ( ) t N ln dm ( ) 1 N t 0 / より t N d 例 )( t 0), t ( ) が共通の値である場合 A 1/ N a a t 0 b e t / b N dt b t / 1 1 e t / b t t M ( ) ln L ln A N ln N ln 1 e dm ( ) 0をみたすが寿命の最尤推定値 d t / b

16 最尤法の直接利用と最小二乗法 最小二乗法を使えないとき= 分布が正規分布でないとき ビンまとめし ヒストグラムをつくると 1ビンあたりに含まれるデータ数が十分大きい場合 正規分布で近似できる この場合最小二乗法が使えるようになる ただし もともとのデータ数が小さい場合は適用付加 最尤法の直接利用複雑なモンテカルロ計算が必要になるような場合 ( 例 :K 中間子の寿命測定 ) も最尤法の直接利用が効果的 M=1/ より最尤法で決めたパラメータ誤差を推定できる しかし 最尤法の直接利用ではあてはめの良さを評価する適当な指標 ( 最小二乗法の χ のような ) がない

17 確率分布 検定 区間推定 いろいろな確率分布 二項分布 ポアッソン分布 正規 ( ガウス分布 ) t 分布 χ 乗分布 統計的検定 仮説の当否を統計的に検証する 区間推定 真の値の範囲を統計的に推定する 相関係数 個のパラメータ間の関連を調べる

18 二項分布 ポアッソン分布 二項分布 n! x PB ( x; n, p) p (1 p) ( n x)! x! x pn np(1 p) (1 p) ポアッソン分布二項分布でp 1の極限 Px ( ; ) x x e x! nx Posson Dstrbuton x 参考 ) 視聴率の誤差について

19 ポアッソン分布の導出その 1 二項分布 n! 1 n! PB ( x; n, p) p (1 p) p (1 p) (1 p) ( n x)! x! x! ( n x)! pn n 1/ p 1 (1 p) lm (1 p) e p0 p0 e lm lm p0 np(1 p) (1 p) x PB ( x; n, p) Pp ( x; ) e x! x n x x x n においてを一定に保ったまま p 1の極限を考える n! x n (for x n) ( n x)! x (1 p) 1 px

20 ポアッソン分布 Posson Dstrbuton ポアッソン分布の例 0 放射線源の1 秒あたりの崩壊数放射線源の測定で1 時間当たりの検出カウント数 1000 人の集団の中で今日が誕生日の人の数 ポアッソン分布の統計誤差 平均値の平方根 ( 複数回の測定ができないとき )1 回の測定値の平方根で置き換えるときもある ポアッソン分布と正規分布 平均値 が大きいとき ( 例えば 0 以上 ) ではポアッソン分布は平均値 分散 の正規分布で近似できる x

21 正規分布 1 Px ( ;, ) exp x Bevngton &Robnson より

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

23 統計的検定 (statstcal 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を伴って 大きい : 仮説は否定できない 危険率 ( 有意水準 )=sgnfcance level

24 いろいろな検定 母平均の検定 : 正規分布 母集団の分散 が既知でない場合 ->t 分布 母平均の差の検定 ->t 分布 母分散の検定 : 分布 母分散の比の検定 :F 分布 相関の有無の検定 : 相関係数の表

25 区間推定 f(t) 例 )n 回の測定の平均値がxと求まったとき母平均の存在する範囲はどのように推定できるか? -t ( /) N-1 1- / +t ( /) N-1 t 母集団の分布は正規分布 (, ) と仮定すると 標本平均は 正規分布 (, /n) に従う ( x ) / s / nは自由度 n 1の t分布に従う 確率 1-となる区間は -t ( /) ( t ( /) N-1 x ) / s / n N-1 変形して x-t N-1( /) s / n x t N-1( /) s / n が信頼係数 100 (1- ) % での母平均 の信頼区間 nが大きいときにはt 分布のかわりに正規分布を使い x-z( /) s / n x z( /) s / n で近似するときもある 信頼区間 =confdence nterval 信頼係数 =confdence level

26 信頼区間の推定 正規分布の場合 -<x-< にくる確率 68.3% -<x-< にくる確率 95.5% -3<x-<3 にくる確率 99.7% -1.96<x-<1.96 にくる確率 95% -.58<x-<.58 にくる確率 99%

27 相関係数 二つの測定量 x,y の間に ( 線形 ) 相関があるかどうか 1に近ければ正の相関 -1に近ければ負の相関 ゼロなら相関なし r N x y x y N x x N y y 1/ 1/ r=0.89 r=-0.05 r=-0.95

28 相関係数の検定 Data Reducton and Error Analyss for the Physcal Scences, Bevngton & Robnson より