작성자: admin 작성일시: 2017-11-24 13:21:53 조회수: 628 다운로드: 42
카테고리: 머신 러닝 태그목록:

몬테카를로 베이지안 분석

Rejection Sampling

당장 샘플링 가능한 유사 확률분포를 사용하여 목표 확률분포의 샘플을 생성

  • $p(x)$: 샘플링하고자 하는 목표 확률분포
  • $q(x)$: 샘플링 가능한 유사 확률분포
  • $k$: $kq(x) \geq p(x)$가 되도록 하는 스케일링 상수

$p(z)/kq(z)$의 확률로 $q(x)$의 샘플을 채택

In:
a = 2
b = 6
rv_p = sp.stats.beta(a, b)
rv_q = sp.stats.norm(loc=0.5, scale=0.5)
k = 5
In:
xx = np.linspace(0, 1, 100)
plt.plot(xx, rv_p.pdf(xx), 'r-')
plt.plot(xx, rv_q.pdf(xx) * k, 'g:')
plt.show()
In:
plt.plot(xx, rv_p.pdf(xx) / (rv_q.pdf(xx) * k), 'r-')
plt.show()
In:
np.random.seed(1)
x_q0 = rv_q.rvs(int(1e4))
x_q = x_q0[(x_q0 >= 0) & (x_q0 <= 1)]
crits = rv_p.pdf(x_q) / (rv_q.pdf(x_q) * k)
coins = np.random.rand(len(x_q))
x_p = x_q[coins < crits]
In:
sns.distplot(x_q, kde=False)
plt.show()
In:
sns.distplot(x_p, kde=False)
plt.show()
In:
y = np.random.rand(len(x_q)) * rv_q.pdf(x_q)
plt.plot(x_q, y, 'bs', ms=1)
ids = coins < crits
plt.plot(x_q[ids], y[ids], 'ro', ms=4)
plt.show()

Expectation Approximation Using Sampling

$$ \text{E}[f(X)] = \int f(x)p(x)dx $$
$$ \text{E}[f(X)] \approx \dfrac{1}{N} \sum_{i=1}^N f(x_i) $$

Importance Sampling

$$ \begin{eqnarray} \text{E}[f(X)] &=& \int f(x)p(x)dx \\ &=& \int f(x)\dfrac{p(x)}{q(x)} q(x) dx \\ &\approx & \dfrac{1}{N} \sum_{i=1}^N f(x_i)\dfrac{p(x_i)}{q(x_i)} \\ \end{eqnarray} $$
In:
mean = a / (a + b)
mean
Out:
0.25
In:
np.mean(x_p)
Out:
0.25317146797551493
In:
np.mean(x_q0 * rv_p.pdf(x_q0) / rv_q.pdf(x_q0))
Out:
0.25016824587488934
In:
len(x_p), len(x_q0)
Out:
(1982, 10000)

Conditional Probability and Markov Chain

$$ p(x^{(t + 1)}) = \sum_{x^{(t)}} p(x^{(t + 1)} \mid x^{(t)}) p(x^{(t)}) $$
  • detailed balance 조건을 만족하면 특정한 분포로 수렴

MCMC(Markov Chain Monte Carlo)

  • Markov Chain의 최종 수렴 분포가 원하는 분포 $p(x)$가 되도록 하는 Markov Chain을 만든다.

Metropolis-Hastings Sampling

  1. (step 0) $x^{(0)}$ 생성
  2. (step 1) 샘플링 가능한 $q(x^{\ast} \mid x^{(t)})$ 분포로 $x^{\ast}$ 생성
  3. (step 2) 다음 확률로 $x^{(t+1)} = x^{\ast}$ 선택, 선택되지 않으면 $x^{(t+1)} = x^{(t)}$ $$ \min \left( 1, \dfrac{p(x^{\ast}) q(x^{(t)} \mid x^{\ast})}{p(x^{(t)}) q(x^{\ast} \mid x^{(t)})} \right)$$
  4. 충분한 수의 샘플이 모일때까지 (step 1) ~ (step 2) 반복

Hamiltonian Monte Carlo

  • gradient 정보를 사용하여 리젝트 확률 감소

PyMC3

  • MCMC Bayesian in Python
  • Theano backend
In:
import pymc3 as pm
In:
cov = np.array([[1., 1.5], [1.5, 4]])
mu = np.array([1, -1])
In:
with pm.Model() as model:
    x = pm.MvNormal('x', mu=mu, cov=cov, shape=(1, 2))
    step = pm.Metropolis()
    trace = pm.sample(5000, step)
100%|██████████| 5500/5500 [00:02<00:00, 2663.97it/s]
In:
pm.traceplot(trace, skip_first=500)
plt.show()
In:
plt.scatter(trace['x'][500:, 0, 0], trace['x'][500:, 0, 1])
plt.show()