PyMCによるベイズ推論
本記事では以下のレポジトリの Ch1_Introduction_PyMC3.ipynb の内容を一部抜粋して実行してみる．
まず，データをダウンロードして表示する．
import numpy as np import matplotlib.pyplot as plt import os # 'data' ディレクトリを作成する makedirs('data', exist_ok=True) # サンプルデータ（ユーザが受信するメッセージ数）をダウンロードして 'data/txtdata.csv' に配置する from urllib.request import urlretrieve urlretrieve('https://git.io/vXTVC', 'data/txtdata.csv') # 可視化 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) 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)
このような計数データはポアソン分布を使ってモデリングする． 日目のメッセージ数を
とすると，
ここで，パラメータをどうやって決めるか．
グラフ後半のほうで数値が大きくなるようにみえる．
そこで，ある日
を境に，パラメータ
が突然大きくなると仮定する．
この2つの未知数について，ベイズ推論による推定を試みる．
ベイズ推論を使うには，2つの
に対して事前分布を決める必要がある．
正の実数のための確率分布として，指数分布を使うのがちょうどよい．
指数分布のパラメータをとおくと，
とかける．
については，変化点が何日目かを判断するのは難しい．そこで，ここでは一様分布を使う．
以上の確率変数 について，PyMCを使って推定する．
まず，に対応するPyMC変数を作成する．
import pymc3 as pm import theano.tensor as tt with pm.Model() as model: alpha = 1.0/count_data.mean() # count_data: メッセージ受信数 lambda_1 = pm.Exponential("lambda_1", alpha) lambda_2 = pm.Exponential("lambda_2", alpha) tau = pm.DiscreteUniform("tau", lower=0, upper=n_count_data - 1)
次に，関数 lambda_ を定義する．実際にはこれを確率変数 とみなすことができる．
with model: idx = np.arange(n_count_data) # Index lambda_ = pm.math.switch(tau > idx, lambda_1, lambda_2)
次に， count_data と lambda_ を受け取り，この例題のデータ生成モデルであるポアソン分布のオブジェクトを生成し，変数 observation に代入する．
with model: observation = pm.Poisson("obs", lambda_, observed=count_data)
以上のコードに対して，マルコフ連鎖モンテカルロ法(MCMC)を使った学習(事後分布からのサンプリング)を行う．
with model: step = pm.Metropolis() trace = pm.sample(10000, tune=5000,step=step)
出力は以下のようになる．
Multiprocess sampling (4 chains in 4 jobs) CompoundStep >Metropolis: [tau] >Metropolis: [lambda_2] >Metropolis: [lambda_1] Sampling 4 chains, 0 divergences: 100%|██████████| 60000/60000 [00:08<00:00, 7404.58draws/s] The number of effective samples is smaller than 25% for some parameters.
サンプル値のヒストグラムをプロットして，事後分布の形状を確かめる．
lambda_1_samples = trace['lambda_1'] lambda_2_samples = trace['lambda_2'] tau_samples = 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", density=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", density=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")
図から分かることは以下．
- 推定値の不確実さ
- 分布の幅が広ければ，事後信念はあまり確信できるものではない
- パラメータの妥当な値
はおよそ18，
はおよそ23
- 2つの
の事後分布が明らかに異なっている
- ユーザのメッセージ受信に変化があった可能性が高い
の分布
- 45日目にユーザが振る舞いを変えた確率が50%程度
最後に，サンプルデータにメッセージ数の期待値を重ねて表示する．
figsize(12.5, 5) N = tau_samples.shape[0] expected_texts_per_day = np.zeros(n_count_data) for day in range(0, n_count_data): ix = day < tau_samples 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")
