AIエンジニア/データサイエンティストが使うベイズ統計モデル
Python3 PyMC3 によるMCMC（Markov chain Monte Carlo）
スペースシャトル「チャレンジャー号」の悲劇の原因はＯリング
頻度主義統計学ともいう非ベイズ統計学では解けない問題をベイズモデルとMCMCで解く。
1986年1月28日，米国はスペースシャトル「チャレンジャー号」の打ち上げに失敗した。原因は，ロケットブースターの溶接部のＯリングの破損であると結論づけられた。
スペースシャトルの過去の打ち上げデータから，そのＯリングの破損状況とその日の外気温のデータがある。このＯリングは過去に幾度も破損しているので，このＯリングの破損がチャレンジャー号の爆発の原因であるとは，別の工学的ストーリーが必要である。
ここでの関心事は，その日の外気温によりＯリングが破損する確からしさを求めることができるか，換言すれば，気温$t$を変数とした，破損発生確率$p(t)$を求めることである。
In [1]:
from IPython.display import Image
Image("./image_common/book_Davidson-Pilon.jpg")
Out[1]:
『Pythonで体験するベイズ推論』 ーPyMCによるMCMC入門ー キャメロン デビッドソン=ピロン (著), 玉木 徹 (翻訳)
In [2]:
%matplotlib inline
import pymc3 as pm
import numpy as np
import theano.tensor as tt
from IPython.core.pylabtools import figsize
import matplotlib.pyplot as plt
import scipy.stats as stats
破損発生と外気温の関係
In [3]:
figsize(12.5, 3.5)
np.set_printoptions(precision=3, suppress=True)
challenger_data = np.genfromtxt("data/challenger_data.csv", skip_header=1,
usecols=[1, 2], missing_values="NA",
delimiter=",")
#drop the NA values
challenger_data = challenger_data[~np.isnan(challenger_data[:, 1])]
#plot it, as a function of tempature (the first column)
print("Temp (F), O-Ring failure?")
print(challenger_data)
plt.scatter(challenger_data[:, 0], challenger_data[:, 1], s=75, color="k",
alpha=0.5)
plt.yticks([0, 1])
plt.ylabel("Damage Incident?")
plt.xlabel("Outside temperature (Fahrenheit)")
plt.title("Defects of the Space Shuttle O-Rings vs temperature");
ロジスティック関数
外気温が下がると，破損発生の確率は明らかに上昇している。外気温と破損発生の間には明確なしきい値があるようには見えない。そこで，気温$t$を変数とした，破損発生確率$p(t)$を考える。
気温の横軸が上昇するにつれて，確率が1から0に滑らかに変化するのは，二値分類に使うロジスティック回帰のロジスティック関数がよい。
$$p(t) = \frac{1}{ 1 + e^{ \;\beta t + \alpha } } $$パラメータ$\alpha$と$\beta$を変えて，ロジスティック関数をプロットする。
パラメータ$\alpha$を変えると，横に移動するので，バイアスと呼ばれている。
In [4]:
def logistic(x, beta, alpha=0):
return 1.0 / (1.0 + np.exp(np.dot(beta, x) + alpha))
x = np.linspace(-4, 4, 100)
plt.plot(x, logistic(x, 1), label=r"$\beta = 1$", ls="--", lw=1)
plt.plot(x, logistic(x, 3), label=r"$\beta = 3$", ls="--", lw=1)
plt.plot(x, logistic(x, -5), label=r"$\beta = -5$", ls="--", lw=1)
plt.plot(x, logistic(x, 1, 1), label=r"$\beta = 1, \alpha = 1$",
color="#348ABD")
plt.plot(x, logistic(x, 3, -2), label=r"$\beta = 3, \alpha = -2$",
color="#A60628")
plt.plot(x, logistic(x, -5, 7), label=r"$\beta = -5, \alpha = 7$",
color="#7A68A6")
plt.legend(loc="lower left");
パラメータ$\alpha$と$\beta$には，正規分布
パラメータ$\alpha$と$\beta$には，制約がないので，正規分布がよい。
$$ f(x | \mu, \tau) = \sqrt{\frac{\tau}{2\pi}} \exp\left( -\frac{\tau}{2} (x-\mu)^2 \right) $$$\mu$ と $\tau$ を変えて，正規分布をプロットしてみる。
In [5]:
import scipy.stats as stats
nor = stats.norm
x = np.linspace(-8, 7, 150)
mu = (-2, 0, 3)
tau = (.7, 1, 2.8)
colors = ["#348ABD", "#A60628", "#7A68A6"]
parameters = zip(mu, tau, colors)
for _mu, _tau, _color in parameters:
plt.plot(x, nor.pdf(x, _mu, scale=1./_tau),
label="$\mu = %d,\;\\tau = %.1f$" % (_mu, _tau), color=_color)
plt.fill_between(x, nor.pdf(x, _mu, scale=1./_tau), color=_color,
alpha=.33)
plt.legend(loc="upper right")
plt.xlabel("$x$")
plt.ylabel("density function at $x$")
plt.title("Probability distribution of three different Normal random \
variables");
正規分布の期待値
$$ E[ X | \mu, \tau] = \mu$$正規分布の分散
$$Var( X | \mu, \tau ) = \frac{1}{\tau}$$ベイズモデルの構築
次にMCMCのモデリングをします。MCMCの典型的な書き方です。
withでモデル構築をします。as modeでモデルに名前を付けて，with model:でモデルの構築を続けます。
尤度はベルヌーイ分布
コイントスと同じ2項分布の1回投げがベルヌーイ分布になります。
$$ \text{Defect Incident, $D_i$} \sim \text{Ber}( \;p(t_i)\; ), \;\; i=1..N$$observed = がモデルの出力で，その右辺が尤度です。その第2引数が事前分布，第3引数が入力データです。
尤度の第2引数$p$を出力しているのが，事前分布のモデルです。$\alpha$と$\beta$を正規分布から出力して，それをロジスティック関数に入力しています。ロジスティック関数はベタなスクラッチで書かれています。
pm.sampleがMCMCを実行しています。trace = は事後分布のサンプルです。その次の行では事後分布のサンプルをスライスしています。
pm.sampleのstartは，この場合には，pm.find_MAP()であり，初期値をMAP推定法(最大事後確率推定法)で決めています（詳細省略）。
pm.sampleのstepは，MCMCの種類であり，この場合には，Metropolis法を選んでいます。
In [6]:
import pymc3 as pm
temperature = challenger_data[:, 0]
D = challenger_data[:, 1] # defect or not?
#notice the`value` here. We explain why below.
with pm.Model() as model:
beta = pm.Normal("beta", mu=0, tau=0.001, testval=0)
alpha = pm.Normal("alpha", mu=0, tau=0.001, testval=0)
p = pm.Deterministic("p", 1.0/(1. + tt.exp(beta*temperature + alpha)))
In [7]:
# connect the probabilities in `p` with our observations through a
# Bernoulli random variable.
with model:
observed = pm.Bernoulli("bernoulli_obs", p, observed=D)
# Mysterious code to be explained in Chapter 3
start = pm.find_MAP()
step = pm.Metropolis()
trace = pm.sample(120000, step=step, start=start)
burned_trace = trace[100000::2]
$\alpha$ と $\beta$ の事後分布
In [8]:
alpha_samples = burned_trace["alpha"][:, None] # best to make them 1d
beta_samples = burned_trace["beta"][:, None]
figsize(12.5, 6)
#histogram of the samples:
plt.subplot(211)
plt.title(r"Posterior distributions of the variables $\alpha, \beta$")
plt.hist(beta_samples, histtype='stepfilled', bins=35, alpha=0.85,
label=r"posterior of $\beta$", color="#7A68A6", normed=True)
plt.legend()
plt.subplot(212)
plt.hist(alpha_samples, histtype='stepfilled', bins=35, alpha=0.85,
label=r"posterior of $\alpha$", color="#A60628", normed=True)
plt.legend();
破損の事後確率の期待値
In [9]:
t = np.linspace(temperature.min() - 5, temperature.max()+5, 50)[:, None]
p_t = logistic(t.T, beta_samples, alpha_samples)
mean_prob_t = p_t.mean(axis=0)
In [10]:
figsize(12.5, 4)
plt.plot(t, mean_prob_t, lw=3, label="average posterior \nprobability \
of defect")
plt.plot(t, p_t[0, :], ls="--", label="realization from posterior")
plt.plot(t, p_t[-2, :], ls="--", label="realization from posterior")
plt.scatter(temperature, D, color="k", s=50, alpha=0.5)
plt.title("Posterior expected value of probability of defect; \
plus realizations")
plt.legend(loc="lower left")
plt.ylim(-0.1, 1.1)
plt.xlim(t.min(), t.max())
plt.ylabel("probability")
plt.xlabel("temperature");
95%信用区間（95%CI ; credible interval）
In [11]:
from scipy.stats.mstats import mquantiles
# vectorized bottom and top 2.5% quantiles for "confidence interval"
qs = mquantiles(p_t, [0.025, 0.975], axis=0)
plt.fill_between(t[:, 0], *qs, alpha=0.7,
color="#7A68A6")
plt.plot(t[:, 0], qs[0], label="95% CI", color="#7A68A6", alpha=0.7)
plt.plot(t, mean_prob_t, lw=1, ls="--", color="k",
label="average posterior \nprobability of defect")
plt.xlim(t.min(), t.max())
plt.ylim(-0.02, 1.02)
plt.legend(loc="lower left")
plt.scatter(temperature, D, color="k", s=50, alpha=0.5)
plt.xlabel("temp, $t$")
plt.ylabel("probability estimate")
plt.title("Posterior probability estimates given temp. $t$");
破損が発生する事後確率
事故が発生した当日の外気温が華氏31度だった。このときの破損が発生する事後確率をプロットする。
Ｏリングが破損する確率は非常に高かったと言わざるを得ないです。
In [12]:
figsize(12.5, 2.5)
prob_31 = logistic(31, beta_samples, alpha_samples)
plt.xlim(0.995, 1)
plt.hist(prob_31, bins=1000, normed=True, histtype='stepfilled')
plt.title("Posterior distribution of probability of defect, given $t = 31$")
plt.xlabel("probability of defect occurring in O-ring");
