tak0kadaの何でもノート

発声練習、生存確認用。

医学関連は 医学ノート へ

「Probabilistic Programming and Bayesian Methods for Hackers」を読んだ - (0 ~ 2章)

読書メモ 統計 数学 python MCMC

MCMCの勉強のためにProbabilistic Programming and Bayesian Methods for Hackersを読んだ。翻訳は終わっていないが、日本語版もある。

第0章

ベイズ推定を数理的にとくのは難しいが、コンピュータの計算力を使うと簡単に解けるYOという話。計算の際にはPyMCを使う(RでのWinBUGSやrstanに相当する)。

章の構成は、

  • 第1章: ベイズ推定入門
  • 第2章: もう少しPyMCについて
  • 第3章: MCMCの扉を開く
  • 第4章: 語られぬ偉大な定理
  • 第5章: 手足を失う損失は?
  • 第6章: Getting our prior-ities straight

のようになっている。

第1章

第1章はイントロダクション的な話を中心にしていて、ベイズ統計と頻度主義の違いについてや、確立分布についての説明、PyMCの使用例がある(第2章でまとめる)。

ベイズ統計での確率は個人の「信念」に基づいていることが特徴。ベイズ主義でもデータの量が増えると頻度主義と同様の値に収束していく。頻度主義の考え方は誤っているわけではなく、最小二乗回帰、LASSO回帰、最尤推定は高速で強力な手法である。ベイズ的な方法論は自由なモデリングができるという利点を活かして上のような手法で解決できない問題を解くのに使う。また、ベイズ主義が必要でなくなるほどデータがたくさんあって余るということはない。たくさんあっても結局は更に細かい問題を解くことになるからである。

第2章

例1 メール数問題

f:id:kuyata:20140818204855p:plain

「自分にやってくるメールの数を数えてデータをプロットしたところ、前半と後半でメールの数が違うように感じる。どうモデリングすればいいだろうか?」という問題について考える。

  • メールの数は単位時間あたりのイベントの生起確率を表すポアッソン分布に従うと考える。
    • i日目のメールの数を C_{i} = \text{Poisson}(\lambda)と表すことにする。
  • この \lambdaは観測できない値であるが、この \lambdaがある日を境にして変化したと考えると、
    • $$\lambda = \begin{cases} \lambda_{1} &\text{if } t \lt τ\\ \lambda_{2} &\text{if } t \geq τ \end{cases}$$
  • ベイズ推定を行うにはパラメータの事前分布を決める必要がある。
    • \lambda \sim \text{Exp}(\alpha)
    • τ \sim \text{DiscreteUniform}(1, 70)
  • \alphaについては、事前分布の決定はそれほど影響がないはずと考えて、
    • $\dfrac{1}{\alpha} = E[\lambda|\alpha] \approx \dfrac{1}{N}\displaystyle\sum C_{i}$によって決め打ちしてしまう。

ここまで考えたところをグラフィカルに表示すると、 graphical model

以下コーディング

モデルの記述

from __future__ import division
import pymc as pm
import numpy as np

alpha = 1 / count_data.mean()
# 変数の名前とPyMCに渡す「名前」は分かりやすさのために同じにしておくのが原則
lambda_1 = pm.Exponential("lambda_1", alpha)
lambda_2 = pm.Exponential("lambda_2", alpha)
tau = pm.DiscreteUniform("tau", lower=0, upper=n_count_data)

@pm.deterministic
def lambda_(tau=tau, lambda_1=lambda_1, lambda_2=lambda_2):
    out = np.empty(n_count_data)
    out[:tau] = lambda_1
    out[tau:] = lambda_2
    return out

#以下のコードが通るかは不明
#lambda_ = np.array([lambda_1 if t < tau else lambda_2
#                     for t in range(n_count_data)])

observation = pm.Poisson("obs", lambda_, value=count_data, observed=True)

model = pm.Model(observation, lambda_1, lambda_2, tau)

MCMCする

mcmc = pm.MCMC(model)
mcmc.sample(40000, 10000, 1)

値を取り出す

lambda_1_samples = mcmc.trace("lambda_1")[:]
lambda_2_samples = mcmc.trace("lambda_2")[:]
tau_samples = mcmc.trace("tau")[:]

プロット

f:id:kuyata:20140818212328p:plain

  • 変数には親子関係が成立する。今回を例にすると、 \lambdaは \alphaの子供になる。
  • PyMCにはdeterministicとstochasticの2通りがある。deterministicなものは引数を全て与えてしまえば値が決定するのに対して、stochasticなものは値が決定できずランダムなものである。
  • stochasticな変数はobservedという引数をとる(デフォルトはFalse)。Trueにすると、値が今の値で固定される。

例2 A/Bテスト

A/Bテストは2つの手続きの効果の違いを統計的に評価する方法のひとつ。

# 以後pymcとnumpyはインポートしてあるものとする
# ベルヌーイ分布からのサンプリング
p_true_A, p_true_B = 0.05, 0.04
N_A, N_B = 1500, 750
observed_A = pm.rbernoulli(p_true_A, N_A)
observed_B = pm.rbernoulli(p_true_B, N_B)

# モデルの記述
p_A = pm.Uniform("p_A", 0, 1)
p_B = pm.Uniform("p_B", 0, 1)

@pm.deterministic
def delta(p_A=p_A, p_B=p_B):
    return p_A - p_B
obs_A = pm.Bernoulli("obs_A", p_A, value=observed_A, observed=True)
obs_B = pm.Bernoulli("obs_B", p_B, value=observed_B, observed=True)

# mcmc
mcmc = pm.MCMC([p_A, p_B, delta, obs_A, obs_B])
mcmc.sample(20000, 1000)

# 取り出し
p_A_sample = pm.trace("p_A")[:]
p_B_sample = pm.trace("p_B")[:]
delta_sample = pm.trace("delta")[:]

例3 テストでカンニングした生徒の割合

カンニングしたか正直に答えてもらうのは難しいので、コインを使ってアンケート調査をする。コインを振って表なら正直に答える。裏ならもう一回コインを振ってその目が表なら「した」裏なら「していない」と答える。こうすればサンプルの半分は捨ててしまうことになるがバイアスがかかっていない答えを得ることができる。

# サンプル作成
N = 100
p = pm.Uniform("freq_cheating", 0, 1)
true_ans = pm.Bernoulli("true_ans", p, size=N)
first_coin = pm.Bernoulli("first_coin", 0.5, size=N)
second_coin = pm.Bernoulli("second_coin", 0.5, size=N)

@pm.deterministic
def obs_proportion(t_a=true_ans, f=first_flips, s=second_flips):
    obs = f * t_a + (1 - f) * s
    return obs.sum() / float(N)

X = 35
observation = pm.Binomial("observation", N, obs_proportion, value=X, observed=True)

model = pm.Model([p, true_ans, first_coin, second_coin, obs_proportion, observation])
mcmc = pm.MCMC([model])
mcmc.sample(40000, 15000)

# これでも似たような答えが出る(エスパーっぽい)
@pm.derterministic
def p_skewed(p=p):
    return 0.5 * p + 0.25

observation = pm.Binomial("number_cheaters", 100, p_skewed, value=X, observed=True)
mcmc = ...

例4 チャレンジャー号の事故について

O-ringという部品の不具合で打ち上げに失敗した。それ以前の24フライト中23でデータを得られていてそのうち7回でO-ringに障害があった。

f:id:kuyata:20140819164241p:plain

見たところ特に不具合が生じるか生じないかの境界はないようである。そこでロジスティック関数 p(\text{temp}) = \dfrac{1}{1 + e^{\beta \text{temp} + \alpha}}を用いてモデリングすることにする。

\alpha、 \betaの事前分布はそれぞれ特に制限する範囲もないので正規分布にする。 N(\mu, τ) = N(0, 0.001)と広めの分布をとっておく。

# data, tempはロードされているものとする
# 事前分布
beta = pm.Nornal("beta", 0, 0.001, value=0)
alpha = pm.Normal("alpha", 0, 0.001, value=0)
@pm.deterministic
def p(a=alpha, b=beta, t=temp):
    return 1 / (1 + np.exp(a + b * t))

observed = pm.Bernoulli("observed", p, value=data, observed=True)
model = pm.Model([observed, alpha, beta])
# 次の章で説明されるよくわからないコード
map_ = pm.MAP(model)
map_.fit()
mcmc = pm.MCMC(model)
mcmc.sample(120000, 100000, 2)

95%CIでプロット

import scipy.stats.mstats import mquantiles

qs = mquantiles(p_t, [0.025, 0.975], axis=0)
plt.fill_between(t[:, 0], *qs, alpha=0.7, color="blue")
...

f:id:kuyata:20140819185129p:plain

この章はこの後モデリングの正当性について触れているが視覚的に評価するだけで数学的な評価は避けていたのでパス。

  • 補足
  • MCMCするためのモデルからサンプリングすることが可能。
simulated = pm.Bernoulli("simulated", p)
N = 10000

mcmc = pm.MCMC([...])
mcmc.sample(N)
  • matplotlibにseparation_plotなるものがあるらしい