[이재우] Bayesian Parameter Estimation

[jingwoo4710]

Introduction

Bayesian parameter estimation 이라고 오늘 공부할 주제라고 딱 들으면 누구나 하기 싫어질것 같다. 왠지 많은양의 Bayesian이 쓴 논문을 읽어야 할 것 같고, 외국사람이니 또 영어로 읽을 생각하니 벌써 머리가 아프다. 나도 하기 싫다. 뭐든 공부라고 생각이 들면, 반발심으로 하기 싫다는 말이 먼저 머릿속에 떠오른다. 그래서 나는 Bayesian parameter estimation을 공부해보는 것이 아닌 누군가 엄청 똑똑한 사람이 많은양의 논문을 읽고 생각하고 고민하고 이해해서 쓴글을 그저 따라하면서 효율적으로 이해해보려고 한다. Hands On Bayesian Statistics with Python, PyMC3 & ArviZ 라는 제목으로 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.

이를 통해서 우리가 예측할 것은 기차푯값을 예측함을 알 수 있다.

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를 통해 잘 불러왔는지 확인한다.

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의 값을 처리하였는가이다.

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을 한다면 좀 더 명확한 근거가 있는 방법을 사용해보도록 해야겠다.

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를 설정해도 상관없다.

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)을 내포하고 있다.

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 한다는 사실을 알면 결국 그 끝에는 같은 결과를 마주하게 되기 때문이다.

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

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

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 를 확인 할 수 있다.

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) 이다. 원문에서의 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은 예측값이 고정 되어 있고, 대신 그 예측값의 범위가 조정하며 불확실성을 다룸을 알 수 있다.

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를 통해서 시뮬레이션을 통해 우리가 예측한 값을 확인하는 방법이다.

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)

_, 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과 같이 Bayesian estimation을 따라해보았다. 기찻값의 표는 Gaussian distribution을 따르며, 평균은 64의 근접 할 것이며, 표준편차는 약 24 정도 될 것이다. 지금까지 Bayesian estimation 과정을 따라오면서 Point estimation과 비교를 많이 하며 그 특성을 이해 할 수 있었다.

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

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

[johnnykoo84]

@jingwoo4710 잘 읽었습니다. 궁금한 점은 towardssdatascience 원문 글과 어떤 관계가 있는 것인가요? 1) 번역 2) work flow 및 의견을 참고하여 새로운 의견 제시 3) 기타?

[jingwoo4710]

3) 기타 번역 + 이 글을 통해 배운 Bayesian estimation 의 특징을 소개 하려했습니다.

[johnnykoo84]

@jingwoo4710 원문을 읽어보면 상당히 많은 부분을 인트로에 할애하고 있는 것을 보실 수 있을 겁니다. @jingwoo4710 님의 경우 원문 링크 언급 후, 바로 본론으로 치고 들어가는데요. 이 때 독자는 무슨 일인지도 모르고 멱살 잡고 끌려가는 느낌이 들 순 있습니다. body 를 쓰시는 능력에 비해 재우님은 독자를 꼬시는 인트로 작성에도 신경을 써 주셔도 좋겠습니다.

동시에, 이게 큰 이슈가 될수도, 작은 이슈가 될 수도 있습니다. 글 마지막에 레퍼런스로, 레퍼런스 링크를 넣게 되면, 제가 첫번째 질문 드린 것처럼, 이 글을 자세히 읽는 사람은 의아해 할 수 있습니다. 이게 저자의 글인지, 아니면 번역된 것인지, 번역 + 나의 의견 및 소개를 더한 것인지. 그렇다면 이런 부분을 명확하게 표현하는 것이 좋은 연습입니다. 동시에, 어떤 부분은 해당 원글에서 인용하거나 가져온 것이고, 어떤 부분은 나의 의견이며 나의 주장인 것인지 구분도 필요합니다. 이 글이 외국 영어 글이라서 문제가 될 가능성이 거의 없지만 (이슈가 안 될 가능성이 0은 아닙니다) 그럼에도 불구하고 원글자가 문제를 제기하고 이 글을 사용하신 @jingwoo4710 님은 업계에서 문제가 될 수 있습니다. 예를 들어 이 글을 tensorflow korea 에 올리게 된다면 댓글로 문제를 제기할 수 있겠죠. 100% 번역이라면, 차라리 번역이라고 표현을 해야 하는데 @jingwoo4710 님의 글은 다소 애매한 부분이 있습니다.

@jingwoo4710 님이 앞으로 많은 글들을 쓰시게 될텐데, 이렇게 처음에 이런 마이너한 부분에 대해 명확히 표현법을 짚고 넘어가는 것은 상당히 중요합니다.

[jingwoo4710]

@johnnykoo84 님의 피드백을 수용하여 글을 다시 써보았습니다. 표현하는것에 있어서 약점이라고 생각하는 저에게 디테일한 피드백 정말 감사드립니다. 계속 발전하는 모습 보여주도록 노력하겠습니다.

[johnnykoo84]

@jingwoo4710 다시 읽어보겠습니다. 빠른 반영이 놀랍네요!!! @jingwoo4710 님의 하루 하루 발전 속도가 기대됩니다.

[johnnykoo84]

@jingwoo4710 훨씬 읽기 편하고 따라가기 수월합니다. 박수를 보냅니다!!!!!!