원문 : http://jakevdp.github.io/blog/2014/03/11/frequentism-and-bayesianism-a-practical-intro/
이 포스팅은 5파트로 나눠 포스팅하겠다.
데이터사이언티스트가 통계학에 대해서 처음으로 들었을 것들중 하나는, 통계학에는 2가지 다른 접근인 빈도주의와 베이지안주의가 있다는것이다. 이것들의 중요성에도 불구하고, 많은 연구자들은 둘 사이의 특징들과 결과에 대한 다른 실용적 접근법들에 대해 배울 기회가 없었다. 이 포스팅의 목적은 빈도주의와 베이지안주의의 접근법들에 대한 철학적인 면과 실용적인 면을 함께 보여, 데이터 분석의 종류들을 이해하는데 무리가 없도록 하는것이다.
Frequentism vs. Bayesianism : Philosophical Debate
기본적으로, 빈도론자들과 베이지안들 사이에서의 논쟁은 확률의 정의에 대한 것이다.
빈도주의자들에게, 확률은 반복된 측정값들 중 하나의 일부케이스의 측면이라는 관점에서 의미를 가진다. 즉 만약 (별의 유속이 시간에 영향을 받지 않는다는 가정하에) 한 별로부터 광량자속(photon flux) F를 측정한다면, 이것을 측정하고 또 측정하여 매 시행마다 측정 도구의 통계적 오차 때문에(측정오차) 약간씩 다른 답을 얻을 것이다. 많은 수의 측정값의 한계를 가진 상태에서, 주어진 어느 값의 빈도는 그 값을 측정할 확률을 가리킨다. 빈도론자들에게, 확률들은 기본적으로 확률들의 빈도들에 관계가 있다. 예를들어 엄격한 빈도주의론자의 입장에서 보면, 별의 진정한(true) flux의 확률에 대해 이야기하는것은 의미가 없다는 것을 가리킨다. 왜냐하면 true flux는 (정의에 따르면) 단 하나의 고정된 값이기 때문에 고정된값에 대한 한 빈도분포에 대해 이야기 하는것은 이치에 맞지 않다는 것이다.
베이지안주의자들에게, 확률의 개념은 상태(statements)에 대한 확신(certainty)의 정도를 다루는 것으로 확장시켰다. 한 베이지안이 어떤 확률로 한 별의 flux F 를 측정한다고 하자 P(F) , 이 확률은 수많은 반복된 실험을 한다는 제약에서 나온 빈도들로부터 확실하게 측정될 수 있다. 그러나 이것은 중요한 것이 아니다. 확률은 측정값이 무엇이 될지에 대한 나의 앎의 상태(statement of my knowledge)이다. 베이지안들에게 확률은 한 사건에 대한 우리의 지식(앎)과 기본적으로 관련이 있다. 예를 들어 베이지안의 입장에서 보면, 우리는 별의 참된(true) flux가 주어진 범위에 있을 확률에 대해 의미있는 말를 할 수 있다. 그 확률은 사전정보 그리고/혹은 이용가능한 데이터에 기반한 값에 대한 우리의 지식/믿음을 나타낸다.
놀라운것은 이 아주 미묘한 철학적 차이가 실제로 통계적 데이터분석을 할때 상당히 다른 방향의 접근을 하게 한다.
밑에 이에 대한 예를 들어보겠다.
Frequentist and Bayesian Approaches in Practice : Counting Photons
수학적 공식이 약간 들어가겠지만 깊이 들어가지는 않을 것이다. 만약 더 알고 싶다면, textbook의 4-5 챕터를 보길 바란다.
The problem : Simple Photon Counts
망원경으로 하늘에 있는 한 별로부터 온 불빛을 관측한다고 해보자. 당분간, 별의 true flux는 시간에 따라 바뀌지 않는다고 가정하자. 다시말해 고정된 값 
측정오차에 대해 불필요한 부분을 제거한다. 즉 우리는 오차가 가우시안 분포를 따른다는 가정을 한다. 빈도주의론적 관점에서,
아래 python으로 작성한 예제가 있다. 포아송분포가 측정과정에 좋은추정이다.
In [1]:
# Generating some simple photon count data
import numpy as np
from scipy import stats
np.random.seed(1) # for repeatability
F_true = 1000 # true flux, say number of photons measured in 1 second
N = 50 # number of measurements
F = stats.poisson(F_true).rvs(N) # N measurements of the flux
e = np.sqrt(F) # errors on Poisson counts estimated via square root이제 측정된 데이터를 간단하게 시각화해보자.
%matplotlib inline
import matplotlib.pyplot as plt
fig, ax = plt.subplots()
ax.errorbar(F, np.arange(N), xerr=e, fmt='ok', ecolor='gray', alpha=0.5)
ax.vlines([F_true], 0, N, linewidth=5, alpha=0.2)
ax.set_xlabel("Flux");ax.set_ylabel("measurement number");
이 각 측정들은 포아송 분포로 부터 추정된 다른 오차값을 가진다. 우리는 이미 true flux
Frequentist Approach to Simple Photon Counts
우리는 전통적 빈도주의론자가 maximum likelihood approach를 한다고 하자. 한 관찰치
이 식은 "
우리는 각 데이터 포인트에 대한 확률들의 곱(결합확률)을 계산하여 likelihood function을 만들 수 있다.
여기서, 
우리가 하고 싶은 것은, 우도가 최대가 되는 
모든 오차들
즉, 직관에 따라, 
이러한 결과들을 단순한 계산들로 얻을 수 있다. 작은 예제데이터셋으로 평가를 해보자.
w = 1. / e ** 2
print("""
F_true = {0}
F_est = {1:.0f} +/- {2:.0f} (based on {3} measurements)
""".format(F_true, (w * F).sum() / w.sum(), w.sum() ** -0.5, N))
flux 50개를 예측한다. 결과는 우리의 예측치는 약 0.4%의 오차를 가지고 입력값과 일치한다.
Bayesian Approach to Simple Photon counts
당신이 기대했듯, 베이지안 접근은 확률로 시작해서 확률로 끝난다. 이것은 기본적으로 우리가 계산하고자 하는것이 질문, 즉 이 사건에서 모수들에 대한 우리의 지식이라는것을 알려준다.
이 문제에 대한 공식이 근본적으로 확률은 
이 결과를 계산하기 위해, 베이지안들은 다음으로 기본적인 확률법칙인 베이즈 이론(Bayes'Theorem)을 적용한다.
베이지안 이론이 베이지안이라는 그들의 이름을 얻은 곳이지만, 베이지안 이론은 이 법칙자체가 아닌, 
이 수식에 나온 각각의 용어들을 보자.
는 우도(likelihood)이다. 이는 위의 빈도론자의 접근법에서,
에 대응한다.
는 모델 사전확률이다. 이는 데이터 D의 적용이전에, 모델에 대해 우리가 아는 것이다.
는 데이터 확률, 실제로 이는 단지 normalization term 에 해당된다.
만약 우리가 사전확률 
But What About the Prior?
여러분도 알아챘겠지만, 위에서 뭔가를 대충 얼버무리고 갔다. 바로 사전확률
한 빈도론자는 사전확률은 참사전정보(true prior information)을 알 수 없을때, 사전확률은 문제가 된다고, 지적할 것이다. 위에서 본 flat prior같이 noninformative prior를 사용하는것이 직관적인 판단으로 보이나, 여기에는 아주 약간의 미묘한 문제들이 있다. 이것은 많은 상황에서, 진정한 noninformative prior는 존재하지 않는다는 것으로 밝혀졌다. 빈도론자들은 결과에 반드시 편향을 일으키는 사전확률의 주관적 선택은 통계분석에서 설자리가 없다는 것을 지적했다.
한 베이지안은 빈도주의가 이 문제를 풀수 없고 단지 문제를 회피한다고 맞받아쳤다. 빈도주의는 사전확률의 몇몇 (암묵적) 선택에 대한 베이지안 접근법의 특수한 케이스로 볼 수 있다. 즉 다시말해, 베이지안은 비록 사전확률에 대한 선택이 약간의 주관성을 포함할지라도, 암묵적인 선택을 형시적으로 만드는 것이 더 낫다고 말한다.
Photon Counts : the Bayesian approach
이 확률에 대한 가치관의 문제를 잠깐 접어놓고, 어떻게 베이지안 결과가 일반적으로 나오는지 실제로 알아보자. 여기서 다뤘던 one parameter problem 에 대해서,
여기서 MCMC에 대한 이론에 대해 자세히 다루지는 않을 것이다. 대신 emcee package를 사용하여 MCMC 접근법을 적용하는 예제를 보이겠다. 사후 확률 분포로부터 뽑힌 한 셋의 데이터들을 생성하여 우리가 찾는 답을 얻기 위해, 이렇게 뽑힌 데이터들을 사용하는 것이 목표임을 명심해라.
MCMC를 실행하기 위해서, 적절하게 normalize할 필요가 없는 , 사전확률 
def log_prior(theta):
return 1 # flat prior
def log_likelihood(theta, F, e):
return -0.5 * np.sum(np.log(2 * np.pi * e ** 2)
+ (F - theta[0]) ** 2 / e ** 2)
def log_posterior(theta, F, e):
return log_prior(theta) + log_likelihood(theta, F, e)
연속적인 데이터들에 대한 임의시작추측(random starting guesses)을 생성하는것을 포함해, 이제 우리는 문제를 설정했다.
ndim = 1 # number of parameters in the model
nwalkers = 50 # number of MCMC walkers
nburn = 1000 # "burn-in" period to let chains stabilize
nsteps = 2000 # number of MCMC steps to take
# we'll start at random locations between 0 and 2000
starting_guesses = 2000 * np.random.rand(nwalkers, ndim)
import emcee
sampler = emcee.EnsembleSampler(nwalkers, ndim, log_posterior, args=[F, e])
sampler.run_mcmc(starting_guesses, nsteps)
sample = sampler.chain # shape = (nwalkers, nsteps, ndim)
sample = sampler.chain[:, nburn:, :].ravel() # discard burn-in points만약 이것이 모두 올바르게 작동한다면, 배열 sample은 사후 확률(분포)로부터 뽑힌 50000개의 연속데이터를 담고 있어야 한다. 그래프를 그래 확인해 보자.
# plot a histogram of the sample
plt.hist(sample, bins=50, histtype="stepfilled", alpha=0.3, normed=True)
# plot a best-fit Gaussian
F_fit = np.linspace(975, 1025)
pdf = stats.norm(np.mean(sample), np.std(sample)).pdf(F_fit)
plt.plot(F_fit, pdf, '-k')
plt.xlabel("F"); plt.ylabel("P(F)")
결국 (정규) 사후 분포로 부터 한 표본 데이터들을 뽑은 것이다. 이 사후확률분포의 평균과 표준편차들은 위의 빈도론자의 최대우도추정(Maximum likelihood estimate)결과와 당연히 같다.
In [7]:
print("""
F_true = {0}
F_est = {1:.0f} +/- {2:.0f} (based on {3} measurements)
""".format(F_true, np.mean(sample), np.std(sample), N))
우리가 예측 했듯이,, 이 단순한 문제에 대해 베이지안 접근은 빈도주의자의 접근과 같은 결과를 냈음을 확인했다.
Discussion
이제 , 베이지안 방법이 불필요하게 복잡하다는 인상을 가지게 되었는 지도 모른다. 이 상황에서는 확실이 그렇다. Affine Invariant Markov Chain Monte Carlo Ensemble sampler를 일차원 정규분포를 만들기 위해 사용하는 것은 비치발리볼을 부수기위해 Death star(혹시 모르는 사람을 위해 스타워즈에 비유한것)를 사용하는 것과 마찬가지이다. 그럼에도 나는 여기서 이것을 사용했다. 왜냐하면 이것은 다차원들에서 복잡한 사후확률들을 다루고, 분석적 우도 접근법이 불가능한, 더 복잡한 상황들에서 좋은 결과를 낼수 있는 하나의 접근법을 보여주기 때문이다.
첨언을 하자면, 여러분은 한가지 교묘한 속임수를 눈치 챘을 것이다. 끝에서, 우리는 사후확률분포 표본들을 만드는데 빈도주의 접근법을 사용했다! 위에서 표본 평균과 표준편차를 계산했을때, 사후 확률분포를 얻기 위해, 확실히 빈도주의 기술을 사용했다. 순수 베이지안 결과는사후확률분포 있는 그대로 내놓는 것이다(즉 이것의 대표적 표본). 다시말해 순수 베이지안주의에서, 한 문제에 대한 정답은 error bar를 가진 단 하나의 값이 아니다. 즉 정답은 모델 파라미터들에 대한 사후확률 분포이다.
Adding a Dimension : Exploring a more sophisticated model
간단하게 약간 더 복잡한 상황을 봐보고 빈도론자들과 베이지안주의 결과들을 다시한번 비교해 보자. 위에서 우리는 별이 정적이라고 가정했다. 이제 우리는 약간의 확률적 변동을 가진 한 객체가진다고 가정하자. 다시말해 시간에 의해 변동이 된다,
우리는 단순히 모수가 2개인 가우시안 모델을 제안한다. 
이제, 우리는 다시 N개의 관찰치들을 해당 오차와 함께 고려할 것이다. 우리는 관찰치들을 아래 방식처럼 생성할 수 있다.
np.random.seed(42) # for reproducibility
N = 100 # we'll use more samples for the more complicated model
mu_true, sigma_true = 1000, 15 # stochastic flux model
F_true = stats.norm(mu_true, sigma_true).rvs(N) # (unknown) true flux
F = stats.poisson(F_true).rvs() # observed flux: true flux plus Poisson errors.
e = np.sqrt(F) # root-N error, as aboveVarying Photon Counts : The Frequentist Approach
결과 우도는 오차 분포를 가진 내부적 분포의 convolution 이다.
위와 유사하게, 우리는 평균에 대한 최고의 예측치를 찾기위해, 이 우도를 최대화할 수 있다. (Maximum likelihood estimation)
여기에 한가지 문제가 있다. 평균의 최적값은 표준편차의 최적값에 의존한다는 것이다. 결과들이 연관되어 있기 때문에, 빈도주의자들의 결과에 도달 하기 위해 더 이상 직관적인 분석 방법을 쓸 수 없다.
그럼에도 불구하고, 우리는 최대 우도값을 결정하기 위해, 수적 최적화기법들을 사용할 수 있다. 여기에서, Scipy의 optimize submodule에 있는 최적화 routines를 사용할 것이다.
def log_likelihood(theta, F, e):
return -0.5 * np.sum(np.log(2 * np.pi * (theta[1] ** 2 + e ** 2))
+ (F - theta[0]) ** 2 / (theta[1] ** 2 + e ** 2))
# maximize likelihood <--> minimize negative likelihood
def neg_log_likelihood(theta, F, e):
return -log_likelihood(theta, F, e)
from scipy import optimize
theta_guess = [900, 5]
theta_est = optimize.fmin(neg_log_likelihood, theta_guess, args=(F, e))
print("""
Maximum likelihood estimate for {0} data points:
mu={theta[0]:.0f}, sigma={theta[1]:.0f}
""".format(N, theta=theta_est))
이 최대 우도값은, 우리 모델에 중요한 부분인, 평균과 표준편차 모수에 대한 최고의 추정치를 준다. 그러나 이것은 오직 반쪽짜리 답이다. 우리는 이 답에 대해 얼마나 신뢰하는지를 결정할 필요가 있다. 즉 우리는 평균과 표준편차에 대한 error bar를 계산할 필요가 있다.
빈도주의 패러다임에서 오차들을 결정하는데 몇몇 접근들이 있다. 위에서와 같이, 우리는 최대 우도에 대해 일반 추정을 적합할수 있고 공분산 행렬을 내놨다. ( 여기서, 우리는 분석적이라기 보다는 수치적으로 적합해야만 했다.) 대안적으로, 우리는 
이 모든것은 유용한 기술들일것이나, 각 기술들은 그 기술만의 가정과 기술간 미묘한 차이들이 있다. 여기서, 간략하게, 우리는 astroML 패키지에서 찾은 기본 부트스트랩 resampler을 사용할 것이다.
from astroML.resample import bootstrap
def fit_samples(sample):
# sample is an array of size [n_bootstraps, n_samples]
# compute the maximum likelihood for each bootstrap.
return np.array([optimize.fmin(neg_log_likelihood, theta_guess,
args=(F, np.sqrt(F)), disp=0)
for F in sample])
samples = bootstrap(F, 1000, fit_samples) # 1000 bootstrap resamplings위의 MCMC 베이지안 사후확률에 대한 방식과 유사한 방식으로, 우리는 모수들에 대한 오차들을 결정하는 표본평균과 표준편차를 계산할 것이다.
mu_samp = samples[:, 0]
sig_samp = abs(samples[:, 1])
print " mu = {0:.0f} +/- {1:.0f}".format(mu_samp.mean(), mu_samp.std())
print " sigma = {0:.0f} +/- {1:.0f}".format(sig_samp.mean(), sig_samp.std())
부스트램 샘플링의 디테일한 것들에서 큰 맥락은 각 모델 파라미터의 평균과 표준편차를 단순히 찾는 것으로 나오지 않는 non-Gaussian 또는 상관관계가 있는 오차들이 잠재한다는 것이다. 그럼에도 불구하고 이 방법이 이 문제에 빈도주의 접근법의 기본적인 생각을 준다는 것을 믿는다.
Varying Photon Counts : The Bayesian Approach
이 문제에 베이지안 접근법은 , 이전 문제에서의 접근법과 거의 같다. 일단 위의 코드를 약간 바꿔 설정을 해보자.
def log_prior(theta):
# sigma needs to be positive.
if theta[1] <= 0:
return -np.inf
else:
return 0
def log_posterior(theta, F, e):
return log_prior(theta) + log_likelihood(theta, F, e)
# same setup as above:
ndim, nwalkers = 2, 50
nsteps, nburn = 2000, 1000
starting_guesses = np.random.rand(nwalkers, ndim)
starting_guesses[:, 0] *= 2000 # start mu between 0 and 2000
starting_guesses[:, 1] *= 20 # start sigma between 0 and 20
sampler = emcee.EnsembleSampler(nwalkers, ndim, log_posterior, args=[F, e])
sampler.run_mcmc(starting_guesses, nsteps)
sample = sampler.chain # shape = (nwalkers, nsteps, ndim)
sample = sampler.chain[:, nburn:, :].reshape(-1, 2)
이제 우리는 표본들을 가졌다. 우리는 1 그리고 2표준편차를 나타내는 trace들과 contours을 그리기 위해 astroML 로부터 편리한 routine을 사용할것이다.
from astroML.plotting import plot_mcmc
fig = plt.figure()
ax = plot_mcmc(sample.T, fig=fig, labels=[r'$\mu$', r'$\sigma$'], colors='k')
ax[0].plot(sample[:, 0], sample[:, 1], ',k', alpha=0.1)
ax[0].plot([mu_true], [sigma_true], 'o', color='red', ms=10);
빨간점은 우리의 문제설정에서 부터 있는 기본참값을 나타낸다. 그리고 contour들은 68% 그리고 95% 신뢰구간을 가진 표준편차들을 나타낸다. 다른말로, 이 분석에 기반해, 우리는 안쪽 등고선(contour)안에 모델이 들어간 68% 신뢰수준 그리고 바깥쪽 등고선 안쪽 모델에 있는 95% 신뢰구간안에 위치를 한다,
주의할 것은 2 표준편차 안에서, 표준편차 
다른 주의할 점은, 이 사후확률분포가 확실히 가우시안분포는 아니라는 것이다. 즉 수직 방향에서, 비대칭성이 약간 보일 수 있다. 그것은 빈도주의에서 쓰인 가우시안 추정이 결론적으로 참 불확실성을 반영하지 못한다는 것이다. 이것은 빈도주의자체와는 문제가 되지 않는다. ( 즉 빈도주의론에서 non-Gaussianity를 다루는 확실한 방법들이 있다.) 그러나 대다수의 응용 빈도주의 기술들이 정규분포의 암묵적 혹은 명시적 가정을 한다. 베이지안 접근들은 이러한 가정들을 일반적으로 요구 하지 않는다.
(사전확률에 대한 부가 설명 : 
There are interesting arguments to be made that the Jeffreys Prior would be more applicable. Here I believe the Jeffreys prior is not suitable, because is not a true scale factor (i.e. the Gaussian has contributions from as well). On this question, I'll have to defer to others who have more expertise. Note that subtle — some would say subjective — questions like this are among the features of Bayesian analysis that frequentists take issue with).
Conclusion
이 포스팅을 통해 얼마나 빈도주의와 베이지안주의의 기저에 깔린 철학적 차이들이 간단한 문제들에 근본적으로 다른 접근들을 하고 있는지를 보였다. 그럼에도 불구하고 이 둘은 매우 비슷한 혹은 똑같은 결과를 낼 수 있다.
차이점들을 요약하자면
- 빈도주의는 실제 혹은 가정된 사건들의 빈도들과 연관된 확률들을 고려한다.
- 베이지안주의는 앎의 정도를 측정한 것을 확률로 간주한다.
- 빈도주의는 일반적으로 점추정과 최우도 접근법의 사용을 통해 분석을 진행한다.
- 베이지안 분석은 몇몇 버전의 MCMC 샘플링을 통해 혹은 직접 사후확률을 계산한다.
'ML & DM > Bayesian ML' 카테고리의 다른 글
| Frequentism and Bayesianism II: When Results Differ (0) | 2016.05.12 |
|---|---|
| Bayesian machine learning (2) | 2016.03.30 |
