[이재우] Bayesian Parameter Estimation

# Introduction

Bayesian parameter estimation 이라고 오늘 공부할 주제라고 딱 들으면 누구나 하기 싫어질것 같다. 왠지 많은양의 Bayesian이 쓴 논문을 읽어야 할 것 같고, 외국사람이니 또 영어로 읽을 생각하니 벌써 머리가 아프다. 나도 하기 싫다. 뭐든 공부라고 생각이 들면, 반발심으로 하기 싫다는 말이 먼저 머릿속에 떠오른다. 그래서 나는 Bayesian parameter estimation을 공부해보는 것이 아닌 누군가 엄청 똑똑한 사람이 많은양의 논문을 읽고 생각하고 고민하고 이해해서 쓴글을 그저 따라하면서 효율적으로 이해해보려고 한다. [*Hands On Bayesian Statistics with Python, PyMC3 & ArviZ*](https://towardsdatascience.com/hands-on-bayesian-statistics-with-python-pymc3-arviz-499db9a59501) 라는 제목으로 *Susan Li*님께서 쓴 글을 번역을 기본으로 하며, 직관적으로 어려운 것들은 나의 말로 다시 설명하여 이해해보자. 신기하게도 이 글에서도 다음과 같은 문구로 시작함을 알 수 있다.

* *I believe that for the things we have to learn before we can do them, we learn by doing them. And nothing in life is so hard that we can’t make it easier by the way we take it.*

주구장창 논문만 읽으면서 이해하는 것은 누구에게나 어렵구나라는 위안과 함께 시작하자.
# Warm-up
이 글에서 Bayesian estimation 을 위해서 사용한 2가지 라이브러리가 있다. 본문에서는 다음과 같이 소개한다.

* *PyMC3 is a Python library for probabilistic programming with a very simple and intuitive syntax. ArviZ, a Python library that works hand-in-hand with PyMC3 and can help us interpret and visualize posterior distributions.*

Bayesian estimation을 하려면 이 두가지의 library를 이용해야함을 알고 넘어가자.

이 글에서는 자주 Bayesian estimation 과 다른 방법을 비교하는 글이 다음과 같이 자주나온다.

* *Bayesian models are also known as probabilistic models because they are built using probabilities. And Bayesian’s use probabilities as a tool to quantify uncertainty. Therefore, the answers we get are distributions not point estimates.*

언뜻 보면 어렵지만, 가장 단순하게 해석하면 Bayesian estimate 와 point estimate를 비교한 것이라고 할 수 있고, 불확실성을 다루기 위해서 분포 또는 확률이라는 개념을 사용하는 것 같다. 그러면서 다음과 같은 기본 단계를 설명해준다.
* *Step 1: Establish a belief about the data, including Prior and Likelihood functions.*
* *Step 2, Use the data and probability, in accordance with our belief of the data, to update our model, check that our model agrees with the original data.*
* *Step 3, Update our view of the data based on our model.*

해석을 하라면 할 순 있지만 익숙하지 않은 단어를 이해하지 못한 상태에서의 해석은 의미 없으므로, `belief`, `Prior` 과 `Likelihood function` 이것이 무엇일까 라는 궁금을 가지면 될 것 같다.

# Data
본문에서는 다음과 같이 시작한다.

* *Since I am interested in using machine learning for price optimization, I decide to apply Bayesian methods to a Spanish High Speed Rail tickets pricing data set*.

이를 통해서 우리가 예측할 것은 기차푯값을 예측함을 알 수 있다.
```python
from scipy import stats
import arviz as az
import numpy as np
import matplotlib.pyplot as plt
import pymc3 as pm
import seaborn as sns
import pandas as pd
from theano import shared
from sklearn import preprocessing
print('Running on PyMC3 v{}'.format(pm.__version__))
data = pd.read_csv('https://raw.githubusercontent.com/susanli2016/Machine-Learning-with-Python/master/data/renfe_small.csv')
data.head(2)
```
데이터를 읽어올 때, 가장 먼저 `head`를 통해 잘 불러왔는지 확인한다.

```py
data.isnull().sum()/len(data)
```
> insert_date 0.000000
origin 0.000000
destination 0.000000
start_date 0.000000
end_date 0.000000
train_type 0.000000
price 0.119467
train_class 0.003993
fare 0.003993
dtype: float64

위의 결과를 통해서 우리는 데이터 세트의 결측치가 존재함을 확인했다. 여기서 우리가 주목해봐야 할 것은 어떻게 이 NaN의 값을 처리하였는가이다.

```py
data['train_class'] = data['train_class'].fillna(data['train_class'].mode().iloc[0])
data['fare'] = data['fare'].fillna(data['fare'].mode().iloc[0])
data['price'] = data.groupby('fare').transform(lambda x: x.fillna(x.mean()))
```
본문에서는 따로 NaN 의 처리에 대해서 깊은 설명은 존재하지 않는다. 하지만, 설명을 하지 않아도 왜 그렇게 사용하였는지는 간단히 알 수 있다. 예를들면, **train_class** 열차의 등급을 나타내고, 이 열의 해당하는 NaN 값은 최대 빈도인 `mode`를 사용했다. `mode`를 사용한 이유는 아무래도 높은 등급에 타는 사람들보다 일반 등급에 기차를 타는 사람인 확률이 높으므로 사용했을 것 같다. 이처럼 categorical data의 imputation 은 모두 mode로 처리하였고, 유일한 numerical value인 price 만 `mean`을 사용했다. 이러한 방법을 사용한 근거는 정확히 나와있지 않지만, 대게는 그럴 가능성이 높다가 맞는것 같다. Imputation에 대해서는 다양한 의견이 존재하고 이것은 맞고 저것은 틀리다라고 말할 수 없다. 내가 만약에 imputation을 한다면 좀 더 명확한 근거가 있는 방법을 사용해보도록 해야겠다.
```py
az.plot_kde(data['price'].values, rug=True)
plt.yticks([0], alpha=0);
```

`az.plot_kde` 는 주어진 data에서의 observed 된 값들의 distribution을 보여준다. 원문에서는 위의 그래프를 다음과 같이 설명하였다.

* *The KDE plot of the rail ticket price shows a Gaussian-like distribution, except for about several dozens of data points that are far away from the mean.*

나는 덜 똑똑해서 단번에 Gaussian 분포를 예상 할 순 없었다. 하지만, 여기서 확인할 수 있는 사실은 기차푯값 데이터의 분포는 약간 right skewed distribution처럼 보기도하고 만약 오른쪽 끝 부분을 자르고 왼쪽 부분의 데이터가 더 있다고 가정하면 Normal distribution처럼 보일 수 있다. 참고로, Gaussian 분포는 Normal 분포와 되게 비슷하다. 그러니 가장 안전한, Normal 분포를 사용한건 아닌지 감히 예측해본다. 따라서, 기차푯값의 분포는 Gaussian 분포라고 가정하자.

# Modeling

나의 짧은 학사 경험으로는 Modeling이 가장 복잡하다. 이유는 내가 왜 이 modle을 사용하는지 이 글을 접하는 모든 사람에게 타당성을 제공 할 수 있어야 하기 때문이다. 원문을 확인해보자.
* *We would instantiate the Models in PyMC3 like this:*
*Model specifications in PyMC3 are wrapped in a with-statement.*
* *Choices of priors:*
* *μ, mean of a population. Normal distribution, very wide. I do not know the possible values of μ, I can set priors reflecting my ignorance. From experience I know that train ticket price can not be lower than 0 or higher than 300, so I set the boundaries of the uniform distribution to be 0 and 300. You may have different experience and set the different boundaries. That is totally fine. And if you have more reliable prior information than I do, please use it!*
* *σ, standard deviation of a population. Can only be positive, therefore use HalfNormal distribution. Again, very wide.*
* *Choices for ticket price likelihood function:*
* *y is an observed variable representing the data that comes from a normal distribution with the parameters μ and σ.*
* *Draw 1000 posterior samples using NUTS sampling.*

이 부분의 *Step 1* 에 해당한다. 먼저, **PyMC3** 라는 library를 통해서 modeling 단계를 시작 한다. 모델 설정을 위해서는 `prior`와 `likelihood`를 설정해야한다. 앞에서 우리가 몰랐던 단어들을 알 수 있는 기회이다. `prior`는 `likelihood function`을 정의하기 위해서 사용되는 parameter의 개념이고, `prior`는 각각의 분포를 가지고 있으며 우리의 `belief`를 기반으로 설정을 한다. 그리고 `likelihood function`은 우리가 위에서 예측한 기차푯값의 분포를 의미하고 정의된 `prior`를 바탕으로 정의 된다. Gaussian distribution을 정의하기 위해서는 두 개의 파라미터 μ, σ가 필요하다. 따라서 우리가 μ, σ 값을 예측할 수 있게 되면 기찻값의 Population Distribution을 예측할 수 있게 된다. 여기서는 prior 설정을 경험을 바탕으로 설정하였지만, prior 설정은 절대적인 기준이 있는것이 아닌, 더 나은 정보들이 존재하면 그것을 사용하여 prior를 설정해도 상관없다.

```py
with pm.Model() as model_g:
μ = pm.Uniform('μ', lower=0, upper=300)
σ = pm.HalfNormal('σ', sd=10)
y = pm.Normal('y', mu=μ, sd=σ, observed=data['price'].values)
trace_g = pm.sample(1000, tune=1000)
```
>Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Sequential sampling (2 chains in 1 job)
NUTS: [σ, μ]
Sampling chain 0, 0 divergences: 100%|██████████| 2000/2000 [00:06<00:00, 328.64it/s]
Sampling chain 1, 0 divergences: 100%|██████████| 2000/2000 [00:04<00:00, 445.79it/s]

원문에서 *y*를 다음과 같이 설명했다. *The y specifies the likelihood. This is the way in which we tell PyMC3 that we want to condition for the unknown on the knows (data).* 내가 이해한 바와 이 분의 설명이 같음을 다시 확인 할 수 있다.

지금까지 우리가 가정한 사실들을 바탕으로 PyMC3를 사용하여 Model 설정하였다. 그리고 우리가 세운 모델로 1,000개의 표본 뽑기를 2번 했다는 결과도 확인할 수 있다. **Model 설정에서 중요한 사실은 parameter 들마다 각각의 distribution을 갖고 있다는 것이다.** 이것이 Bayesian estimation의 point estimation 과 비교했을때의 특징임을 알 수 있다. 지금 우리가 가진 데이터로 어떠한 결론을 낼 수 있다고 하더라도, 나중에 다른 데이터를 받아들여서 만약 그 결과가 바뀌게 될 수도 있다는 사실(=Uncertainty)을 내포하고 있다.

```py
az.plot_trace(trace_g);
```

위의 그래프의 설명은 본문의 내용이 너무나 직관적으로 설명을 해주고 있어서 그대로 가져와서 사용할 것이다.
* *On the left, we have a KDE plot, — for each parameter value on the x-axis we get a probability on the y-axis that tells us how likely that parameter value is.*
* *On the right, we get the individual sampled values at each step during the sampling. From the trace plot, we can visually
get the plausible values from the posterior.*
* *The above plot has one row for each parameter. For this model, the posterior is bi-dimensional, and so the above figure is showing the marginal distributions of each parameter.*

간단히 말하면,

* 왼쪽은 2 chain of sample의 `kde` distribution을 각 parameter에 대해 보여주고 있다.

* 오른쪽은 각 표본 뽑기 마다 실질적인 parameter의 값을 step by step 보여주고 있다.

* 각 행은 하나의 parameter에 대해서 marginal distribution을 보여준다.

여기서 말하는 *marginal distribution*은 여러 parameter의 관계없이 하나의 parameter에 대해서만 보여주는 distribution이다.
각각의 parameter는 비슷한 shape를 보여주고 크게 눈에 띄는 이상한 모양이나 값이 없음을 확인할 수 있다. 이 말은, 더 많은 표본 뽑기(=chain)을 만들게 되면 각 parameter의 가장 peak 인부분으로 converge 함을 알 수 예상 할 수 있다. 여기서 중요한 사실은 convergence이다. 앞서 prior 설정은 크게 상관없다는 이유를 여기서 알 수 있다. 출발을 어디서 하든, 만약 우리가 추론하려는 parameter가 converge 한다는 사실을 알면 결국 그 끝에는 같은 결과를 마주하게 되기 때문이다.

```py
az.plot_joint(trace_g, kind='kde', fill_last=False);
```

위의 estimate는 2개의 parameter 를 예측 하기 때문에, 2 dimension prediction 이고, 위의 그래프 처럼 joint graph 를 통해서 서로 다른 2개의 parameter 는 independent 함을 알 수 있다. 왜냐하면 각 각의 peak에서 서로의 best case 를 확인 할 수 있기 때문이다.

```py
az.summary(trace_g)
```
mean | sd | hpd_3% | hpd_97% | mcse_mean | mcse_sd | ess_mean | ess_sd | ess_bulk | ess_tail | r_hat
| -- | -- | -- | -- | -- | -- | -- | -- | -- | -- | --
62.377 | 1.420 | 59.679 | 64.981 | 0.032 | 0.023 | 1925.0 | 1925.0 | 1928.0 | 1569.0 | 1.0
23.581 | 1.033 | 21.678 | 25.583 | 0.026 | 0.018 | 1630.0 | 1626.0 | 1641.0 | 1334.0 | 1.0

`az.summary` 는 간단한 descriptive statistics 를 확인 할 수 있다.

```py
az.plot_posterior(trace_g);
```

다음은 원문의 그래프 해석이다.

* *Unlike Frequentist inference, in Bayesian inference, we get the entire distribution of the values.*
* *Every time ArviZ computes and reports a HPD, it will use, by default, a value of 94%.*
* *Please note that HPD intervals are not the same as confidence intervals.*
* *Here we can interpret as such that there is 94% probability the belief is between 63.8 euro and 64.4 euro for the mean ticket price.*

첫번째로, *Frequentist inference* 와 *Bayesian inference*의 비교인데 이를 통해서 우리는 아까 *point estimation*과 *Frequentist inference*은 같은 의미를 가진다고 볼 수 있다. 다음으로 주목해야 할 부분은 *[HPD( Highest Posterior Density)](https://en.wikipedia.org/wiki/Credible_interval)* 이다. 원문에서의 *Wikipedia* 링크에서 찾은 설명이다.

* *Bayesian intervals treat their bounds as fixed and the estimated parameter as a random variable, whereas frequentist confidence intervals treat their bounds as random variables and the parameter as a fixed value. Also, Bayesian credible intervals use (and indeed, require) knowledge of the situation-specific prior distribution, while the frequentist confidence intervals do not.*

다시 한번, `Bayesian` 과 `frequentist`를 비교 함을 알 수 있다. 여기서 *Bayesian interval*은 *HPD*를 말하고, *frequentist confidence interval*은 익숙한 그 *Confidence Interval*(CI)을 의미한다. 지금까지 *Bayesian estimation*과 *Point estimation*의 차이를 간략히 하면, 예측하는 parameter 의 값을 변하게 두는지 인 것 같다. *Bayesian estimation*은 예측값이 변할 수 있기 때문에, `prior`라는 예전의 상황과 `posterior`라는 미래의 상황을 모두 담아두기 위해서 각 parameter 마다 분포를 설정하였다. 반면에, *Point estimation*은 예측값이 고정 되어 있고, 대신 그 예측값의 범위가 조정하며 불확실성을 다룸을 알 수 있다.

```py
pm.gelman_rubin(trace_g)
```
Gelman Rubin test의 결과를 본문에서는 다음과 같이 설명한다. *We can verify the convergence of the chains formally using the Gelman Rubin test. Values close to 1.0 mean convergence.* 다시 한번, convergence를 확인해서 우리의 결과의 타당성을 입증했다.

# Model Validation
앞에서 우리는 Gelman Rubin test를 통해서 각 parameter의 값이 converge 함을 확인했다. 다음으로는, 우리가 추론한 값의 타당성을 얻기 위해서, 설정한 모델의 validation을 확인 해야 한다. 원문에서는 *PyMC3*의 `pm_sample_posterior_predictive`을 다음과 같이 설명하였다.
* *Posterior predictive checks (PPCs) are a great way to validate a model. The idea is to generate data from the model using parameters from draws from the posterior.*
* *Now that we have computed the posterior, we are going to illustrate how to use the simulation results to derive predictions.*
* *The following function will randomly draw 1000 samples of parameters from the trace. Then, for each sample, it will draw 25798 random numbers from a normal distribution specified by the values of μ and σ in that sample.*

`pm_sample_posterior_predictive`의 방법을 설명하기 앞서, 우리는 `posterior`의 개념을 이해할 필요가 있다. `posterior`는 앞의 우리가 설정한 `prior`를 기반으로 우리가 가진 데이터를 통해서 예측된 값이다. 따라서 우리가 예측한 `posterior`를 통해서 시뮬레이션을 통해 우리가 예측한 값을 확인하는 방법이다.

```py
ppc = pm.sample_posterior_predictive(trace_g, samples=1000, model=model_g)
np.asarray(ppc['y']).shape
```
> /usr/local/lib/python3.6/dist-packages/pymc3/sampling.py:1247: UserWarning: samples parameter is smaller than nchains times ndraws, some draws and/or chains may not be represented in the returned posterior predictive sample
"samples parameter is smaller than nchains times ndraws, some draws "
100%|██████████| 1000/1000 [00:03<00:00, 329.56it/s]
(1000, 25798)

```py
_, ax = plt.subplots(figsize=(10, 5))
ax.hist([y.mean() for y in ppc['y']], bins=19, alpha=0.5)
ax.axvline(data.price.mean())
ax.set(title='Posterior predictive of the mean', xlabel='mean(x)', ylabel='Frequency');
```


그래프를 보면 실제 우리의 데이터의 기찻값의 평균과 예측된 parameter 들로 표본 평균의 결과가 거의 일치 한다고 볼 수 있다. 따라서 모델의 타당성까지 확인 할 수 있다.

# 결과
우리는 [*Hands On Bayesian Statistics with Python, PyMC3 & ArviZ*](https://towardsdatascience.com/hands-on-bayesian-statistics-with-python-pymc3-arviz-499db9a59501)과 같이 **Bayesian estimation**을 따라해보았다. 기찻값의 표는 Gaussian distribution을 따르며, 평균은 64의 근접 할 것이며, 표준편차는 약 24 정도 될 것이다. 지금까지 **Bayesian estimation** 과정을 따라오면서 **Point estimation**과 비교를 많이 하며 그 특성을 이해 할 수 있었다.

1. 값의 불확실성에 대해서 열려있다. 따라서, 앞으로 변화하는 상황에 대해 좀 더 유동적인 예측임을 알 수 있었다.
2. 모든 예측 파라미터에 대해서 각각의 분포가 존재한다.

아무리 우리가 예측한 값으로 converge 한다는 사실을 알아도, 100% 확신을 할 수는 없다. 어디까지나 모든 통계의 테스트와 같이 통계적 유의미를 갖는 결과일 뿐, 절대적인 진리나 사실이 아님을 다시 한 번 상기하면서 이 글을 마친다.