Welcome to Bayesian Methods for Hackers. The full Github repository is available at github/Probabilistic-Programming-and-Bayesian-Methods-for-Hackers. The other chapters can be found on the project's homepage. We hope you enjoy the book, and we encourage any contributions!
あなたは優秀なプログラマーだ.しかし,だれしも書いたコードにバグはある.実装するのが非常に難しいアルゴリズムのコードをなんとか書いた後,そのコードが正しいかどうかを簡単な例題でテストしようと考えた.OK,テストにパスした.次にもっと難しい問題でコードをテストした.今度もパスした.そして,なんともっともっと難しい問題でも,パスした! だからあなたは,このコードにはバグはないかもしれないと思い始めてしまった...
もしこのように考えたことがあるのなら,おめでとう! あなたはベイズ的に考る人の仲間入りだ.ベイズ推論とは, 新しい証拠が得られるたびに自分の考えを改めるというものである. ベイズ的に考える人は,ある結果が必ず起こるとは考えない.たぶん起こるだろうと考えるのである. 上の例のように,普通は,プログラムに100%まったくバグがないとは考えない. ないと言い切るには,実際にはありえないような場合も含めて,すべての場合についてチェックしなければならないだろう. それよりも,たくさんの問題にたいしてテストして,それら全てにパスしたら, プログラムにバグは「たぶんないだろう」と思うのである.しかし「まったくない」とは言い切れない. ベイズ推論も同じである.情報が得られたら信念を更新する.すべての可能な場合をチェックしなければ,絶対に,とは言わないのである.
ベイス推論が伝統的な統計的推論と異なるのは, 「不確実」なものは不確実なままにするという点である. 不確実なままということに,最初はダメな方法だと思うかもしれない. 統計とはランダムな現象から確実さを引き出すものではなかったのか? これを理解するには,ベイズ的に考えることが必要になる.
ベイズ的な考え方,つまりベイズ主義(Bayesian)では, 確率を「ある出来事がどのくらい信頼できるか」を表す指標と解釈する. つまり,ある事象が生じるということを, どのくらい確かだと思っているのか表すものと考える. すぐあとで見るように,実際にこれが確率を解釈する自然な方法である.
この解釈をもっと分かりやすくするために, 確率のもう一つの解釈を考えてみよう. それは「頻度主義」(Frequentist)というものであり, もっと古典的な統計学である. 頻度主義では,確率を「長期間における事象の頻度」とみなす (だから「頻度主義」という名前で呼ばれている). たとえば,飛行機事故の確率を頻度主義で考えれば, 「長期間における飛行機事故の頻度」になる. この考え方は,多くの場合,事象の確率として意味がある. しかし長期間にわたって事象が発生しないような場合には, 理解することが難しくなる. 例えば大統領選挙の結果の確率を計算しようとしても, ある特定の選挙は1回きりしか行われないのだ! 頻度主義でこの問題を避けるには, 他のすべての選挙も考慮して, これらの発生する頻度で確率を定義することになる.
一方のベイズ主義では,もっと直感的に考える. ベイズ主義では,確率を, ある事象が発生する信念(belief) もしくは確信(confidence)の度合いとみなす. 確率とは,思っていることを要約したものであるだけなのである. ある人が,ある事象の信念を0だと思っている場合, その事象が発生するとは考えていないことになる. 反対に,ある事象の信念を1だと思っている場合, その事象が必ず発生すると考えていることになる. 信念を0から1の実数値で表せば,それを使って 他の結果に重みを付けることができる. この定義を使えば,飛行機事故の確率の例を うまく表現することができる. 飛行機事故の頻度が得られた時に,他の情報が何もなければ, ある人の信念はその頻度に一致するべきである. 同様に,確率が信念であるという定義を使えば, 大統領選挙の結果の確率(信念)というものを考えてもよいだろう. つまり,ある候補者Aが当選することをどのくらい確信しているのか, を表しているのである.
上のパラグラフでは,一般的な信念(確率)ではなく, 「ある人」の信念(確率)というものを説明していたことに注意して欲しい. これが面白いのは,人によって信念が違うということを定義が許していることになる点である. これは日常では見かける.この世界について持っている情報は人それぞれ違うので, ある人の信念は別の人の信念とは違うのである. 信念が違っているということは,「誰かが間違っている」ということではない. 以下の例で,各個人が持っている信念と確率との関係を考えてみてみよう.
人間にとっては,確率を信念のように扱うということは自然なやり方である.この世の中で生きていくために,いつもこのようなやり方をしているし,真実というものが完全ではないという例もたくさん見ている.それでも信念から何か情報がえられないかと考えてもいる.もし頻度主義のように考えるとするなら,かなり訓練しなければならないだろう.
従来の確率論の記法に従って, ある事象$A$が生じるという信念を$P(A)$と表し, 事前確率(prior probability)と呼ぶことにする.
偉大な経済学者であり思想家であるジョン・メイナード・ケインズ曰く, 「事実が変わったならば,わたしは考えを改める.あなたはどうしますか?」 これは証拠が得られた後に信念を更新す,ベイズ主義的なやり方である. もし,その証拠が最初に思っていた信念と相反することであったとしても, 証拠を無視することはできない. この更新された信念を$P(A |X )$と表し, 証拠$X$が与えられた時の$A$の確率である,と解釈する. 事前確率に対応して,この更新された信念を事後確率(posterior probability)と呼ぶ. 例えば上記の例では,証拠$X$が得られた後の事後確率(事後信念と言ってもよい)は次のようになる.
1. $P(A): \;\;$ コインの表が出る確率は50パーセントである.$P(A | X):\;\;$ コインの結果をみて表が出ていたとする.この情報を$X$とする.明らかに,表の事後確率は1.0で,裏の事後確率は0.0である.
2. $P(A): \;\;$ この複雑で巨大なプログラムにはたぶんバグがある.$P(A | X): \;\;$ プログラムはすべてのテスト$X$にパスした.たぶんバグがあるかもしれないが,その可能性は非常に小さいだろう.
3. $P(A):\;\;$ ある患者が病気にかかっている.$P(A | X):\;\;$ 血液検査の結果,$X$という証拠が得られたので,いくつかの病気の可能性は排除してもよいだろう.
これらの例から明らかなように,新しい証拠$X$が得られたとしても, 事前の信念を完全に否定することはなく,新しい証拠を事前確率の重みとして使っている (つまり,ある信念にはより大きい重み,つまり確信度を与えるのである).
事前にある事象がどのくらい生じるのかということを考えても,それは非常に不確実である. したがって,どんな結果を予想したとしても,間違っている可能性が高い. データや証拠や情報が得られれば,信念を更新して, 間違っている可能性がもっと少ない予測をすることができるようになる. コインの裏表を予測すること例では,正しい予測ができるようになる.
頻度主義とベイス主義がプログラミング言語の関数だったら,統計的な問題を入力すると,ユーザーに返される結果は同じではないだろう. 頻度主義の推論関数の戻り値は,推定値を表す数値である(標本平均などの要約統計量であることが多い). 一方でベイズ主義の推論関数は,「確率」を返す.
例えばデバッグの例題であれば,頻度主義関数の引数に 「このプログラムはテスト$X$のすべてをパスしたんだ.このプログラムにはバグがないかな?」を渡すと, 戻り値は「バグはありません」だろう. しかしベイズ主義関数に 「プログラムを書くといつもバグがあるんだ. このプログラムはテスト$X$のすべてをパスしたんだ.このプログラムにはバグがないかな?」という引数を渡すと, 「バグはありません」と「バグがあります」の答えのそれぞれに確率が返される.
バグがない確率は0.8,バグがある確率は0.2です.
この戻り値は頻度主義関数の戻り値とはまったく違うものである. ベイズ主義関数は引数に「プログラムを書くといつもバグがあるんだ」という情報を追加していることに気がついてほしい. これが事前情報(prior)である. この事前情報パラメータを引数に与えることで,今の状況についての信念をベイズ主義関数に伝えている. これを与えるかどうかはユーザーの自由だが,与えない場合には別の結果が得られることになる. その例は後で見ることにしよう.
証拠をたくさん手に入れることができれば,事前の信念は,その多数の証拠にかき消されてしまう. これは想像できるだろう. 例えば,あなたが「今日,太陽が爆発するんじゃないか」という事前信念を持っていたとすれば,日に日にその信念は揺らいでいき, そして,どんな推論でもいいから自分の間違いを正してくれ,少なくともこの信念をもっとマシなものにしてくれ, と思うようになる(かもしれない). そしてベイズ推論は,その信念を正してくれる.
$N$を手に入る証拠の数とする.もし無限個の証拠が手に入れば,つまり$N \rightarrow \infty$ならば, ベイズ推論の結果は頻度主義の結果と(多くの場合)一致する. したがって$N$が大きくなれば,統計的推論は客観的なものになる. 反対に$N$が小さければ,推論は不安定なものになる. 頻度主義の推定値は分散も信頼区間も大きくなる. そんな時にはベイズ推論の出番である. 事前分布を引数にとり,結果に(推定値ではなく)確率を出力する. これは,$N$の小さいデータセットに対する統計的推論の不安定さを反映した,不確実さを表すものになっている.
$N$が非常に大きい場合は,頻度主義とベイズ主義は似たような推論結果を出してくるので,二つの区別はつかなくなるだろう. そのため,少ない計算で済む頻度主義を用いたくなるかもしれない. もしそんな状況にあるのであれば,そうする前に以下の [Andrew Gelman (2005)][1]の文章を読んでほしい.
サンプル数が大きい場合,というものは存在しない.もし$N$が小さすぎて十分に正確な推定値を得ることができないのであれば,データをもっと増やす(もしくはもっと多くの仮定を使う)必要がある.しかし,もし$N$が「十分に大きい」のであれば,データを分割してもっと多くの情報を得ることができるだろう(例えば世論調査の場合には,全国区での良い推定値が得られたら,次は男女別,地域別,年齢別の推定値を得ることもできるだろう).$N$が十分であることはない.もし「十分」だとしたら,あなたはもうすでにもっと多くのデータを必要とする次の問題に取り組んでいるのだ.
間違ってはいない.
頻度主義の方法は今でも多くの分野で有用であり,最先端で使われれいる. 最小二乗回帰やlasso回帰,EMアルゴリズムなどのツールはどれも優れていて処理も速い. ベイズ主義の手法は,それらの手法を補うものである. それらの手法が適用できない問題を解いたり, もっと柔軟なモデル化で隠れた構造を解き明かしたりするのである.
逆説的に聞こえるかもしれないが,ビッグデータで予測したり解析したりする問題には, 実際には比較的単純なアルゴリズムが使われている[2][4]. ビッグデータを用いた予測の難しさは,アルゴリズムにあるのではない. ビッグデータを保存し読み出すストレージや ビッグデータに対して実行する時の計算量が大変なのである. (上述のGelmanの文章を読んで「自分は本当にビッグデータを持っているのだろうか?」と考えてみてほしい)
解析するのがもっと難しい問題は,「ミディアムなデータ」の場合であり, 特に問題となるのは「スモールデータ」の場合である. Gelmanの文章を借りるなら,ビッグデータの問題が「十分にビッグ」で実際には解けないのであれば, 「それほど十分にビッグではない」データを扱えばよいのである.
計算するべき信念は,ベイズ的に考えた確率と解釈することができる. ここで,ある事象$A$について「事前」信念を持っているとしよう (例えば,テストを実行する前に,プログラムにバグがありそうかどうかについての信念).
次は,得られた証拠を使おう.バグありプログラムの例を使えば, プログラムはテスト$X$にパスしたので,その情報を取り入れて信念を更新したい. この更新された新しい信念を「事後」信念と呼ぶことにする. 以下の式を使えば,信念を更新することができる. この式は,発見者のトーマス・ベイズにちなんで,ベイズの定理と呼ばれている.
\begin{align} P( A | X ) = & \frac{ P(X | A) P(A) } {P(X) } \\\\[5pt] & \propto P(X | A) P(A)\;\; (\propto \text{は比例を表す} ) \end{align}この公式はベイズ推論だけのものではない.ベイズ推論以外でも使われている数学的事実である. ベイズ推論では,単にこの式を使って 初期の事前確率$P(A)$と更新後の事後確率$P(A | X)$を結びつけているだけである.
統計学のテキストであれば,コイン投げの問題を扱っていない本はない. ちょっと変わったやり方でこの問題を扱ってみよう. あなたは,コインの表が出る確率がよくわからないとする(本当は50%). 何らかの比率(ここでは$p$とする)で表裏がでるということについては信じているが, その$p$がどのくらいなのかについては,まったく情報を持っていない.
ではコイン投げを初めて,表$H$が出たのか裏$T$が出たのかを記録することにする. ここでちょっと考えてみよう. データが増えるにつれて,推論結果はどのように変わっていくのだろうか? もっと正確に言えば,データが少ない時とデータが多い時とで,事後確率はどのように違うのだろうか?
以下のコードは, (コイン投げの)データが増えるたびに更新される事後確率の系列をプロットするものである.
"""
本書ではmatplotlibのグラフのスタイルを変更するために,matplotlibrcファイルをカスタマイズしている.
本書を実行して,本書のスタイルを使いたいのであれば,以下の2つの方法がある.
1. 本書の style/ ディレクトリにあるrcフィアルで,自分の環境のmatplotlibrcを書き換える.
http://matplotlib.org/users/customizing.htmlを参照.
2. スタイルはbmh_matplotlibrc.jsonファイルにもある.これを使って以下のコードを実行すれば,
本書にだけスタイルを適用することができる.
import json
s = json.load( open("../styles/bmh_matplotlibrc.json") )
matplotlib.rcParams.update(s)
"""
# 以下のコードは読み飛ばして構わない.ここではあまり重要ではないし,
# まだ説明していない進んだ内容も含んでいる.その下のグラフを見て!
%matplotlib inline
from IPython.core.pylabtools import figsize
import numpy as np
from matplotlib import pyplot as plt
figsize(11, 9)
import scipy.stats as stats
dist = stats.beta
n_trials = [0, 1, 2, 3, 4, 5, 8, 15, 50, 500]
data = stats.bernoulli.rvs(0.5, size=n_trials[-1])
x = np.linspace(0, 1, 100)
# 分かっている人へ:ここでは二項分布の共役事前分布を使っている.
for k, N in enumerate(n_trials):
sx = plt.subplot(len(n_trials) / 2, 2, k + 1)
# u"$p$, 表が出る確率"
plt.xlabel("$p$, probability of heads") \
if k in [0, len(n_trials) - 1] else None
plt.setp(sx.get_yticklabels(), visible=False)
heads = data[:N].sum()
y = dist.pdf(x, 1 + heads, 1 + N - heads)
plt.plot(x, y, label="observe %d tosses,\n %d heads" % (N, heads)) # u"%d回投げて, \n 表は%d回"
plt.fill_between(x, 0, y, color="#348ABD", alpha=0.4)
plt.vlines(0.5, 0, 4, color="k", linestyles="--", lw=1)
leg = plt.legend()
leg.get_frame().set_alpha(0.4)
plt.autoscale(tight=True)
plt.suptitle("Bayesian updating of posterior probabilities", # u"ベイズ推論による事後確率の更新"
y=1.02,
fontsize=14)
plt.tight_layout()
事後確率は曲線で表されている. 不確実さは,この曲線の広がり具合に比例している. 上のグラフを見れば分かるように, データが増えるたびに事後確率の曲線は右へ左へと動き回る. 最終的に,データがたくさん手に入れば(たくさんコインを投げたら), 事後確率曲線は,真の確率である$p=0.5$に次第に集まってくる.
なお,曲線のピークの位置は0.5ではないし,そうである理由もない.$p$の値については何も知らない,という前提だったのだから. 実際のとこr,コイン投げの結果が極端な,たとえば8回投げて表が1回しかなかったような場合には, 事後確率曲線のピークは0.5から非常に離れているだろう (事前情報がないのだから,8回投げて表が1回しか出ないコインにイカサマはないと,どのくらい確信できるだろう?). もっとデータが増えれば,確率はもっと$p=0.5$に近くなるだろう.
次の例で,数学がベイズ推論でどのように使われるのか見てみよう.
「プログラムにはバグがない」という事象を$A$とする. 「このプログラムがすべてのデバッグテストにパスする」という事象を$X$とする. とりあえず,バグがないという事前確率$P(A)$を変数$p$にしておこう. つまり$P(A) = p$とする.
今から考えるのはこの$P(A|X)$だ. つまり,「デバッグテスト$X$をパスした時に,バグがない」確率である. 上の公式を使うために,いくつか計算しなければならない.
では$P(X | A)$とは何だろう? これは,「バグがない時にすべてのテスト$X$をパスする」確率である.明らかに,これは1である.バグがなければ,どんなテストにもパスするからだ.
それよりも厄介なのは$P(X)$である.事象$X$が起きる可能性は2つある. 実はバグがある(これを$\sim A$や$\lnot A$と書いてnot $A$と読む)にも関わらず事象$X$が起こっているのか,それともバグがないから事象$X$が起きているのか,である. すると,$P(X)$は以下のように解釈できる.
すでに$P(X|A)$は計算してある.しかし$P(X | \sim A)$をどうするかは,主観的である. プログラムはテストをパスしたが,それでもバグがあるのだ. しかしバグがある確率は小さくなっている. これは実行したテストの数や,テストがどれだけ精巧なのかにも依存する. ここでは控えめに考えて,$P(X|\sim A) = 0.5$としよう.すると以下のようになる.
\begin{align} P(A | X) & = \frac{1\cdot p}{ 1\cdot p +0.5 (1-p) } \\\\ & = \frac{ 2 p}{1+p} \end{align}これが事後確率である.これを事前確率パラメータ$p \in [0,1]$の関数としてみたら, どんな形をしているだろう?
figsize(12.5, 4) # グラフの縦横サイズを12.5:4にする
p = np.linspace(0, 1, 50) # 0から1までを50点に分割
plt.plot(p, 2 * p / (1 + p), color="#348ABD", lw=3) # 事後確率をプロット.色は青系,線幅は3
# plt.fill_between(p, 2*p/(1+p), alpha=.5, facecolor=["#A60628"]) # コメントを外して試してみよう
plt.scatter(0.2, 2 * (0.2) / 1.2, s=140, c="#348ABD") # p=0.2のところに点を描画.色は青系,サイズは140
plt.xlim(0, 1) # x軸の範囲を(0,1)に設定
plt.ylim(0, 1) # y軸の範囲を(0,1)に設定
plt.xlabel("Prior, $P(A) = p$") # x軸ラベル u"事前確率$P(A) = p$"
plt.ylabel("Posterior, $P(A|X)$, with $P(A) = p$") # y軸ラベル u"$P(A) = p$の時の事後確率$P(A|X)$"
plt.title("Are there bugs in my code?") # グラフのタイトル u"プログラムにバグがあるか?"
<matplotlib.text.Text at 0x10850ee90>
事前確率$p$が小さい時は,テスト$X$にパスしたという証拠が非常に効いている. ここで事前確率の値を一つ決めてみよう. 私は優秀なプログラマーなので(自分ではそう思っている),0.20でも現実的だろう. つまり,20%の確率でバグのないプログラムを書くことができる,というわけだ. もっとも現実的には,この事前確率は プログラムがどれだけ複雑で大規模なのかにもよるのだが, とりあえず0.20としておこう. すると,プログラムにはバグがないという更新された信念は0.33となる.
ここで事前確率は確率である,ということを思い出しておこう. $p$はバグがない事前確率で,$1-p$はバグがある事前確率である
同様に,事後確率も確率である. $P(A | X)$は「すべてのテストをパスしてバグがない」確率, $1-P(A|X)$は「すべてのテストをパスしてバグがある」確率である. この事後確率はどんな値だろうか? 以下のグラフは,事前確率と事後確率を計算したものである.
figsize(12.5, 4)
colours = ["#348ABD", "#A60628"]
prior = [0.20, 0.80]
posterior = [1. / 3, 2. / 3]
plt.bar([0, .7], prior, alpha=0.70, width=0.25,
color=colours[0], label="prior distribution", # u"事前確率"
lw="3", edgecolor=colours[0])
plt.bar([0 + 0.25, .7 + 0.25], posterior, alpha=0.7,
width=0.25, color=colours[1],
label="posterior distribution", # u"事後確率"
lw="3", edgecolor=colours[1])
plt.xticks([0.20, .95], ["Bugs Absent", "Bugs Present"]) # u"バグがない", u"バグがある
plt.title("Prior and Posterior probability of bugs present") # u"バグがある事前確率と事後確率"
plt.ylabel("Probability") # u"確率"
plt.legend(loc="upper left"); # 凡例の位置は左上
<matplotlib.legend.Legend at 0x1085247d0>
証拠となる事象$X$が得られた後では,バグがないという確率が大きくなっている. テストの数を増やせば,バグがない(確率1)ということを確信できるだろう.
これはベイズ推論とベイズ則の非常に単純な例である. 残念ながら,もう少し複雑なベイズ推論を実行するための数学は, 非常に調整された例題でなければ,もっともっと難しくなってしまう. あとで見るように,この手の数学的な解析は実際には必要ない. いろいろなモデル化のためのツールを知るほうが先である. 次の節は「確率分布」を扱う.もしよく知っているなら, 読み飛ばして(もしくは斜め読みして)もよい. よく知らなければ,非常に重要なのでよく理解してほしい.
確率分布とは何かを簡単におさらいしよう. $Z$を確率変数とする. $Z$が取る値それぞれに確率を与えるのが確率分布である. グラフで描けば,確率分布は曲線で書くことができて, その曲線の高さに比例して確率が大きくなる. すでにこの章の最初の図で,確率分布の曲線の例を説明している.
確率変数には以下のように3種類ある.
もし$Z$が離散なら,確率分布は確率質量関数と呼ばれる. これは$Z$が値$k$を取る確率を$P(Z=k)$で表す. 確率質量関数は確率変数$Z$を完全に決定する.つまり 確率質量関数が分かれば$Z$がどのように振る舞うのかが分かるのである. この先よく登場する有名な確率質量関数がいくつかあるが, 必要に応じて紹介する.一番最初に紹介する有用な確率質量関数は,ポワソン分布である. $Z$の確率質量関数が以下の式の時,$Z$はポワソン分布に従うと言う.
$$P(Z = k) =\frac{ \lambda^k e^{-\lambda} }{k!}, \; \; k=0,1,2, \dots $$$\lambda$は分布のパラメータで,分布の形状を決める. ポワソン分布の場合,$\lambda$は任意の正の実数である. $\lambda$を大きくすると大きな値の確率が高くなり, $\lambda$を小さくすると小さな値の確率が高くなる. そのため,$\lambda$はポワソン分布の「強度」と考えてもよい.
任意の正の実数である$\lambda$とは異なり,上の式での$k$の値は正の整数,つまり0, 1, 2, ... でなければならない.これは重要な事である.人数をモデル化しようとするなら,4.25人とか5.612人というのは意味が無いのだから.
もし確率変数$Z$がポワソン分布に従うなら,それを以下のように書く.
$$Z \sim \text{Poi}(\lambda) $$ポワソン分布の便利な性質の一つは,期待値が分布パラメータに等しいということである.
$$E[ \;Z\; | \; \lambda \;] = \lambda $$この先この性質を利用するので覚えておいてほしい. 以下のグラフは,いくつかの$\lambda$の値について確率質量関数をプロットしたものである. このグラフで注意してほしいことは2つ.1つ目は$\lambda$を大きくすれば, 大きな値の確率が高くなること.2つ目は,グラフの横軸は15で終わっているが, 分布はそうではない.すべての正の整数に対して確率が割り当てられている.
figsize(12.5, 4)
import scipy.stats as stats
a = np.arange(16)
poi = stats.poisson
lambda_ = [1.5, 4.25]
colours = ["#348ABD", "#A60628"]
plt.bar(a, poi.pmf(a, lambda_[0]), color=colours[0],
label="$\lambda = %.1f$" % lambda_[0], alpha=0.60,
edgecolor=colours[0], lw="3")
plt.bar(a, poi.pmf(a, lambda_[1]), color=colours[1],
label="$\lambda = %.1f$" % lambda_[1], alpha=0.60,
edgecolor=colours[1], lw="3")
plt.xticks(a + 0.4, a)
plt.legend()
plt.ylabel("probability of $k$") # u"$k$の確率"
plt.xlabel("$k$")
plt.title("Probability mass function of a Poisson random variable; differing \
$\lambda$ values") # u"いくつかの$\lambda$に対するポワソン分布の確率質量関数"
<matplotlib.text.Text at 0x108588a90>
連続確率変数は,確率質量関数ではなく確率密度分布関数で表される. これは単なる名前の問題のように見えるかもしれないが, 密度関数と質量関数はまったく違うものなのである. 連続確率変数の例は,以下の式で表される指数分布である.
$$f_Z(z | \lambda) = \lambda e^{-\lambda z }, \;\; z\ge 0$$Like a Poisson random variable, an exponential random variable can take on only non-negative values. But unlike a Poisson variable, the exponential can take on any non-negative values, including non-integral values such as 4.25 or 5.612401. This property makes it a poor choice for count data, which must be an integer, but a great choice for time data, temperature data (measured in Kelvins, of course), or any other precise and positive variable. The graph below shows two probability density functions with different $\lambda$ values.
When a random variable $Z$ has an exponential distribution with parameter $\lambda$, we say $Z$ is exponential and write
$$Z \sim \text{Exp}(\lambda)$$Given a specific $\lambda$, the expected value of an exponential random variable is equal to the inverse of $\lambda$, that is:
$$E[\; Z \;|\; \lambda \;] = \frac{1}{\lambda}$$a = np.linspace(0, 4, 100)
expo = stats.expon
lambda_ = [0.5, 1]
for l, c in zip(lambda_, colours):
plt.plot(a, expo.pdf(a, scale=1. / l), lw=3,
color=c, label="$\lambda = %.1f$" % l)
plt.fill_between(a, expo.pdf(a, scale=1. / l), color=c, alpha=.33)
plt.legend()
plt.ylabel("PDF at $z$")
plt.xlabel("$z$")
plt.ylim(0, 1.2)
plt.title("Probability density function of an Exponential random variable;\
differing $\lambda$");
This question is what motivates statistics. In the real world, $\lambda$ is hidden from us. We see only $Z$, and must go backwards to try and determine $\lambda$. The problem is difficult because there is no one-to-one mapping from $Z$ to $\lambda$. Many different methods have been created to solve the problem of estimating $\lambda$, but since $\lambda$ is never actually observed, no one can say for certain which method is best!
Bayesian inference is concerned with beliefs about what $\lambda$ might be. Rather than try to guess $\lambda$ exactly, we can only talk about what $\lambda$ is likely to be by assigning a probability distribution to $\lambda$.
This might seem odd at first. After all, $\lambda$ is fixed; it is not (necessarily) random! How can we assign probabilities to values of a non-random variable? Ah, we have fallen for our old, frequentist way of thinking. Recall that under Bayesian philosophy, we can assign probabilities if we interpret them as beliefs. And it is entirely acceptable to have beliefs about the parameter $\lambda$.
Let's try to model a more interesting example, one that concerns the rate at which a user sends and receives text messages:
You are given a series of daily text-message counts from a user of your system. The data, plotted over time, appears in the chart below. You are curious to know if the user's text-messaging habits have changed over time, either gradually or suddenly. How can you model this? (This is in fact my own text-message data. Judge my popularity as you wish.)
figsize(12.5, 3.5)
count_data = np.loadtxt("data/txtdata.csv")
n_count_data = len(count_data)
plt.bar(np.arange(n_count_data), count_data, color="#348ABD")
plt.xlabel("Time (days)")
plt.ylabel("count of text-msgs received")
plt.title("Did the user's texting habits change over time?")
plt.xlim(0, n_count_data);
Before we start modeling, see what you can figure out just by looking at the chart above. Would you say there was a change in behaviour during this time period?
How can we start to model this? Well, as we have conveniently already seen, a Poisson random variable is a very appropriate model for this type of count data. Denoting day $i$'s text-message count by $C_i$,
$$ C_i \sim \text{Poisson}(\lambda) $$We are not sure what the value of the $\lambda$ parameter really is, however. Looking at the chart above, it appears that the rate might become higher late in the observation period, which is equivalent to saying that $\lambda$ increases at some point during the observations. (Recall that a higher value of $\lambda$ assigns more probability to larger outcomes. That is, there is a higher probability of many text messages having been sent on a given day.)
How can we represent this observation mathematically? Let's assume that on some day during the observation period (call it $\tau$), the parameter $\lambda$ suddenly jumps to a higher value. So we really have two $\lambda$ parameters: one for the period before $\tau$, and one for the rest of the observation period. In the literature, a sudden transition like this would be called a switchpoint:
$$ \lambda = \begin{cases} \lambda_1 & \text{if } t \lt \tau \cr \lambda_2 & \text{if } t \ge \tau \end{cases} $$If, in reality, no sudden change occurred and indeed $\lambda_1 = \lambda_2$, then the $\lambda$s posterior distributions should look about equal.
We are interested in inferring the unknown $\lambda$s. To use Bayesian inference, we need to assign prior probabilities to the different possible values of $\lambda$. What would be good prior probability distributions for $\lambda_1$ and $\lambda_2$? Recall that $\lambda$ can be any positive number. As we saw earlier, the exponential distribution provides a continuous density function for positive numbers, so it might be a good choice for modeling $\lambda_i$. But recall that the exponential distribution takes a parameter of its own, so we'll need to include that parameter in our model. Let's call that parameter $\alpha$.
\begin{align} &\lambda_1 \sim \text{Exp}( \alpha ) \\\ &\lambda_2 \sim \text{Exp}( \alpha ) \end{align}$\alpha$ is called a hyper-parameter or parent variable. In literal terms, it is a parameter that influences other parameters. Our initial guess at $\alpha$ does not influence the model too strongly, so we have some flexibility in our choice. A good rule of thumb is to set the exponential parameter equal to the inverse of the average of the count data. Since we're modeling $\lambda$ using an exponential distribution, we can use the expected value identity shown earlier to get:
$$\frac{1}{N}\sum_{i=0}^N \;C_i \approx E[\; \lambda \; |\; \alpha ] = \frac{1}{\alpha}$$An alternative, and something I encourage the reader to try, would be to have two priors: one for each $\lambda_i$. Creating two exponential distributions with different $\alpha$ values reflects our prior belief that the rate changed at some point during the observations.
What about $\tau$? Because of the noisiness of the data, it's difficult to pick out a priori when $\tau$ might have occurred. Instead, we can assign a uniform prior belief to every possible day. This is equivalent to saying
\begin{align} & \tau \sim \text{DiscreteUniform(1,70) }\\\\ & \Rightarrow P( \tau = k ) = \frac{1}{70} \end{align}So after all this, what does our overall prior distribution for the unknown variables look like? Frankly, it doesn't matter. What we should understand is that it's an ugly, complicated mess involving symbols only a mathematician could love. And things will only get uglier the more complicated our models become. Regardless, all we really care about is the posterior distribution.
We next turn to PyMC, a Python library for performing Bayesian analysis that is undaunted by the mathematical monster we have created.
PyMC is a Python library for programming Bayesian analysis [3]. It is a fast, well-maintained library. The only unfortunate part is that its documentation is lacking in certain areas, especially those that bridge the gap between beginner and hacker. One of this book's main goals is to solve that problem, and also to demonstrate why PyMC is so cool.
We will model the problem above using PyMC. This type of programming is called probabilistic programming, an unfortunate misnomer that invokes ideas of randomly-generated code and has likely confused and frightened users away from this field. The code is not random; it is probabilistic in the sense that we create probability models using programming variables as the model's components. Model components are first-class primitives within the PyMC framework.
B. Cronin [5] has a very motivating description of probabilistic programming:
Another way of thinking about this: unlike a traditional program, which only runs in the forward directions, a probabilistic program is run in both the forward and backward direction. It runs forward to compute the consequences of the assumptions it contains about the world (i.e., the model space it represents), but it also runs backward from the data to constrain the possible explanations. In practice, many probabilistic programming systems will cleverly interleave these forward and backward operations to efficiently home in on the best explanations.
Because of the confusion engendered by the term probabilistic programming, I'll refrain from using it. Instead, I'll simply say programming, since that's what it really is.
PyMC code is easy to read. The only novel thing should be the syntax, and I will interrupt the code to explain individual sections. Simply remember that we are representing the model's components ($\tau, \lambda_1, \lambda_2$ ) as variables:
import pymc as pm
alpha = 1.0 / count_data.mean() # Recall count_data is the
# variable that holds our txt counts
with pm.Model() as model:
lambda_1 = pm.Exponential("lambda_1", alpha)
lambda_2 = pm.Exponential("lambda_2", alpha)
with model:
tau = pm.DiscreteUniform("tau", lower=0, upper=n_count_data)
In the code above, we create the PyMC variables corresponding to $\lambda_1$ and $\lambda_2$. We assign them to PyMC's stochastic variables, so-called because they are treated by the back end as random number generators. We can demonstrate this fact by calling their built-in random()
methods.
print "Random output:", tau.random(), tau.random(), tau.random()
Random output:
--------------------------------------------------------------------------- AttributeError Traceback (most recent call last) <ipython-input-14-04009b63e469> in <module>() ----> 1 print "Random output:", tau.random(), tau.random(), tau.random() AttributeError: 'FreeRV' object has no attribute 'random'
@pm.deterministic
def lambda_(tau=tau, lambda_1=lambda_1, lambda_2=lambda_2):
out = np.zeros(n_count_data)
out[:tau] = lambda_1 # lambda before tau is lambda1
out[tau:] = lambda_2 # lambda after (and including) tau is lambda2
return out
--------------------------------------------------------------------------- AttributeError Traceback (most recent call last) <ipython-input-13-55d04f3ec163> in <module>() ----> 1 @pm.deterministic 2 def lambda_(tau=tau, lambda_1=lambda_1, lambda_2=lambda_2): 3 out = np.zeros(n_count_data) 4 out[:tau] = lambda_1 # lambda before tau is lambda1 5 out[tau:] = lambda_2 # lambda after (and including) tau is lambda2 AttributeError: 'module' object has no attribute 'deterministic'
This code creates a new function lambda_
, but really we can think of it as a random variable: the random variable $\lambda$ from above. Note that because lambda_1
, lambda_2
and tau
are random, lambda_
will be random. We are not fixing any variables yet.
@pm.deterministic
is a decorator that tells PyMC this is a deterministic function. That is, if the arguments were deterministic (which they are not), the output would be deterministic as well.
observation = pm.Poisson("obs", lambda_, value=count_data, observed=True)
model = pm.Model([observation, lambda_1, lambda_2, tau])
The variable observation
combines our data, count_data
, with our proposed data-generation scheme, given by the variable lambda_
, through the value
keyword. We also set observed = True
to tell PyMC that this should stay fixed in our analysis. Finally, PyMC wants us to collect all the variables of interest and create a Model
instance out of them. This makes our life easier when we retrieve the results.
The code below will be explained in Chapter 3, but I show it here so you can see where our results come from. One can think of it as a learning step. The machinery being employed is called Markov Chain Monte Carlo (MCMC), which I also delay explaining until Chapter 3. This technique returns thousands of random variables from the posterior distributions of $\lambda_1, \lambda_2$ and $\tau$. We can plot a histogram of the random variables to see what the posterior distributions look like. Below, we collect the samples (called traces in the MCMC literature) into histograms.
# Mysterious code to be explained in Chapter 3.
mcmc = pm.MCMC(model)
mcmc.sample(40000, 10000, 1)
--------------------------------------------------------------------------- AttributeError Traceback (most recent call last) <ipython-input-15-05abc05ea0f7> in <module>() 1 # Mysterious code to be explained in Chapter 3. ----> 2 mcmc = pm.MCMC(model) 3 mcmc.sample(40000, 10000, 1) AttributeError: 'module' object has no attribute 'MCMC'
lambda_1_samples = mcmc.trace('lambda_1')[:]
lambda_2_samples = mcmc.trace('lambda_2')[:]
tau_samples = mcmc.trace('tau')[:]
figsize(12.5, 10)
# histogram of the samples:
ax = plt.subplot(311)
ax.set_autoscaley_on(False)
plt.hist(lambda_1_samples, histtype='stepfilled', bins=30, alpha=0.85,
label="posterior of $\lambda_1$", color="#A60628", normed=True)
plt.legend(loc="upper left")
plt.title(r"""Posterior distributions of the variables
$\lambda_1,\;\lambda_2,\;\tau$""")
plt.xlim([15, 30])
plt.xlabel("$\lambda_1$ value")
ax = plt.subplot(312)
ax.set_autoscaley_on(False)
plt.hist(lambda_2_samples, histtype='stepfilled', bins=30, alpha=0.85,
label="posterior of $\lambda_2$", color="#7A68A6", normed=True)
plt.legend(loc="upper left")
plt.xlim([15, 30])
plt.xlabel("$\lambda_2$ value")
plt.subplot(313)
w = 1.0 / tau_samples.shape[0] * np.ones_like(tau_samples)
plt.hist(tau_samples, bins=n_count_data, alpha=1,
label=r"posterior of $\tau$",
color="#467821", weights=w, rwidth=2.)
plt.xticks(np.arange(n_count_data))
plt.legend(loc="upper left")
plt.ylim([0, .75])
plt.xlim([35, len(count_data) - 20])
plt.xlabel(r"$\tau$ (in days)")
plt.ylabel("probability");
Recall that Bayesian methodology returns a distribution. Hence we now have distributions to describe the unknown $\lambda$s and $\tau$. What have we gained? Immediately, we can see the uncertainty in our estimates: the wider the distribution, the less certain our posterior belief should be. We can also see what the plausible values for the parameters are: $\lambda_1$ is around 18 and $\lambda_2$ is around 23. The posterior distributions of the two $\lambda$s are clearly distinct, indicating that it is indeed likely that there was a change in the user's text-message behaviour.
What other observations can you make? If you look at the original data again, do these results seem reasonable?
Notice also that the posterior distributions for the $\lambda$s do not look like exponential distributions, even though our priors for these variables were exponential. In fact, the posterior distributions are not really of any form that we recognize from the original model. But that's OK! This is one of the benefits of taking a computational point of view. If we had instead done this analysis using mathematical approaches, we would have been stuck with an analytically intractable (and messy) distribution. Our use of a computational approach makes us indifferent to mathematical tractability.
Our analysis also returned a distribution for $\tau$. Its posterior distribution looks a little different from the other two because it is a discrete random variable, so it doesn't assign probabilities to intervals. We can see that near day 45, there was a 50% chance that the user's behaviour changed. Had no change occurred, or had the change been gradual over time, the posterior distribution of $\tau$ would have been more spread out, reflecting that many days were plausible candidates for $\tau$. By contrast, in the actual results we see that only three or four days make any sense as potential transition points.
We will deal with this question for the remainder of the book, and it is an understatement to say that it will lead us to some amazing results. For now, let's end this chapter with one more example.
We'll use the posterior samples to answer the following question: what is the expected number of texts at day $t, \; 0 \le t \le 70$ ? Recall that the expected value of a Poisson variable is equal to its parameter $\lambda$. Therefore, the question is equivalent to what is the expected value of $\lambda$ at time $t$?
In the code below, let $i$ index samples from the posterior distributions. Given a day $t$, we average over all possible $\lambda_i$ for that day $t$, using $\lambda_i = \lambda_{1,i}$ if $t \lt \tau_i$ (that is, if the behaviour change has not yet occurred), else we use $\lambda_i = \lambda_{2,i}$.
figsize(12.5, 5)
# tau_samples, lambda_1_samples, lambda_2_samples contain
# N samples from the corresponding posterior distribution
N = tau_samples.shape[0]
expected_texts_per_day = np.zeros(n_count_data)
for day in range(0, n_count_data):
# ix is a bool index of all tau samples corresponding to
# the switchpoint occurring prior to value of 'day'
ix = day < tau_samples
# Each posterior sample corresponds to a value for tau.
# for each day, that value of tau indicates whether we're "before"
# (in the lambda1 "regime") or
# "after" (in the lambda2 "regime") the switchpoint.
# by taking the posterior sample of lambda1/2 accordingly, we can average
# over all samples to get an expected value for lambda on that day.
# As explained, the "message count" random variable is Poisson distributed,
# and therefore lambda (the poisson parameter) is the expected value of
# "message count".
expected_texts_per_day[day] = (lambda_1_samples[ix].sum()
+ lambda_2_samples[~ix].sum()) / N
plt.plot(range(n_count_data), expected_texts_per_day, lw=4, color="#E24A33",
label="expected number of text-messages received")
plt.xlim(0, n_count_data)
plt.xlabel("Day")
plt.ylabel("Expected # text-messages")
plt.title("Expected number of text-messages received")
plt.ylim(0, 60)
plt.bar(np.arange(len(count_data)), count_data, color="#348ABD", alpha=0.65,
label="observed texts per day")
plt.legend(loc="upper left");
Our analysis shows strong support for believing the user's behavior did change ($\lambda_1$ would have been close in value to $\lambda_2$ had this not been true), and that the change was sudden rather than gradual (as demonstrated by $\tau$'s strongly peaked posterior distribution). We can speculate what might have caused this: a cheaper text-message rate, a recent weather-to-text subscription, or perhaps a new relationship. (In fact, the 45th day corresponds to Christmas, and I moved away to Toronto the next month, leaving a girlfriend behind.)
1. Using lambda_1_samples
and lambda_2_samples
, what is the mean of the posterior distributions of $\lambda_1$ and $\lambda_2$?
# type your code here.
2. What is the expected percentage increase in text-message rates? hint:
compute the mean of lambda_1_samples/lambda_2_samples
. Note that this quantity is very different from lambda_1_samples.mean()/lambda_2_samples.mean()
.
# type your code here.
3. What is the mean of $\lambda_1$ given that we know $\tau$ is less than 45. That is, suppose we have been given new information that the change in behaviour occurred prior to day 45. What is the expected value of $\lambda_1$ now? (You do not need to redo the PyMC part. Just consider all instances where tau_samples < 45
.)
# type your code here.
PyMC: Bayesian Stochastic Modelling in Python. Journal of Statistical Software, 35(4), pp. 1-81.
from IPython.core.display import HTML
def css_styling():
styles = open("../styles/custom.css", "r").read()
return HTML(styles)
css_styling()