PyMC3를 활용한 샘플링
동전 던지기 예시
prior(사전분포), likelihood(가능도), posterior(사후분포)가 어떤 관계를 갖는지 실제 파이썬 코드로 구현해보며 확인해보자. 파이썬에서 사용할 수 있는 베이즈 최적화 모듈은 PyMC가 있다. 파이썬 3에 최적화된 PyMC3를 사용한다. PyMC2와는 어느정도 차이가 있으므로, 주의하며 사용하자. 동전 던지기 (coin flipping) 예시는 결과가 0아니면 1만 나오는 Binomial 사건이므로, prior로 beta prior를 택했다.
PyMC3 모듈을 사용해, 사후분포와 가능도가 변할 때 어떻게 사후분포가 변하는 지 확인할 수 있다.
아래의 모든 예시들은 PyMC3의 pm.find_MAP() 최적화 접근과 pm.Metropolis() 샘플링 방법론을 사용했다. 이 방법론들에 대한 구체적인 이야기는 이전 이론 노트 페이지에 설명 돼 있다. 추후에 코드와 함께 또 다시 설명할 예정이다. 이번 예시에서는 사전분포, 가능도, 사후분포들의 관계만 시각적으로 확인한다.
MCMC는 사전분포와 가능도의 곱인 \(P(\theta)P(D \mid \theta)\)를 계산하고, random walk의 반복을 통해 사후분포를 추정한다.
n = 100
h = 2
alpha = 1
beta = 1
niter = 1000
동전을 던지는 전체 시행 횟수는 100, 그 중에서 앞면은 극단적으로 적게 2번이 나왔다고 가정해보자. alpha와 beta는 prior 분포를 결정한다. 구체적인 이야기는 아래 그래프를 보며 더 설명하겠다. sampling iteration 횟수는 1000번으로 정했다.
plt.figure(figsize = (8,6))
plt.hist(trace['p'], 15, histtype = 'step', normed = True, label = 'post');
x = np.linspace(0,1,100)
plt.plot(x, stats.beta.pdf(x, alpha, beta), label = 'prior');
plt.legend(loc='best');
plt.show()
alpha = beta = 1 인 탓에, 사전분포는 uniform 분포 형태를 취하고 있으며, 이는 데이터가 우리에게 주어지기 전에 어떠한 정보도 포함하고 있지 않음을 얘기한다. 이 때문에, 그래프가 보여주는 사후분포는 우리의 데이터를 따라 매우 왼쪽으로 편향 돼있다. 데이터에서 100번의 시도 중 2번만 앞면이 나왔으니, 데이터를 보고 나서 우리의 믿음인 사후분포는 $\theta$값이 상당히 작을 것이라고 믿게 된 것이다. 이제 여기서 사전분포가 정보를 포함하도록 변경해보자.
n = 100
h = 2
alpha = 50
beta = 50
niter = 1000
alpha = beta = 50 으로 변경이후, 사전분포는 $p$ 값이 0.5 부근에 존재할 것이라는 믿음을 어느정도의 확실성으로 보여준다.
사전분포가 0.5 부근에서 꽤 뾰족하게 나타나고 있다. 알파 베타 값이 각각 20으로 같았다면, 이 분포는 조금 덜 뾰족할 것이다. 즉, 불확실성이 높았을 것이다.
주어진 데이터는 변하지 않았는데, 사전분포의 정보가 가운데 값에서 어느정도의 믿음을 보여준 탓에, 우리의 사후 분포도 이전 케이스보다 조금 더 오른쪽으로 이동하게 됐다. 이번에는 사전분포는 그대로 두고, 데이터를 바꿔보자. 어느정도 예상 가능한 결과가 나올 것이다.
n = 100
h = 98
alpha = 50
beta = 50
niter = 1000
사전분포는 그대로 alpha = beta = 50, 주어진 데이터는 100회 중 98번이 앞면이 나온다. 사전분포는 그대로지만, 주어진 데이터가 나올 가능도는 $p$ 값이 1에 가까운 쪽으로 치우진다. 그 결과로 사후 분포 또한 오른쪽으로 이동하게 됐다.
정규분포의 평균과 표준편차 값을 추론해보자
$\ $ 정규분포 \(X \thicksim N(\mu, \sigma^2)\) 에 대하여 평균과 표준편차에 대한 추정을 진행해보자.
$\ $ 총 100개의 데이터 셋을 생성한다. 평균은 10, 표준편차는 2로 잡고, 이 정규분포에서 총 100개의 데이터 셋을 생성한다. (observed Data Set)
우리가 추정하고자 하는 $\mu$와 $\sigma$에 대해, 가능도 함수는 Normal distribution이므로, pm.Normal()을 가능도 함수로 설정한다. 샘플링 함수는 pm.NUTS()를 사용했으며, 이에 대한 설명은 잠시 스킵한다.
# generate observed data
N = 100
_mu = np.array([10])
_sigma = np.array([2])
y = np.random.normal(_mu, _sigma, N)
niter = 1000
with pm.Model() as model:
# define priors
mu = pm.Uniform('mu', lower=0, upper=100, shape=_mu.shape)
sigma = pm.Uniform('sigma', lower=0, upper=10, shape=_sigma.shape)
# define likelihood
y_obs = pm.Normal('Y_obs', mu=mu, sd=sigma, observed=y)
# inference
start = pm.find_MAP()
step = pm.Slice()
trace = pm.sample(niter, step, start, random_seed=123, progressbar=True)
plt.figure(figsize=(10,4))
plt.subplot(1,2,1);
plt.hist(trace['mu'][-niter/2:,0], 25, histtype='step');
plt.subplot(1,2,2);
plt.hist(trace['sigma'][-niter/2:,0], 25, histtype='step');
N = 100
_mu = np.array([10])
_sigma = np.array([2])
y = np.random.normal(_mu, _sigma, N)
niter = 1000
with pm.Model() as model:
## prior distributions
mu = pm.Uniform('mu', lower = 0, upper = 20, shape = _mu.shape)
sigma = pm.Uniform('sigma', lower=0, upper = 3, shape = _sigma.shape)
## likelihood
y_obs = pm.Normal('Y_obs', mu = mu, sd=sigma, observed = y)
# inference
start = pm.find_MAP()
step = pm.NUTS()
trace = pm.sample(niter, step, start, random_seed = 123, progressbar = True)
pm.traceplot(trace);
plt.subplot(1,2,1);
plt.hist(trace['mu'], bins = 50, alpha = 0.5, color = 'orange');
plt.title("Inference for Mu", fontsize = 15)
plt.subplot(1,2,2);
plt.title("Inference for Sigma", fontsize = 15)
plt.hist(trace['sigma'], bins = 50, alpha = 0.5, color = 'orange');
선형회귀 모델의 파라미터들을 추론해보자
$\ $우리가 친숙하게 알고있는 선형 회귀 모델의 파라미터들을 추정해보자.
다음과 같이 error term 을 포함한 식으로 표현할 수 있다.
이는 다음의 확률분포에서 $y$ 값에 대한 샘플링으로 이해할 수 있다.
이제, PyMC3 를 이용해서 해당 회귀 모델의 파라미터들인 \(a,b,\sigma\)에 대한 추론을 진행하자. 각 파라미터들에 사전분포는 다음과 같이 정의하고 시작한다.
Likelihood 함수는 정규분포 확률함수, 샘플러는 마찬가지로 pm.NUTS 를 사용했다.
# observed data
n = 11
_a = 6
_b = 2
x = np.linspace(0, 1, n)
y = _a*x + _b + np.random.randn(n)
n_iter = 100
with pm.Model() as model:
a = pm.Normal('a', mu=0, sd=20)
b = pm.Normal('b', mu=0, sd=20)
sigma = pm.Uniform('sigma', lower=0, upper=20)
y_est = a*x + b # simple auxiliary variables
likelihood = pm.Normal('y', mu=y_est, sd=sigma, observed=y)
# inference
start = pm.find_MAP()
step = pm.NUTS() # Hamiltonian MCMC with No U-Turn Sampler
trace = pm.sample(niter, step, start, random_seed=123, progressbar=True)
pm.traceplot(trace);
