FC2ブログ
11«1.2.3.4.5.6.7.8.9.10.11.12.13.14.15.16.17.18.19.20.21.22.23.24.25.26.27.28.29.30.31.»01

たまひびとらの絵本の実

読書好きな姉妹と弟と父母の読んだ本

Pythonで体験するベイズ推論 

bayesianmethodsforhackers.jpg
Pythonで体験するベイズ推論
Pythonで体験するベイズ推論

本書はpymc2で書かれているのでpymc2をインストールして、githubのpymc2の方からダウンロードしたサンプルコード・データを使う。普段プログラミングをする人たちには当たり前でも、素人にはこういうのがハードルなんだよね。

そもそも理解できていないところ。尤度を求めるには、尤度関数(密度関数と同じ)が必要。でもパラメーターがどんな分布に従うのかわからなければ、どの尤度関数を使えばいいのかわからない。たとえば離散で範囲限定なら二項分布を尤度関数として尤度を計算すればよい。しかし、範囲が限定されているのかわからないケースもあるだろう。連続だとしても、正規分布なのかガンマ分布なのかまた他の密度関数なのか、どうやってわかるんだろう。このリンクによれば、従う分布をわかる必要があるというような書き方。
(後日: ベイズ推論では尤度関数が何になるかを想定して計算してみる、という手順らしい。つまり尤度関数が何になるかは自分である程度アイディアを持つ必要がある。)

まあごちゃごちゃ考えるよりやってみるべき、と思ったので本書の例をちょっとだけ変えて MCMCをやってみたのが①。1日当たり何人の知り合いに会うかを、最初はA町で、どこかの時点からはB町で数える。日数は全部で100日。何人に会うかはポアソン分布(尤度もポアソンで計算しているのだと思う)。A町での人数のポアソン分布のラムダ、B町での人数のポアソン分布のラムダ、移動日、の分布の推定をする、というもの。やってみるとわかるけど、尤度は確率ではないので合計は1にならない。

②は、事前分布が、事前分布×尤度=事後分布になる様子の例。尤度の式は、ポアソン分布で確率の積を求めている。つまり、事前分布で確率がゼロであれば、事後でも確率はゼロとなってしまう。

③では pymc のobservation 文の尤度関数を変えるとどうなるのかを実験。ベルヌーイとポアソンでヒストグラムを作ってみると、やっぱり事後分布に違いがあるみたい。ベルヌーイ尤度の方が、事後分布の分散が大きそう。

後半に出てくる共益事前分布: 事前分布をベータ分布、尤度をベルヌーイ分布とすると、事後分布はベータ分布になる。共益事前分布を使えば、MCMCなしに、事後分布を直接計算できる。他にも同様の組合せがある。事前分布をベータ分布、尤度を二項分布とすると、事後分布はベータ分布。事前分布をディリクレ分布(ベータ分布の多変量バージョン)、尤度を多項分布(二項分布の多変量バージョン)とすると、事後分布はディリクレ分布。共益事前分布は少なくともベータ + 二項分布のケースは簡単。単に確率の掛け算をすれば、掛け算の結果もベータ分布になることがわかる。


import pymc as pm
import numpy as np
from IPython.core.pylabtools import figsize
from matplotlib import pyplot as plt
%matplotlib inline

# create sample data
T = np.random.randint(20, 80)
L1 = np.random.randint(10, 25)
L2 = np.random.randint(5, 10)
print('true L1',L1)
print('true L2',L2)
print('true T',T)
data = np.r_[pm.rpoisson(L1, T), pm.rpoisson(L2, 100-T)]
data_count = 100
data_sum = np.sum(data)

# prior probability
beta = 1.0 / (data_sum/data_count)
L1 = pm.Exponential("L1", beta)
L2 = pm.Exponential("L2", beta)
T = pm.DiscreteUniform("T", lower=0, upper=data_count)

@pm.deterministic
def Seek(T=T, L1=L1, L2=L2):
out = np.zeros(data_count)
out[:T] = L1
out[T:] = L2
return out

# modelling
observation = pm.Poisson("obs", Seek, value=data, observed=True)
model = pm.Model([observation, L1, L2, T])

# sampling
mcmc = pm.MCMC(model)
mcmc.sample(20000,10000)

# show posterior
L1_samples = mcmc.trace('L1')[:]
L2_samples = mcmc.trace('L2')[:]
T_samples = mcmc.trace('T')[:]

figsize(15, 15)
ax = plt.subplot(311)
plt.title(r"""Posterior distributions of the variables L1, L2, T""")
plt.hist(L1_samples, histtype='stepfilled', bins=20, label="posterior distribution of L1", color="blue", density=True)
plt.legend(loc="upper left")
plt.xlim([0, 50])
plt.ylim([0, 1])

ax = plt.subplot(312)
plt.hist(L2_samples, histtype='stepfilled', bins=20, label="posterior distribution of L2", color="blue", density=True)
plt.legend(loc="upper left")
plt.xlim([0, 50])
plt.ylim([0, 1])

plt.subplot(313)
w = 1.0 / T_samples.shape[0] * np.ones_like(T_samples)
plt.hist(T_samples, histtype='stepfilled', bins=10, label="posterior distribution of T", color="red", density=True)
plt.legend(loc="upper left")
plt.ylim([0, 1])
plt.xlim([0, 100])


import scipy.stats as stats
from IPython.core.pylabtools import figsize
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
%matplotlib inline

figsize(15, 15)

N = 1
lambda_1_true = 0.05
lambda_2_true = 0.03

data = np.concatenate([
stats.poisson.rvs(lambda_1_true, size=(N, 1)),
stats.poisson.rvs(lambda_2_true, size=(N, 1))], axis=1)

x = y = np.linspace(0, 5, 10)
likelihood_x = np.array([stats.poisson.pmf(data[:, 0], _x) for _x in x]).prod(axis=1)
likelihood_y = np.array([stats.poisson.pmf(data[:, 1], _y) for _y in y]).prod(axis=1)
L = np.dot(likelihood_x[:, None], likelihood_y[None, :])

fig = plt.figure()
X, Y = np.meshgrid(x, y)

uni_x = stats.uniform.pdf(x, loc=0, scale=5)
uni_y = stats.uniform.pdf(x, loc=0, scale=5)
M = np.dot(uni_y[:, None], uni_x[None, :])

ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, M, cmap='bwr', vmax=1, vmin=-.15)
ax.plot_surface(X, Y, M * L, cmap='plasma')
ax.view_init(azim=200)
plt.title("Uniform Prior Distribution and Posterior Distribution");


import numpy as np
from IPython.core.pylabtools import figsize
import matplotlib.pyplot as plt
import pymc as pm
%matplotlib inline

p_true = 0.4
x = np.random.binomial(1, p_true, size=1000)

p = pm.Uniform('p', 0.0, 1.0)

figsize(15, 5)

obsB = pm.Bernoulli('obsB', p, value=x, observed=True)
mcmcB = pm.MCMC([p, obsB])
mcmcB.sample(10000, 2000)
p_traceB = mcmcB.trace('p')[:]
plt.hist(p_traceB, bins=30, histtype='stepfilled', alpha=0.2, color='red', density=True)

obsP = pm.Poisson('obsP', p, value=x, observed=True)
mcmcP = pm.MCMC([p, obsP])
mcmcP.sample(10000, 2000)
p_traceP = mcmcP.trace('p')[:]
plt.hist(p_traceP, bins=30, histtype='stepfilled', alpha=0.5, density=True)

スポンサーサイト

category: 父の本

cm 0   tb 0   page top

コメント

page top

コメントの投稿

Secret

page top

トラックバック

トラックバックURL
→http://tamahibi.blog122.fc2.com/tb.php/3457-4f4e9b87
この記事にトラックバックする(FC2ブログユーザー)

page top

訪問者数

カテゴリ

最新記事

最新コメント

最新トラックバック