• 検索結果がありません。

数値計算法

N/A
N/A
Protected

Academic year: 2021

シェア "数値計算法"

Copied!
28
0
0

読み込み中.... (全文を見る)

全文

(1)

数値計算法

2011/5/25

林田 清

(大阪大学大学院理学研究科)

(2)

レポート課題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別データで計算しよう

(3)

あてはめの良さ

(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) : として

(4)

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 、  はしかし自由度 の 分布

(5)

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

あてはめの良さの検定

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

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

(6)

• http://cluster.f7.ems.okayama-u.ac.jp/~yan/jscscd/table/chi.htmlに

も同様の表(但しreduced chi-squaredではなくchi-squaredの値)が掲 載されている。

(7)

パラメータの推定誤差

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                                  

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

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

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

参考

(8)

任意関数の最小二乗(カイ二乗)フィット

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

(9)

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

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

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

(10)

グラフの書き方練習

gnuplot

 “端末”で gnuplotとうつと起動する  簡単関数やファイルに書き込んだデータをプロットできる  使い方は help というコマンドで参照できる  インターネットで参照できる日本語のマニュアルもあり  http://lagendra.s.kanazawa-u.ac.jp/ogurisu/manuals/gnuplot-intro/  http://lagendra.s.kanazawa-u.ac.jp/ogurisu/manuals/gnuplot/

(11)

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)

(12)

最小二乗(カイ二乗)フィットのまとめ

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

(13)

最尤法の直接的な利用1

 K0中間子の寿命の測定  K0中間子の生成点は生成 に伴う二次荷電粒子の飛跡 から、崩壊点と運動量は崩 壊後のパイ中間子の飛跡と 運動量の測定から決められ る  点線の領域内で崩壊が起 こった現象だけ取り扱う  観測した崩壊イベントの平 均が寿命の最良推定値にな るか? No

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

Robinson より 参考)

(14)

最尤法の直接的な利用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            

 

離を 出る(崩壊が起こらなかったとして)までの距離を とし、対応する 時間を とする。  は次のように規格化する。 個のイベントについて尤度は これを最大にするような が求めたい答え 参考)

(15)

最尤法の直接的な利用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                            をみたす が寿命の最尤推定値 参考)

(16)

最尤法の直接利用と最小二乗法

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

(17)

確率分布、検定、区間推定

いろいろな確率分布

 二項分布、ポアッソン分布、正規(ガウス分布)  t分布、χ2乗分布 

統計的検定

 仮説の当否を統計的に検証する 

区間推定

 真の値の範囲を統計的に推定する 

相関係数

 2個のパラメータ間の関連を調べる

(18)

二項分布、ポアッソン分布

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

(19)

ポアッソン分布の導出その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      

(20)

ポアッソン分布

 ポアッソン分布の例  放射線源の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 

(21)

正規分布

Bevington &Robinsonより  2 2 1 ( ; , ) exp 2 2 x P x              

(22)

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 、  はしかし自由度 の 分布 参考)

(23)

統計的検定

(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

(24)

いろいろな検定

母平均の検定:正規分布

 母集団の分散2が既知でない場合->t分布  母平均の差の検定->t分布 

母分散の検定:

2

分布

 母分散の比の検定:F分布 

相関の有無の検定:相関係数の表

(25)

区間推定

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)

(26)

信頼区間の推定

正規分布の場合

 -<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%

(27)

相関係数

二つの測定量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

(28)

相関係数の検定

Data Reduction and Error Analysis for the

Physical Sciences, Bevington & Robinson より

参照

関連したドキュメント

This year, the world mathematical community recalls the memory of Abraham Robinson (1918–1984), an outstanding scientist whose contributions to delta-wing theory and model theory

(at least a proof can be reconstructed from it, after correcting a number of misprinted formulae). Other authors have subsequently given different proofs, see for instance [Kn1,

⑥ニューマチックケーソン 職種 設計計画 設計計算 設計図 数量計算 照査 報告書作成 合計.. 設計計画 設計計算 設計図 数量計算

Indeed, when GCD(p, n) = 2, the Gelfand model M splits also in a differ- ent way as the direct sum of two distinguished modules: the symmetric submodule M Sym , which is spanned by

Example Shapes (Young diagrams, left), shifted shapes (shifted Young diagrams, middle) and swivels (right) are

 当図書室は、専門図書館として数学、応用数学、計算機科学、理論物理学の分野の文

(注)本報告書に掲載している数値は端数を四捨五入しているため、表中の数値の合計が表に示されている合計

Ⅰ.連結業績