スライド 1

おしらせ
• 第二回のレポートを、学科事務室にて返却しています
– 各人、都合の良い時に、早めに取りに行ってください
• 課題レポート提出に関して
– 期限：８月４日（金）
– 提出先：物理学科事務室
数値計算法概論
第十一回
データ処理（最小２乗法）
若狭智嗣
粒子物理学講座
復習ー行列の固有値問題ー
•
固有値問題
固有値λと
固有ベクトルx
を求める
•
固有値方程式
λに関する
n次方程式の解
（一般に解くのが困難）
•
物理学でよく出てくる固有値方程式（２重振り子等）
簡単化
対角要素とそれに
隣り合う要素のみが
非０で対称
復習ー３重対角行列の固有値問題の数値解法ー
•
•
•
３重対角実対称行列の解法
– スツルムの定理（定理の説明は省略）
推奨される態度
– ３重実対称行列のような特別な（疎な）行列に対してなら、
固有値を求めるプログラムがすでにあるはずである
• （だから、スツルムの定理を知らなくても解けるはず）
多くの有用な計算プログラムが掲載された本
– 「Numerical Recipes in Fortran 77, Second Edition」
W.H.Press, S.A.Teukolsky, W.T.Vetterling, and B.P. Flannery
•
– 3重対角実対称行列の固有値を求めるプログラムが掲載
既存のプログラム（ライブラリ）を用いる際の注意
– プログラムの動作・制限を正しく理解する
• 自分の必要としているものか？
– 入力・出力を正しく理解する
• 正しく使用する
プログラム(tqli)の説明
説明を読み、正しく使おう
(/home/teacher/z6wt01in/SAMPLE/tqli.f)
プログラムの
機能
上書きされる
上書き
される
Numerical
Recipes
に載っている
プログラム例
行列の次元
（今の場合２）
プログラム例（続き）
対角要素を代入
非対角要素を代入
単位行列を生成
行に関する
ループ
列に関する
ループ
対角要素なら
非対角要素なら
配列の初期化や演算では、ループ文を積極的に用いること
プログラム例（続き）
サブルーチンを呼んで、
固有値と固有ベクトルを求める
固有値
に関する
ループ
書式（フォーマット）を使って、
表示を分かり易くしている
実際の出力と比較してみること
FORTRANのフォーマット文（書式）について
• I2
–
–
• A
–
–
自然数（integer）表示
全体で２桁（I4と書けば全体で４桁）
文字列（ASCIIキャラクタ）
全体の長さは、実際の文字列が収まるように調整
• ‘,’なら１文字
• ‘-th value:’なら１０文字
• F8.4
– 小数点（floating point）表示
– 全体で８桁、小数点以下４桁（F10.5なら全体で１０桁、小数点以下５桁）
• E12.5
– 指数（exponential）表示
– 全体で１２桁、小数点以下５桁
• πなら、0.31416e+01と表示
実行結果の例
  2 2
  2 2


x2   2  1 x1
x2 


2  1 x1
書式を指定することにより、小数点の位置を揃えて見やすくしている
統計誤差と系統誤差
•
•
•
誤差
– 真値と測定量とのズレ
– ものの長さ（真値は１つ）とノギスで測った値（測定数だけ存在）の差
統計誤差（statistical errorまたはrandom error）
– 測定値は真値の周りに分布
• ズレは正負両方ある
• ズレの平均値は（理想的には）ゼロ（分散はある）
– 統計誤差の例
• 放射線の崩壊頻度
• 「正しい」ものさしを使った複数回の測定
系統誤差（systematic error）
– 測定値が真値からずれて分布
• ズレの平均値がゼロではない（正負に偏りがある）
– 系統誤差の例
• 「正しくない」ものさしを使った複数回の測定
– ものさしの目盛りが正しくない→系統的に測定値にズレが生じる
統計誤差と系統誤差
•
•
統計誤差
– 誤差の従う分布関数は一般に分かっている（ポアソン分布・ガウス分布）
– 測定データの組み合わせや物理量の推定は、通常の統計学の手法
• 誤差伝播（AとB両方に誤差があるとき、A/Bの誤差は？）
系統誤差
– 一般には誤差の従う分布関数は分からない
– 便宜的にはガウス型（正規分布）をすると仮定する
– 複数の系統誤差が関与する場合は、
と推測することが多い
各系統誤差の要因の間に相関が無いとした事に相当
分布関数（確率密度関数）
•
•
ポアソン(Poisson）分布
– ヒストグラムに区分されたデータのチャンネル毎のカウント数など
– 分布関数（平均値＝分散＝λ、測定値＝ｘ）
– カウント数が n の場合、その統計誤差は n
ガウス(Gauss)分布
– 分布関数（平均値＝μ、分散＝σ2、測定値＝ｘ）
– ガウス分布＝平均値λが大きいポアソン分布
• λ＝１０でほぼ一致
• λ＝３０でよく一致
正規分布と誤差
•
誤差を伴う測定値の表記法
A
•
– A：平均値μの推定値
– σ：標準偏差の推定値
– １σ誤差表記という
Gauss分布の場合
– 確率変数ｘが μ±σ に含まれる確率
６８．３％
– 確率変数ｘが μ±2σ に含まれる確率
９５．３％
– 確率変数ｘが μ±3σ に含まれる確率
９９．７％
（１０００回測定すると３回、±３σ以内に入らない）
演習（標準偏差と半値幅）ー筆算・電卓ー
•
•
•
ガウス分布（正規分布）を考える
ガウス分布は、平均値μと分散σ2を指定すれば決まるが、
分散を指定する代わりに半値幅(FWHM: Full Width at
Half Maximum)を指定する場合がある
半値幅とは、ガウス分布のピークの値（x = μでの値）の半
分（半値）になる時の x の値を x1 、x2 （x2 > x1）とするとき
FWHM = x2 – x1
で与えられる。
•
FWHMと標準偏差σの間の関係を求めよ
（FWHMはσの何倍か？）
解答例
•
平均０、高さ１のガウス分布を考える
 x2 

G( x )  e xp 
2 
 2 
•
半値になる時のｘを以下の様におく
x  a
•
半値（高さの半分）は１/２なので
 (a )2  1

exp 
2 
 2  2
a   2 ln2
FWHM  2 2 ln 2   2.35
最尤法
•
•
実験により得られた、ある分布に従ってばらついている n 個のデータ
x1 , x2 ,, xn
から、物理量 a を推定する事を考える
– ばらつきはランダムである（統計的である）
– 分布関数 p(xi;a) は物理量aを与えれば計算できる
（ポアソンやガウス）
１回の測定で測定値が xi となる確率が p(xi;a)
– 独立な n 回の測定で上記のデータセットを得る確率は
n
L(a )   p( xi ; a )
Likelihood（確からしさ）関数
i 1
•
最尤法
– a の推定値として最も確からしい値（最尤値）は、Likelihood関数
L(a)が最大になるときの値である。
– a がどのような値のときに、手元のデータセットを得る確率が最大
となるか？
最尤法の例
最尤法の例（続き）
•
Likelihood 関数
n
多くの場合、
対数をとると
見通しが良い
L( )    e xp( t i )
i 1
n
n
i 1
i 1
ln L( )   ln  exp(  t i )  n ln     t i
• L(a) が極大となるλを求めるために、λで微分
dL( ) n n
   ti  0
d
 i 1
n
ˆ 
1


t
i 1
n
i
ˆ  2.50 s
分布関数が指数関数の場合、
最尤値（推定値）は加算平均
例題の場合
（文献値は
2.20μs）
最尤法の誤差
•
確率分布関数を以下のように定義する
•
τの誤差は、確率分布関数から求められる以下の量で表す
– 信頼区間（Confidence Interval）
– 信頼度（Confidence Level）
具体例
•
規格化の
ため必要
演習
•
１１２ページに与えたμ粒子の崩壊時間分布から、以下の積分


0
•
L( )d 
を数値積分法（台形法またはSimpson法）により求めよ。
– 積分区間が[0,∞]であることに注意
– 0では被積分関数（指数関数部分）が発散
• 結果を見ながら、十分小さい数にする（１E-03より大きいこと）
– ∞は数値計算では不可能
• 結果を見ながら、十分大きな数にする（１E+04より小さいこと）
– 区分数は収束を見ながら決める
• ２１６より小さいこと（これ以上大きくすると計算時間がかかりすぎる）
確率分布関数
L( )


0
L(  )d 
を計算し、図示（gnuplot等を用いて）せよ
被積分関数
n
1
 ti 
L( )   e xp  
 
i 1  
非積分関数
の定義
被積分関数
1
 ti 
L( )   e xp  
 
i 1  
n
数値積分部分
積分の始点と終点
途中省略
途中省略
出力部分
指数表示
出力例
215でほぼ収束
確率分布関数のプログラム例
Likelihood
関数の積分値
時間
Likelihood関数
時間0.01μsから
10μsまで
0.1μs間隔で計算
出力は標準出力へ
（リダイレクションを使って
ファイルに出力させる）
結果の例
途中省略
時刻2.5μsで
Likelihood関数が最大
最小二乗法
•
•
最小二乗法とは
– 独立変数 xi に対して、測定データ yi が誤差 σi が n 個得られた場合
– 分布関数が分かっていない場合
• 最尤法が適用できない場合
– データに理論曲線を当てはめる手法として有効
データに k 次の多項式を当てはめる場合
k
f ( x)   a j x j
j 0
•
– データを最も良く再現する aj を求める
何をもって、「最も良く再現した」と判断するか？
– 関数 f(x) とデータ点 yi のズレを、カイ二乗
n
2  
i 1
( yi  f ( xi )) 2
 i2
で表し、この χ2が最小になる aj を求める
n
2( yi  f ( xi )) j
 2

xi  0 ( j  0, , k ) k次多項式の場合
2
a j i  1
i
k次多項式の場合の具体的表式
•
n
2( yi  f ( xi )) j
 2

xi  0 ( j  0, , k )
2
a j i  1
i
(k+1)次連立方程式
1次式（直線）の場合
•
•
係数a、bの誤差は、行列Aの逆行列A-1の対角要素で与えれる
– Data Reduction and Error Analysis for the Physical Science
P.R. Bevington, (McGraw-Hill, 1969)（興味のあるものは読んでみること）
A-1をどのようにもとめるか？
– ２行２列なので、式を書き下す方がよい（Gauss-Jordanを用いない）
最小二乗法と最尤法
•
•
最小二乗法で最小にしているカイ二乗の表式（天下り的に与えられた）はど
のようして決まっているのか？
– 二乗ではなく４乗でないのは何故か？
データがガウス分布で与えられる場合、最小二乗法と最尤法（尤度関数Lの
意味ははっきりしている）は等価である。
– xi における測定値の平均値（推定値）が f(xi ; a) で与えられる場合
ガウス分布
– 第１項目：変数 a に依存しない量（定数）
– 第２項目：χ2の-1/2
• χ2を最小にする ⇔ L を最大にする
再びμ粒子の崩壊（ヒストグラム分布の場合）
•
μ粒子は弱い相互作用により、ある寿命τで、
   e  e  
•
•
と崩壊する。
μ粒子の崩壊が指数分布
に従うとして、データから寿命τを推定し、その誤差を検討する。
例：μ粒子の崩壊数を0.5μs毎に測定
して、右表のようなデータを得た場合
（先程とは、データが異なる）
指数分布の場合の最小二乗法
•
当てはめたい関数（指数関数）
•
両辺のlnをとって
ln f ( t )  t  ln 
•
– １次式に帰着する
測定値のlnをとって、それに対して、１次式を当てはめればよい。
– lnをとっているので、測定値の誤差は
i 
i
yi
演習
•
μ粒子の寿命測定のデータを、
/home/teacher/z6wt01in/SAMPLE/mu_decay.dat
に置く。
– データは、
時間崩壊数崩壊数の誤差（崩壊数の平方根）
の並びである
• 最小二乗法のプログラムを書き、寿命τとその誤差を求めよ。
– サブルーチン”file_open”を用いてファイルを開く
• ５０～５１ページ参照
– ファイルから、データを配列に読み込む
• サブルーチン”read_data”を応用する
• ５１～５２ページ参照
– yiはln(yi)を代入する
– σiはσi/yiを代入する
– 具体的表式に従い、一次式の係数とその誤差を求める
メインプログラムの例
扱うデータ数の
最大数
データ
（x,y,δy）
を格納する配列
行列及び列ベクトル要素
０次の係数と誤差
行列式
１次の係数と誤差
メインプログラムの例（続き）
行列及び列ベクトルの初期化
この部分が
以前と異なる
データ
読み込み
行列及び
列ベクトルの
計算
メインプログラムの例（続き）
最小二乗法による
１次式の係数と
その誤差（逆行列要素）
結果の
出力
read_dataの変更の例
途中省略
log(y)をyの値とする
以下省略
log(y)の誤差はTaylor展開
１次でΔy/y
計算例
μ粒子の平均寿命
はこの逆数なので、
τ＝２．２μｓ

Download Report