「ベイズ統計で実践モデリング」をPythonで実装する #1

【概要】

「ベイズ統計で実践モデリング」という書籍に掲載の事例をPythonで実装してみます
生で書くのではなく、PyMC、Pyroという二つのPPLを使ってを書いていきます

【目次】

はじめに
第三章二項分布を使った推論
おわりに
参考資料

はじめに

本シリーズでは、「ベイズ統計で実践モデリング」に記載の事例をPythonを使って実装します。利用するフレームワークはPyMC3とPyroです。なお、Pyroを使いますが、推論はMCMCを使っていく予定です。

第1回目は書籍の3章を扱います。

書籍の1,2章はベイズ推論の基礎とWinBUGSの説明です。そのため、本シリーズは3章から初めていきます。なお書籍ではRをベースに、WinBUGSをMCMC サンプラーとして、その実装も公開されています。サポートページはこちら。

筆者は特にPyroは初心者なので、間違いやより良い実装など指摘していただけるとめちゃくちゃ助かります。

なお、Pyroの記法については、以下のページでも扱っています。

learning-with-machine.hatenablog.com

【トップに戻る】

第三章二項分布を使った推論

3章では、最も基本的な確率分布のパラメータ推論として、二項分布のパラメータ推論を扱っています。

細かい実装については、例題ごとにnotebookでの実装をGithubで公開しているのでこちらを確認してください。

github.com

3.1 比率の推論

これは最も単純な二項分布のパラメータ推論です。3章は全体に解析解の導出も容易な問題ですが、導入ということですね。notebookはこちら。

二項分布のパラメータ推論なので、グラフィカルモデルとしては次の通りです。

f:id:hippy-hikky:20200922215931p:plain

$k$ は観測される値、 $n$ は試行回数です。これら二つの確率変数は所与のものという意味で青で塗りつぶしています。ここではなんらかの事象が発生する確率を表す確率変数 $\theta$ だけが未知であるというモデルです。

PyMC3でのモデルの定義は次のように書きました。

with pm.Model() as model:
    theta = pm.Beta('theta', 1, 1)
    y = pm.Binomial('y', p=theta, n=n, observed=k)

pyroの場合

def model_binom(n):
    theta = pyro.sample('theta', dist.Beta(torch.tensor(1.), torch.tensor(1.)))
    y = pyro.sample('y', dist.Binomial(total_count=n, probs=theta))
    return y

cond_model = pyro.condition(
    model_binom, 
    data={'y':torch.tensor(float(k))})

データは、n=10, k=5だけとして、推論させた結果は次の通りです。

f:id:hippy-hikky:20200922220426p:plain — 推論結果．

これはPyMC3での推論結果です。Pyroでもだいたい同じような結果になりました。

3.2 2つの比率の差

次は、発生確率 $\theta$ が異なる二つの対象から標本を採取したとして、その二つの対象の発生確率にどの程度の差異があるのかを推論します。notebookはこちら

グラフィカルモデル的には次の通りで、3.1のモデルが二つ結合したものになっています。

f:id:hippy-hikky:20200922221634p:plain

$\delta$ が二つの $\theta$ の差を表しています。

PyMC3では次のようにかけます。

with pm.Model() as model:
    theta1 = pm.Beta('theta1', 1, 1)
    theta2 = pm.Beta('theta2', 1, 1)
    
    y1 = pm.Binomial('y1', p=theta1, n=n1, observed=k1)
    y2 = pm.Binomial('y2', p=theta2, n=n2, observed=k2)
    
    delta = pm.Deterministic('delta', theta1 - theta2)

pm.Deterministicで確率変数の操作を定義することができるので、 $\delta$ をこのように書いています。

pyroの場合は次のようにしました。

def model_2binom(n1, n2):
    theta1 = pyro.sample('theta1', dist.Beta(torch.tensor(1.), torch.tensor(1.)))
    theta2 = pyro.sample('theta2', dist.Beta(torch.tensor(1.), torch.tensor(1.)))
    y1 = pyro.sample('y1', dist.Binomial(total_count=n1, probs=theta1))
    y2 = pyro.sample('y2', dist.Binomial(total_count=n2, probs=theta2))
    return y1, y2

cond_model = pyro.condition(
    model_2binom, 
    data={'y1':torch.tensor(float(k1)), 
          'y2':torch.tensor(float(k2))})

pymcのようにdeterminisiticな変数を定義できるのかと思ったんですが、書けないようでした。delta分布からサンプルするようにすれば良いのかとも思ったんですがそれも違うようです。なので、 $\theta_1, \theta_2$ のMCMCサンプルから別でdeltaに相当する分布を作っています。

結果は次のようになりました。

f:id:hippy-hikky:20200922222320p:plain — 推論結果．左がtheta_1、右がtheta_2．

f:id:hippy-hikky:20200922222450p:plain:w300 — 推論結果．二つの確率の差の事後分布．

3.3 共通の比率の推論

次は同一のモデルから標本を二つ取得した場合の例です。notebookはこちら。

グラフィカルモデルとしては3.1のものとほぼ同じですので省略します。

PyMC3でのモデル定義

with pm.Model() as model:
    theta = pm.Beta('theta', 1, 1)

    for i, (n, k) in enumerate(zip(ns, ks)):
        y = pm.Binomial(f'y_{i}', p=theta, n=n, observed=k)

$n, k$ のセットをリストで渡してくるため、forでループさせています。（この書き方が良いのかはあまり自信がない。。。）

次にPyroの場合です。

def model_binom(params=None):
    theta = pyro.sample('theta', dist.Beta(torch.tensor(1.), torch.tensor(1.)))
    y = pyro.sample('y', dist.Binomial(total_count=n, probs=theta))
    return y

cond_model = pyro.condition(
    model_binom, 
    data={'y':torch.tensor(ks).float(), 
          'n':torch.tensor(ns).float()})

pyroの場合は、モデル（model_binom）自体は3.1のものと全く同じです。与えるデータをリストで渡しているだけが違う点です。

推論結果は次の通りで、与えたデータセットは、nとkのペアが(10, 5)と(10,7)です。標本数が二つに増えているので事後分布の広がりが狭くなっています。

f:id:hippy-hikky:20200922223642p:plain — 推論結果．

3.4 事前/事後予測

確率モデルは、データの生成過程を表現しています。そして、既知の確率変数を条件として固定して、未知の確率変数に対する条件付き確率を推論するのが確率モデルでした。ここでは、3.1のモデルを使って、データが何も与えられていない場合のデータのサンプル（事前予測）とデータを得た後の事後分布を使った予測（事後予測）を行います。 notebookはこちら。

PyMCの場合は事前予測と事後予測は次のように書きます。

# 事前予測
with model:
    sample_pri = pm.sample_prior_predictive(samples=3000)

# 事後予測
with model:
    sample_post = pm.sample_posterior_predictive(trace)

このように、それぞれAPIが用意されています。

pyroの場合、事前予測は $n$ だけを条件として、サンプルしました。

cond_prior = pyro.condition(
    model_binom, 
    data={'n':torch.tensor(n).float()})

N_pri = 3000
prior_sample = [cond_prior().item() for i in range(N_pri)]

事後予測は、 $\theta$ のMCMCサンプルを条件として入力し、サンプルを得ましたが、これで合っているのかはあまり自信がありません。

cond_post = pyro.condition(
    model_binom, 
    data={'n':torch.tensor(n).float(), 
          'theta':mcmc.get_samples()['theta']})
post_sample = cond_post()

推論結果は次のようにまりました（PyMC3の結果）。

f:id:hippy-hikky:20200922224728p:plain — 推論結果．青が事前予測．オレンジが事後予測．

ちなみに、事後分布のために与えたデータは、n=15, k=1です。

事前予測の結果は、 $\theta$ の事前分布から一様分布状であると期待していましたがその通りになりました。事後予測は、0,1のあたりがピークになっていることがわかります。

3.5 事後予測

予測結果と観測されたデータが合致しない場合には、モデルの設定が何かおかしいと思うべきです。ここではモデルの適切さをチェックしてみます。 notebookはこちら。

問題としては3.2の二標本の問題になります。ですが、与えるデータを次のように極端にします。

ns = [10, 10]
ks = [0, 10]

3.1では二項分布のパラメータは共通であるとしていましたが、このデータは共通のパラメータを持ったモデルから生成されたとは考えにくいんじゃないかと思います。無理やり推論すると結果は次のようになりました。

f:id:hippy-hikky:20200922225625p:plain — 推論結果．

中間の0.5付近にピークがある分布です。

この結果に基づいて、3.4の事後予測を行った結果が以下です。

f:id:hippy-hikky:20200922225748p:plain — 事後予測結果．縦軸はy_0，横軸はy_1のサンプル結果．二つのデータがあるので2次元のヒートマップで表現．

この事後予測結果に対して、入力したデータはどこに位置するかと確認してみると、右上端です( $(y_0, y_1) = (0, 10)$ )。事後予測結果からすると非常に稀なケースであることがわかります。このことからも、推論語のモデルから期待されるデータと実際のデータがだいぶ解離していることがわかりますので、モデルの設計を見直した方が良いのかなと想像できます。確率的には0ではないので、このようなデータが得られないことはないのですが。

なお、この図は書籍の図と軸の向きが異なるので書籍と合わせて確認している場合は御注意ください。

3.6 同時分布

3章の最後では、これまで既知としていた試行回数に相当する $n$ も未知であり推論したいパラメータだとします。 notebookはこちら。

グラフィカルモデルとそれぞれの確率分布の定義は以下の通りです。

f:id:hippy-hikky:20200923233444p:plain

試行回数nは正の整数なので、カテゴリ分布からサンプルされるものとします。また、nの最大値だけはわかっているものとします。

こちらのモデルですが、以下のように実装したのですが、これで合っているのか自信がありません。詳しい方指摘いただけると助かります。

n_max = 500
with pm.Model() as model:
    n_ = pm.Categorical('n_', np.ones(n_max)/n_max)
    n = pm.Deterministic('n', n_+max(ks))
    
    theta = pm.Beta('theta', 1, 1)

    for i, k in enumerate(ks):
        y = pm.Binomial(f'y_{i}', p=theta, n=n, observed=k)