モンテカルロ法(2) 誤差解析と不偏推定量 / Simulation 03

シミュレーション工学

Transcript

1. 1 64 モンテカルロ法(2) 誤差解析と不偏推定量 シミュレーション工学 慶應義塾大学大学院理工学研究科基礎理工学専攻物理情報専修 渡辺宙志
2. 2 64 • エラーバーとは何か、どのような性質を持つ かを理解する • 統計誤差と系統誤差について理解する • 系統誤差を除去する(Jackknife法) 本講義の目的 測定と誤差 • 一般に測定値には実験誤差がある • 数値計算においても、測定結果は誤差を伴う • 「誤差」の理解は難しい
3. 3 64 ある回路の電流を6回測定したら、以下のデータ[mA]を得た 1.16, 1.13, 1.12, 1.12, 1.11, 1.08 この回路の電流の観測値は? 1.12 ± 0.01 mA 観測値 エラーバー 大 雑 把 な 意 味 ： 観測値の1.1 ま で は 自 信 があるが、 小 数 点 第 二 位 は 自 信 がな く 、1.11かもし れ ないし1.13かもし れ ない 正確な定義は？
4. 4 64 観測するたびに値が変化する量を確率変数 ෠ 𝑋とみなす この変数の「全ての可能性の集合」を母集団と呼ぶ 観測により母集団から標本を取り出す 標本 {𝑋𝑖 } 標本の集合から母集団の性質を推定するのが目的 母集団 { ෠ 𝑋}
5. 5 64 母集団の特徴量(平均や分散)を知りたい 𝑓(𝑋) ത 𝑋 𝜎 ത 𝑋 = ∫ 𝑋𝑓 𝑋 𝑑𝑋 分布の1次のモーメント 𝜎2 = ∫ 𝑋 − ത 𝑋 2𝑓 𝑋 𝑑𝑋 手元にあるのは𝑁個の標本{𝑋𝑖 } 標本から特徴量を得る関数を推定量(estimator)と呼ぶ 分布の2次のモーメント
6. 6 64 平均値 ത 𝑋 = 1 𝑁 ෍ 𝑖 𝑁 𝑋𝑖 母分散 𝜎2 = 1 𝑁 − 1 ෍ 𝑖 𝑁 𝑋𝑖 − ത 𝑋 2 平均値の 推定値の分散 𝜎ത 𝑋 2 = 1 𝑁(𝑁 − 1) ෍ 𝑖 𝑁 𝑋𝑖 − ത 𝑋 2 推定したい量 そのestimator ※ N で はな く N-1 で 割るのは不偏分散を求めるため
7. 7 64 1.12 ± 0.01 mA 0.01 平均値の推定値の分散の平方根をエラーバーとする ത 𝑋 ± 𝜎ത 𝑋 2 エラーバーの 意 味 は？ 平均値の推定値の標準偏差を誤差とみなす
8. 8 64 • 一般に、観測値の分布はガウス分布 で はない • しかし、観測値が独立同分布に従う確率変数と みなせる時、その期待値はガウス分布に近づ く サイコロの目の母集団の分布は 一様分布だが、 千回振った目の平均値の分布は ガウス分布に近づ く
9. 9 64 標本が多 く なるほど平均値の分布の分散は 小 さ く な る N=10, 100, 1000回サイコロを振った時、出た目の平均の頻度分布
10. 10 64 𝜎2 = 1 𝑁 − 1 ෍ 𝑖 𝑁 𝑋𝑖 − ത 𝑋 2 サンプル数𝑁を増やした時 lim 𝑁→∞ 𝜎2 = const. 母分散のestimatorは一定値に収束する 平均値の推定値の分散のestimatorは1/𝑁の早さ で ゼ ロに収束する 𝜎ത 𝑋 2 = 𝜎2 𝑁 lim 𝑁→∞ 𝜎ത 𝑋 2 = 0
11. 11 64 𝑓(𝑥) 𝜎 𝑥 𝑓 𝑥 = 1 2𝜋𝜎2 exp − 𝑥 − 𝜇 2 2𝜎2 𝜇 平均𝜇、分散𝜎2のガウス分布 平均𝜇： 分布の中心の 位 置 標準偏差𝜎： 分布の幅
12. 12 64 𝑃 𝑎 < ෠ 𝑋 < 𝑏 = න 𝑎 𝑏 𝑓 𝑥 𝑑𝑥 確率変数 ෠ 𝑋の値がaとbの間にある確率が で 与 えら れ る 時、𝑓 𝑥 を ෠ 𝑋の確率密度関数と呼ぶ 確率密度関数が平均𝜇、分散𝜎2のガウス分布 で ある時 𝜇 − 𝜎 < 𝑥 < 𝜇 + 𝜎 の範囲を1シグマの範囲と呼ぶ
13. 13 64 𝑓 𝑥 = 1 2𝜋𝜎2 exp 𝑥 − 𝜇 2 2𝜎2 で あ るとき、 𝑃 𝜇 − 𝜎 < ෠ 𝑋 < 𝜇 + 𝜎 = න 𝜇−𝜎 𝜇+𝜎 𝑓 𝑥 𝑑𝑥 ∼ 0.6827 平均𝜇、分散𝜎2のガウス分布に従う確率変数が、 平均の周りに𝜎の間 で 揺らぐ確率が68.27% ex) テスト で 偏差値40～60 ま で の間の人が68.27%
14. 14 64 0.01 エラーバーを平均値の推定値の標準偏差とする(1シグマの範囲) 同様な実験を繰り返した場合、観測値が エラーバーの間に入る確率が68.27%
15. 15 64 𝜇 − 𝑛𝜎 < 𝑥 < 𝜇 + 𝑛𝜎 同様に「nシグマの範囲」が定義 で きる 1シグマ ： 入る確率 68.27% 外 れ る確率 31.73% 2シグマ ： 入る確率 95.45% 外 れ る確率 4.55% 3シグマ ： 入る確率 99.73% 外 れ る確率 0.27% 5シグマ ： 入る確率 99.9994% 外 れ る確率 0.0005%
16. 16 64 • エラーバーとは観測値を確率変数とみなした時に、 その平均値の分布の推定標準偏差のこと • サンプル数を増やせば増やすほど、エラーバーは 小 さ く なる • 観測値が独立同分布なら、サンプル数を増やして い く と平均値の分布はガウス分布に漸近する • 平均𝜇、分散𝜎2のガウス分布について、以下を「n シグマの範囲」と呼ぶ • ガウス分布に従う確率変数が独立 で あ るなら • 「1シグマの範囲」からは3つに1つは外 れ る • 「5シグマの範囲」から外 れ る確率はほぼゼロ 𝜇 − 𝑛𝜎 < 𝑥 < 𝜇 + 𝑛𝜎
17. 17 64 データがガウス分布に従い、かつ独立 で あるとする 観測量の母集団の分布の平均を「真の値」と呼ぶと • 観測値は「真の値」の上下に均等にばらつ く • 観測値の3つに1つが「真の値」の1シグマの 範囲に入らない • 観測値と「真の値」がエラーバーの2倍離 れ る ことは稀、5倍離 れ ることは ま ずない この知識を活用して「おかしなグラフ」に気づ く ことが で き る
18. 18 64 何かが指数関数的に減衰しているようだが・・・？
19. 19 64 なんとな く こ んな線が見える 計算精度を高 く していったら、データはこの線に収束する で あろうと期待さ れ る線→「真の値」
20. 20 64 もしエラーバーが1シグマの範囲 で 取ら れ ていたら 3つに1つは「真の値」から外 れ ないとおかしい 全てのデータ 点 について 「真の値」にエラーバーがかかっている
21. 21 64 import numpy as np N = 10 np.random.seed(1) for i in range(10): x = i + 0.5 d = np.zeros(N) d += np.exp(-x/3) d += np.random.randn(N)*0.1 y = np.average(d) e = np.std(d) print(f"{x} {y} {e}") 先ほどのデータを生成したコード エラーバーとしてnumpy.stdを その ま ま 使っている 平均値の推定誤差 で はな く 、母集団の標準偏差を求めてし ま っている
22. 22 64 import numpy as np N = 10 np.random.seed(1) for i in range(10): x = i + 0.5 d = np.zeros(N) d += np.exp(-x/3) d += np.random.randn(N)*0.1 y = np.average(d) e = np.std(d)/np.sqrt(N) print(f"{x} {y} {e}") 𝜎 ത 𝑋 = 𝜎 𝑁 こ れ が正しいコード
23. 23 64 線から外 れ て いるデータがある 1シグマの範囲なら「外 れ ているデータ」がないと不 自 然
24. 24 64 ある観測値のサンプル数n依存性のグラフ サンプル数が増えると収束し、かつエラーバーが 小 さ く な るのは もっともらしいが・・・？
25. 25 64 「同じ側」に外 れ てることが 続いている 各データ 点 が 独立なら、「真の値」の両側にばらつ く はず 「真の値」はこの あたりにありそう
26. 26 64 import numpy as np np.random.seed(1) N = 2048 d = np.random.random(N) for i in range(4, 12): n = 2**i dd = d[:n] ave = np.average(dd) err = np.std(dd)/np.sqrt(n) print(f"{n} {ave} {err}") 先ほどのデータを生成したコード 先に全データを作成し、部分配列に ついて誤差を計算している
27. 27 64 異なるデータ 点 で 共通するデータを使っている →データ 点 が独立 で はない
28. 28 64 import numpy as np np.random.seed(1) N = 2048 for i in range(4, 12): n = 2**i dd = np.random.random(n) ave = np.average(dd) err = np.std(dd)/np.sqrt(n) print(f"{n} {ave} {err}") データ 点 ごとに異なる データセットを生成している データを適切に生成するコード
29. 29 64 「真の値」の両側に均等にばらついており、もっともらしい
30. 30 64 よほど複 雑 なデータ で ない限り、ゼロの ま わりを揺らぐデータ に見えるが・・・？
31. 31 64 明らかに5シグマ以上離 れ ている
32. 32 64 import numpy as np N = 1000 v = 0.0 gamma = 0.1 np.random.seed(1) for j in range(10): d = np.zeros(N) for i in range(N): v += np.random.randn()*0.1 v -= gamma * v d[i] = v ave = np.average(d) err = np.std(d) / np.sqrt(N) print(f"{(j+0.5)*N} {ave} {err}") ランジュバン方程式の数値解法 𝑚 𝑑𝑣 𝑑𝑡 = −𝛾𝑣 + ෠ 𝑅 𝑣 1000ステップごとに平均、標準偏差を計算するコード
33. 33 64 速度の時間発展データ 相関時間 速度が「記憶」を失う ま で にそ れ なりの時間がかかる
34. 34 64 「記憶」を忘 れ そうな時間をあけて観測してサンプリングする ※もっとかしこい方法もある
35. 35 64 データのばらつき具合、エラーバーの外 れ 具合、ともにもっともらしい
36. 36 64 エラーバーが大きすぎる 適切なグラフ 原因の例 • 𝑁で割り忘 れ て いる • データに相関がある
37. 37 64 偏りが大きい 適切なグラフ 原因の例 • データに相関がある
38. 38 64 エラーバーが 小 さすぎる 適切なグラフ 原因の例 • データに相関がある
39. 39 64 データがガウス分布に従い、かつ独立 で あるなら • 観測値は「真の値」の上下に均等にばらつ く • 観測値の3つに1つが「真の値」の1シグマの範囲に入らな い • 観測値と「真の値」がエラーバーの2倍離 れ ることは稀、5倍 離 れ ることは ま ずない 逆に • 観測値の全てが「真の値」をエラーバーの範囲に含む • 「真の値」の片側に連続してず れ て いる • 「真の値」と5シグマ以上離 れ ている で あ るなら、何かがおかしい エラーバーがおかしいグラフは、データの相関が原因 で あることが多い
40. 40 64 M個のデータがある そ れ をN個ずつのブロックに分割する 𝑁 𝜇𝑀/𝑁 𝜇1 𝜇2 そ れ ぞ れ のブロック で 期待値を計算する 期待値の期待値を計算する 𝜇 = 𝑁 𝑀 ෍ 𝑖 𝜇𝑖 𝑀 ブロックサイズ𝑁を変えた時 𝜇 は変わるか？
41. 41 64 𝜇 = 𝑁 𝑀 ෍ 𝑖 𝜇𝑖 𝜇𝑖 = 1 𝑁 ෍ 𝑘∈𝑖 𝑋𝑘 ブロックごとの期待値 期待値の期待値 𝜇 = 1 𝑀 ෍ 𝑗 𝑋𝑗 単なる全体の平均になる 𝜇 は𝑁依存性をもたない
42. 42 64 𝜇1 𝜇2 そ れ ぞ れ のブロック で 期待値を計算する ⋯ 1/𝜇1 ⋯ 1/𝜇2 そ れ ぞ れ のブロックの期待値の逆数を計算する 期待値の逆数の期待値を計算する 1/𝜇 = 𝑁 𝑀 ෍ 𝑖 1 𝜇𝑖 1/𝜇 は𝑁依存性を持つか？
43. 43 64 サイコロの目の期待値の逆数は？ 期待値 𝜇 = 1 6 ෍ 𝑘=1 6 𝑘 = 7 2 = 3.5 期待値の逆数 1 𝜇 = 2 7 ∼ 0.286
44. 44 64 サイコロを65536回振る 4個ずつ分割 8個ずつ分割 16個ずつ分割 ・・・ ・・・ ・・・ ・・・ 各ブロック で 期待値𝜇𝑖 を計算し、その逆数の期待値を計算する
45. 45 64 明らかに𝑁依存性がある 厳密解 同じデータセットを使っているのに、ブロックサイズが 小 さいところ で 挙動がおかしい→系統誤差 ず れ ているわりには エラーバーが 小 さい
46. 46 64 誤差(真値からのず れ )には、統計誤差と系統誤差の 二 種類がある 統計誤差 (statistical error) • 我々が制御 で きない要因により値が揺らぐこと(偶然誤差) • 数値計算 で は 乱数や粗視化に起因 • 不確かさ(uncertainty)とも 系統誤差 (bias, systematic error) • 誤差を生む要因が説明 で きるもの • 決定論的なず れ • 数値計算 で は 有限サイズ効果や理論誤差などに起因
47. 47 64 このデータがず れ て いるのは偶然 統計誤差 このデータがず れ て いるのは必然 系統誤差 • この系統誤差はどこから く るのか？ • どうやって減らすか を知るのがこの節の目的
48. 48 64 N回測定して期待値を推定する(こ れ も確率変数) 確率変数 ෠ 𝑋の期待値𝜇の関数の値を推定したい 𝑦 = 𝑔(𝜇) ො 𝜇𝑁 = 1 𝑁 ෍ 𝑖 𝑋𝑖 推定値の期待値は期待値に一致する 𝑔(ො 𝜇𝑁 ) ≠ 𝑔 𝜇 Ƹ 𝜇𝑁 = 𝜇 推定値の関数の期待値は期待値の関数と一致しない
49. 49 64 標本から得ら れ た推定量(estimator)の期待値が 母集団の期待値と一致する時、その推定量を 不偏推定量(unbiased estimator)と呼ぶ ො 𝜇𝑁 = 1 𝑁 ෍ 𝑖 𝑁 𝑋𝑖 例 ： 確率変数 ෠ 𝑋の𝑁個のサンプル{𝑋𝑖 }から母集団{ ෠ 𝑋}の期 待値𝜇と分散𝜎2を求めたい ො 𝜎𝑁 2 = 1 𝑁 ෍ 𝑖 𝑁 (𝑋𝑖 −ො 𝜇𝑁 ) Ƹ 𝜇𝑁 = 𝜇 期待値は一致する(不偏推定量) ො 𝜎𝑁 2 = 𝑁 − 1 𝑁 𝜎2 ≠ 𝜎2 分散は一致しない(不偏推定量 で はない) 期待値の推定値を 使っているのがポイント
50. 50 64 一般に確率変数 ෠ 𝑋 について 𝑔( ෠ 𝑋) ≠ 𝑔 ෠ 𝑋 と 𝑔( ෠ 𝑋) 関数の期待値 期待値の関数 𝑔 ෠ 𝑋 は 一致しない 期待値の関数は、期待値の関数の不偏推定量 で はない
51. 51 64 μ 𝑔(𝑥)を上に凸な関数とし、𝑥 = 𝜇で 接 線 をひ く 𝑦 = 𝑎 𝑥 − 𝜇 + 𝑔(𝜇) 𝑦 = 𝑔(𝑥) ※ Thanks to @genkuroki 上図より明らかに 𝑔 𝑥 ≤ 𝑎 𝑥 − 𝜇 + 𝑔(𝜇) 両辺の期待値を取 れ ば 𝑔 𝑥 ≤ 𝑔 𝜇 = 𝑔( 𝑥 ) 下に凸の場合は符号が逆になる
52. 52 64 𝜀 = Ƹ 𝜇𝑁 − 𝜇 Ƹ 𝜇𝑁 = 1 𝑁 ෍ 𝑖 ෠ 𝑋𝑖 N回の測定 で 得ら れ た期待値の推定量 真の期待値とのず れ 𝑔 Ƹ 𝜇𝑁 − 𝑔 𝜇 = 𝑔 𝜇 + 𝜀 − 𝑔 𝜇 = 𝑔′ 𝜇 𝜀 + 1 2 𝑔′′ 𝜇 𝜀2 + 𝑂(𝜀3) 𝑔 Ƹ 𝜇𝑁 − 𝑔 𝜇 ∼ 1 2 𝑔′′ 𝜇 𝜀2 = 𝑔′′(𝜇)𝜎2 2𝑁 真の値 推定値 推定値と真の値のず れ の期待値 N依存性 期待値の推定値の分散
53. 53 64 𝑔 Ƹ 𝜇𝑁 − 𝑔 𝜇 ∝ 1 𝑁 𝑁個 のサンプルから推定した期待値の関数と、 真の期待値の関数のず れ は 1/𝑁に比例する こ れ を1/Nバイアスと呼ぶ 𝑔 Ƹ 𝜇𝑁 − 𝑔 𝜇 ∼ 𝑔′′(𝜇)𝜎2 2𝑁 関数𝑔 𝑥 の 二 階微分がゼロ(線形) で ある場合はバイアスは生じない 特に𝑔 𝑥 = 𝑥の場合 Ƹ 𝜇𝑁 = 𝜇
54. 54 64 平均0、分散𝜎2のガウス分布に従う確率変数Xを考える ෠ 𝑋2 = 𝜎2 ෠ 𝑋4 = 3𝜎4 2次のモーメント 4次のモーメント 4次と2次のモーメントの比を取ると、分散依存性が消える ෠ 𝑋4 ෠ 𝑋2 2 = 3 尖度(Kurtosis) この量の1/Nバイアスを確認する
55. 55 64 ෠ 𝑋2 𝑁 = 1 𝑁 ෍ 𝑖 ෠ 𝑋𝑖 2 ෠ 𝑋4 𝑁 = 1 𝑁 ෍ 𝑖 ෠ 𝑋𝑖 4 N個のサンプリング(N回の測定) で 得ら れ たデータから 2次と4次のモーメントを推定する ෡ 𝑈𝑁 = ෠ 𝑋4 𝑁 ෠ 𝑋2 𝑁 2 得ら れ たモーメントから尖度を計算する 上記を十分に繰り返して෡ 𝑈𝑁 の期待値 ෡ 𝑈𝑁 を計算する 統計誤差を消し、系統誤差だけを残す
56. 56 64 𝑈𝑁 1/𝑁 理論値 十分なサンプリング回数にも関わらず、真の値からず れ ている(バイアス) 推定値
57. 57 64 不偏推定量 で はあるが、ばらつきのせい で 真の値から ず れ る誤差を統計誤差と呼ぶ Ƹ 𝜇𝑁 = 1 𝑁 ෍ 𝑖 ෠ 𝑋𝑖 Ƹ 𝜇𝑁 − 𝜇 = 𝑂(1/ 𝑁) 不偏推定量 で ない推定量の期待値について、真の値か らのず れ を系統誤差(バイアス)と呼ぶ。 𝑔( Ƹ 𝜇𝑁 ) − 𝑔 𝜇 = 𝑂(1/𝑁) サンプル数を増やすと統計誤差は減るが、 系統誤差は減らせない
58. 58 64 期待値の関数の推定には1/Nバイアスが乗る N無限大極限 で は一致するが、収束が遅い 手持ちのデータから1/Nバイアスを除去したい Jackknifeリサンプリング
59. 59 64 N個のデータがある 𝑁 全部のデータを使って期待値𝜇𝑁 を計算 そ れ を使って関数の推定値𝑈𝑁 = 𝑔(𝜇𝑁 )を計算 1個のデータを捨てる 𝑁 − 1 残りのデータを使って期待値𝜇𝑁−1 を計算 そ れ を使って関数の推定値𝑈𝑁−1 = 𝑔(𝜇𝑁−1 )を計算
60. 60 64 𝑈𝑁 は、真の値𝑈∞ に対して1/Nバイアスがあると仮定 𝑈𝑁 = 𝑈∞ + 𝑎/𝑁 一つデータを捨てて得た𝑈𝑁 のバイアスは 𝑈𝑁−1 = 𝑈∞ + 𝑎/(𝑁 − 1) この2式から𝑈∞ を求めると 𝑈∞ = 𝑁𝑈𝑁 − (𝑁 − 1)𝑈𝑁−1 ※ Thanks to smorita and yomichi
61. 61 64 𝑈𝑁 1/𝑁 𝑁 = ∞ NとN-1から1/N→0外挿を行った
62. 62 64 1個のデータ除外して計算 せっか く のデータを捨てるのはもったいないの で 活用する 𝑈𝑁−1 1 𝑈𝑁−1 2 別のデータ除外して計算 ・ ・ ・ 𝑈𝑁−1 𝑁 𝑈𝑁−1 = 1 𝑁 ෍ 𝑖 𝑈𝑁−1 𝑖 精度の高い「N-1個のデータの推定量」 が得ら れ る
63. 63 64 𝑈𝑁 1/𝑁 理論値 単純な推定値 Jackknifeによるバイアス除去 𝑁𝑈𝑁 − (𝑁 − 1)𝑈𝑁−1 𝑈𝑁
64. 64 64 • 母集団の何かを推定する量を推定量(estimator)と呼ぶ • 誤差には統計誤差と系統誤差(バイアス)がある • その期待値が母集団の期待値に一致する量(バイアスが無い 量)を不偏推定量(unbiased estimator)と呼ぶ • 期待値の関数の単純な推定は不偏推定量を与えない • リサンプリングによりバイアスを除去 で きる • Jackknife法はリサンプリング法の一種