In [1]:
%matplotlib inline
import numpy as np
import math
import matplotlib.pyplot as plt
観測されたデータ。N(5,1)から作られた10個の乱数。
In [2]:
X = np.zeros(10)
for i in range(len(X)):
X[i] = np.random.normal(5,1)
In [3]:
X
Out[3]:
[ベイズ推論]</P>
確率変数$X_1, X_2,..., X_n$が互いに独立に平均がμ、分散が1であるような正規分布に従うとする。
μの事前分布にt分布を仮定する。
１ 初期値μ^(0)を決め、t=1とおく。</p>
２ 現在μ^(t-1)であるとき、次の点μ^tの候補μ'を$μ'$~ $N(μ^{(t-1)}, r^2)$により発生させて、$$ α(μ^{(t-1)}, μ'|x) = min \{ 1, \frac{π(μ')}{π(μ^{(t-1)})}\}$$とおく。
３ (0, 1)上の一様乱数uを発生させて、$$μ^{(t)}=μ' if u≤α(μ^{(t-1)}, μ'|x)のとき$$ $$μ^{(t)}=u^{(t-1)} if μ^{(t-1)}>α(μ^{(t-1)}, μ'|x)のとき$$とする。
４ tをt+1として２にもどる。
事前分布に自由度5のt分布を仮定する。r^2=0.001とする。
In [28]:
class RWMH:
def __init__(self, X):
self.mu = 2
self.freedom = 5.0
self.x_var = np.mean(X)
def prior_dist(self, t):
ft = math.gamma((self.freedom+1.0)/2.0)/(math.sqrt(self.freedom*math.pi)*math.gamma((self.freedom)/2.0))*pow(1+pow(t, 2)/self.freedom, -(self.freedom+1)/2.0)
return ft
def prop_dist(self):
self.mu = np.random.normal(self.mu, 0.5)
return self.mu
def accept(self, mu_new, mu):
return min([1, self.prior_dist(mu_new)/self.prior_dist(mu)])
def simulate(self):
mu = np.zeros(110000)
mu[0] = self.mu
for i in range(1,110000):
mu_new = self.prop_dist()
u = np.random.uniform()
if u <= self.accept(mu_new, mu[i-1]):
mu[i] = mu_new
else:
mu[i] = mu[i-1]
self.mu = mu[i]
return mu
In [29]:
rwmh = RWMH(X)
In [30]:
mu = rwmh.simulate()
In [31]:
plt.plot(mu)
Out[31]:
In [32]:
plt.hist(mu)
Out[32]:
In [16]:
plt.hist(mu[10000:])
Out[16]:
In [18]:
mu.mean()
Out[18]:
$\mu$の事前分布に無情報事前分布を仮定する。
In [10]:
class RWMH:
def __init__(self, X):
self.mu = 5
self.x_var = np.mean(X)
def prior_dist(self, t):
ft = np.random.normal(0, 10000)
return ft
def prop_dist(self):
self.mu = np.random.normal(self.mu, 0.001)
return self.mu
def accept(self, mu_new, mu):
return min([1, self.prior_dist(mu_new)/self.prior_dist(mu)])
def simulate(self):
mu = np.zeros(11000)
mu[0] = self.mu
for i in range(1,11000):
mu_new = self.prop_dist()
u = np.random.uniform()
if u <= self.accept(mu_new, mu[i-1]):
mu[i] = mu_new
else:
mu[i] = mu[i-1]
self.mu = mu[i]
return mu
In [26]:
rwmh = RWMH(X)
In [27]:
mu = rwmh.simulate()
In [28]:
plt.plot(mu)
Out[28]:
In [29]:
plt.hist(mu)
Out[29]:
In [30]:
plt.hist(mu[1000:])
Out[30]:
In [ ]: