狠狠撸

第 2 章例題: スペースシャトル「チャレンジャー号」
の悲劇
市東亘
西南学院大学経済学部
July 31, 2019
講義ノート: https://courses.wshito.com/semi2/2019-bayes-AI
市東亘 (西南学院大学経済学部) 第 2 章例題: スペースシャトル「チャレンジャー号」の悲劇 July 31, 2019 1 / 25

概観
概観
テキスト「Python で体験するベイズ推論」pp.63–74
例題: スペースシャトル「チャレンジャー号」の悲劇
? 分析の背後にある統計モデルを理解する．
? PyMC3 による実装方法を学ぶ．

分析の目的
分析の目的
目的
? 外気温が与えられた時の O リングの破損確率を推定したい．
? 標本データは 23 回のフライトデータから得られた外気温（華氏）と O
リング不良の有無．https://git.io/vXknD

分析の目的
説明したい事象と確率分布
? カンニングのケース
カンニング総数 ? 2 項分布
? チャレンジャーのケース
個別の O リングが破損するか否か ? ベルヌーイ分布
? ベルヌーイ分布
確率 p で 1 を，確率 1 ? p で 0 をとる離散確率分布．
f(x) = px(1 ? p)1?x for x = 1 or 0

分析の目的
統計モデルの構築
与えられるデータ
温度 O リング破損
この確率は外気温に依存すると考えられる．
ロジスティック関数
ロジスティック関数の形状は分からない
のでαとβも推定すべきパラメータとなる
ti Di
Di ～ Bernoulli(p)
Di ～ Bernoulli(p(ti))
p(ti) =
1
1 + e?(βti+α)

分析の目的
なぜロジスティック関数なのか？
? Di ～ Bernoulli(p(ti)) の p は 0 から 1 の値をとる確率．
? 温度 ti を 0 から 1 の値にマッピングする関数が必要．
? シグモイド曲線
6 4 2 0 2 4 6
0.0
0.2
0.4
0.6
0.8
1.0

分析の目的
ロジスティック曲線
温度が高い方が破損確率が高いケース
p =
1
1 + e?(βti+α)
β = 1, α = 0
6 4 2 0 2 4 6
0.0
0.2
0.4
0.6
0.8
1.0

分析の目的
温度が低い方が破損確率が高いケース
p =
1
1 + e?(βti+α)
β = ?1, α = 0
6 4 2 0 2 4 6
0.0
0.2
0.4
0.6
0.8
1.0

分析の目的
水平方向のシフト
p =
1
1 + e?(βti+α)
β = ?1, α = 50
40.0 42.5 45.0 47.5 50.0 52.5 55.0 57.5 60.0
0.0
0.2
0.4
0.6
0.8
1.0

分析の目的
変曲点の移動
p =
1
1 + e?(βx+α)
β = ?0.93, α = 50
40.0 42.5 45.0 47.5 50.0 52.5 55.0 57.5 60.0
0.0
0.2
0.4
0.6
0.8
1.0

分析の目的
統計モデル
f(α, β|Di, ti) ∝ f(Di, ti|α, β)f(α)f(β)
? f(Di, ti|α, β): Di ～ Bernoulli(p(ti))
where p(ti) = 1/(1 + e?(βti+α))
? f(α): α ～ N(0, 1000), ? = 0, σ2 = 1000
? f(β): β ～ N(0, 1000), ? = 0, σ2 = 1000
? Precision: τ = 1/σ2 = 1/1000 = 0.001

分析の目的
PyMC3 による実装
コード 1 データの読み込み
1 import numpy as np
2
3 # データは， Date,Temperature,Damage Incident のフォーマット．欠損値は NA が記されている
4 data = np.genfromtxt("challenger_data.csv", skip_header=1,
5 usecols=[1,2], missing_values="NA", delimiter=",")
6 data = data[~np.isnan(data[:, 1])]
7 data.shape
Out[1]: (23, 2)

分析の目的
PyMC3 による実装
コード 2 モデル: チャレンジャー号の悲劇
1 import numpy as np
2 import pymc3 as pm
3
4 temperature = data[:, 0] # 外気温
5 D = data[:, 1] # 破損の有無
6
7 with pm.Model() as model:
8 beta = pm.Normal("beta", mu=0, tau=0.001, testval=0) # βの事前分布
9 alpha = pm.Normal("alpha", mu=0, tau=0.001, testval=0) # αの事前分布
10 # ロジスティック関数
11 p = pm.Deterministic("p", 1.0/(1. + np.exp(-beta*temperature - alpha)))
12 D = pm.Bernoulli("D", p=p, observed=D) # 尤度
13
14 start = pm.find_MAP()
15 trace = pm.sample(60000, start=start) # 初期値を指定して MCMC を開始

分析の目的
生成サンプルの視覚化
pm.traceplot(trace["alpha"])
pm.traceplot(trace["beta"])

分析の目的
破損発生確率の視覚化

分析の目的
破損発生確率の視覚化
コード 3 破損発生確率の視覚化コード
1 temperature = data[:, 0]
2 D = data[:, 1]
3
4 # t は 50x1．p は 100000 x 50
5 t = np.linspace(temperature.min()-5, temperature.max()+5, 50)[:, None]
6 p = logistic(t.T, beta=trace["beta"][:, None], alpha=trace["alpha"][:, None])
7
8 mean_p = p.mean(axis=0)
9
10 fig, ax = plt.subplots(figsize=(12.5, 4))
11 ax.grid(True)
12 ax.plot(t, mean_p, lw=3, label="破損の平均事後確率")
13 ax.plot(t, p[100, :], ls="--", label="事後分布からのサンプル")
14 ax.plot(t, p[180, :], ls="--", label="事後分布からのサンプル")
15 ax.scatter(x=temperature, y=D, s=50, color="k", alpha=0.5)
16
17 ax.legend(loc="lower left")
18 ax.set_title("破損発生確率の事後期待値と 2 つのサンプリング値", fontsize=18)
19 ax.set_ylabel("確率", fontsize=16)
20 ax.set_xlabel("外気温", fontsize=16)

分析の目的
破損発生確率の視覚化コード解説
? 5 行目: グラフを 50 個の点を接続して描くために，外気温の最小値か
ら最大値まで 50 個の均等な値を変数 t に用意．
? 6 行目: MCMC で生成した 100000 個の α と β のサンプルデータ毎
に，50 個の t の値に対応した確率 p の値を計算する．p は
100000 × 50 の配列．
? 7 行目: 100000 個の確率データの平均を計算．
? 13–14 行目: 100000 個のサンプルデータからインデックス番号 100 と
180 の p の値をそれぞれ点線でプロット．

分析の目的
破損発生確率の 95%信用区間

分析の目的
コード 4 破損発生確率の 95%信用区間コード
1 from scipy.stats.mstats import mquantiles
2
3 qs = mquantiles(p, [0.025, 0.975], axis=0)
4 fig, ax = plt.subplots(figsize=(12.5, 4))
5 ax.grid(True)
6
7 ax.fill_between(t[:, 0], *qs, alpha=0.7, color="#7A68A6") # 第１引数は 1 次元でなけれ
8 ax.plot(t, qs[0], label="95%信用区間", color="#7A68A6", alpha=0.7) # 信用区間下側境界
9 ax.plot(t, qs[1], color="#7A68A6", alpha=0.7) # 信用区間上側境界線
10 ax.plot(t, mean_p, lw=1, ls="--", color="k", label="破損の平均事後確率")
11 ax.scatter(x=temperature, y=D, s=50, color="k", alpha=0.5)
12
13 ax.legend(loc="lower left")
14 ax.set_xlim(t.min(), t.max())
15 ax.set_title("破損発生確率の事後期待値と 95%信用区間", fontsize=18)
16 ax.set_ylabel("確率", fontsize=16)
17 ax.set_xlabel("外気温", fontsize=16)

分析の目的
破損発生確率の 95%信用区間コード解説
? mquantiles() は指定された割合にデータを分割し，その境界になる
データの値を返す関数．
? 3 行目: 引数 p は 100000 × 50 の配列．axis=0 はデータ p を列方向
に見ていく．したがって，50 の列ごとに 100000 個のサンプリング
データから 2.5%と 97.5%に入るデータの境界値（2 × 50 の配列）を
返す．
? 7 行目: fill_between() の引数のアスタリスクは配列のアンパッキン
グを表す．つまり 2 × 50 の配列が，50 個の要素を持つ 2 つの配列と
して渡される．

分析の目的
事故当時の外気温における破損確率
チャレンジャー号が打ち上げ失敗した日の気温華氏 31 度で，O リングが破
損する確率を推定してみる．

分析の目的
コード 5 破損発生確率の 95%信用区間コード
1 fig, ax = plt.subplots(figsize=(12.5, 2.5))
2 ax.grid(True)
3
4 # 華氏 31 度の時の事後確率サンプルデータを計算
5 prob_31 = logistic(31, beta=trace["beta"][:, None],
6 alpha=trace["alpha"][:, None])
7 # 1000 分割のヒストグラム．確率 0.001 ごとに頻度集計
8 ax.hist(prob_31, bins=1000, density=True, histtype="stepfilled")
9 ax.set_xlim([0.990, 1])
10 ax.set_ylabel("密度")
11 ax.set_xlabel("O リングで破損が起こる確率")
12 ax.set_title("外気温が華氏 31 度のときの破損が発生確率の事後分布")

分析の目的
事故当時の外気温における破損確率コード解説
? 5 行目: 100000 個の α と β のサンプルデータを利用して，事故があっ
た日の気温華氏 31 度における O リング破損の事後確率のサンプル
データを生成．
? 8 行目: この 100000 個の事後確率サンプルデータを 1000 個の区間に
分割し，その頻度からヒストグラムを作成．確率 0 から 1 まで 1000 分
割すると 0.001 毎にグラフが描かれる．
? 9 行目: ヒストグラムは 0.990 から 1 までのデータのみ描画．

狠狠撸

2019年演習II．第2章例題: スペースシャトル「チャレンジャー号」の悲劇

More Related Content

2019年演習II．第2章例題: スペースシャトル「チャレンジャー号」の悲劇

狠狠撸

2019年 演習II．第2章 例題: スペースシャトル「チャレンジャー号」の悲劇

More Related Content

2019年 演習II．第2章 例題: スペースシャトル「チャレンジャー号」の悲劇

2019年演習II．第2章例題: スペースシャトル「チャレンジャー号」の悲劇

2019年演習II．第2章例題: スペースシャトル「チャレンジャー号」の悲劇