１ 初期値$μ^{(0)}$を決め、t=1とおく。</p>
２ 現在μ^(t-1)であるとき、次の点μ^tの候補μ'をμ'~N(x_bar, 1/n)により発生させて、$$ α(μ^{(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として２にもどる。
In [5]:
%matplotlib inline
import numpy as np
import math
import matplotlib.pyplot as plt
観測されたデータ。N(5,1)から作られた10個の乱数。
In [6]:
X = np.zeros(10)
for i in range(len(X)):
X[i] = np.random.normal(5,1)
In [7]:
X
Out[7]:
[ベイズ推論]</P>
確率変数X_1, X_2,..., X_nが互いに独立に平均がμ、分散が1であるような正規分布に従うとする。
μの事前分布にt分布を仮定する。
In [28]:
class IMH:
def __init__(self, X):
self.mu_0 = 1
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):
mu = np.random.normal(self.x_var, 1.0/10.0)
return 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_0
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]
return mu
In [29]:
imh = IMH(X)
In [30]:
mu = imh.simulate()
In [31]:
plt.plot(mu)
Out[31]:
In [32]:
plt.hist(mu)
Out[32]:
In [33]:
plt.hist(mu[1000:])
Out[33]:
$\mu$の事前分布に無情報事前分布を仮定する。
In [22]:
class IMH:
def __init__(self, X):
self.mu_0 = 1
self.x_var = np.mean(X)
def prior_dist(self, t):
ft = np.random.normal(0, 10000)
return ft
def prop_dist(self):
mu = np.random.normal(self.x_var, 1.0/10.0)
return 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_0
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]
return mu
In [23]:
imh = IMH(X)
In [24]:
mu = imh.simulate()
In [25]:
plt.plot(mu)
Out[25]:
In [26]:
plt.hist(mu)
Out[26]:
In [27]:
plt.hist(mu[1000:])
Out[27]:
In [ ]: