이 포스트는 Hoff, First Course in Bayesian Statistics (2009) 및 학부 통계전산 강의안을 참고하여 정리한 것임을 밝힙니다.

1. 왜 샘플링이 필요할까?

베이지안 통계학을 접하다 보면 빠지지 않는 내용이 MCMC, Gibbs sampling 등 다양한 sampling 방법에 관한 것이다. 기본적으로 통계학의 목적은 모집단으로부터 추출한 샘플을 통해 모집단의 특성을 유추하는 것이므로, 통계학에서 샘플링 방법론은 당연히 중요할 수 밖에 없다. 그런데 유독 베이지안 통계학에서 위와 같은 샘플링 방법들이 자주 등장하는 이유는 무엇일까?

지난 포스트에서 베이지안들의 궁극적인 목적이 posterior distribution임을 살펴본 바 있다. FCB 5장을 보면 normal model에서 mean 과 variance에 대해 joint (posterior) distribution을 구하는 내용이 나온다. 헌데 conjugate prior를 사용해준다 치더라도, $\mu$와 $\sigma^2$의 joint distribution은 우리가(학부생 따위가) 알고 있는 분포도 아니고, 단순히 pdf 만 가지고는 감을 잡기가 어렵다. 이런 경우 $\mu$와 $\sigma^2$의 샘플들을 마구마구 생성한 다음, 이 친구들의 plot을 그려보면 어느정도 분포에 관해 유추해볼 수 있을 것이다. 실제 posterior distribution은 마음처럼 예쁘게 나와주지 않는 경우가 대부분이므로, sampling을 통한 분포 유추는 베이지안들에게 꼭 필요하다고 볼 수 있다.

2. 깁스 샘플링이란 무엇인가?

그럼 이제 깁스 샘플링이 무엇인지 알아보자. 깁스 샘플링의 유용성은 아래와 같다.

Useful when sampling from the joint denstiy $f$ is difficult, but the full conditionals $f_{i|-i}(x_i|x_j,j\neq i)$ are available and sampling from them is easy

즉, 구하고자 하는 것이 joint distribution인데 이것이 어려울 때, full conditional distribution을 활용하여 우회하여 구하자는 것이다.

예를 들어, $\theta$, $\phi$ 라는 두 모수의 joint posterior distribution을 구하고 싶다고 하자.

$$ P(\theta, \phi | y) = P(\phi|y)P(\theta|\phi,y) = P(\theta|y)P(\phi|\theta,y) $$

이므로,

$P(\theta|\phi,y)$ 와 $P(\phi|\theta,y)$라는 full conditional distribution으로부터 $P(\theta, \phi | y)$라는 joint density를 따르는 $(\theta, \phi)$ 를 샘플링할 수 있다. ( joint density를 각각 $\theta$, $\phi$ 의 함수로 보면 $P(\phi|y)$, $P(\theta|y)$ 는 상수로 취급하여 무시할 수 있다. )

따라서 깁스 샘플링의 과정은 다음과 같다.

Step1. Initialize $\theta^{(0)}$, $\phi^{(0)}$ Step2. Sample $\theta^{(i)}$ from $P(\theta|\phi^{(i-1)},y)$ Step3. Sample $\phi^{(i)}$ from $P(\phi|\theta^{(i)},y)$ and so on

• 이를 살펴보면 $(\theta, \phi)$ 샘플은 그 이전 단계의 샘플에만 영향을 받는데 (conditionally independent) , 이는 Markov chain의 특성이다. 따라서 깁스 샘플링은 크게 보면 MCMC에 속한다.

이렇게 얻은 sampling distribution은 샘플 수가 많아지면 target distribution에 수렴한다.

$$ P(\theta^{(s)}, \phi^{(s)} \in A) \rightarrow \int_Ap(\theta, \phi) d\theta d\phi \ as \ s \rightarrow \infin $$

참고) sampling distribution은 sample의 분포가 아니라 statistic의 분포이다. (A sampling distribution is a probability distribution of a statistic obtained from a larger number of samples drawn from a specific population.)

3. 샘플링을 이용한 추정과 표준오차

이렇게 샘플을 생성했다면, 분포를 plotting 하는 것 외에도 다양한 추정을 할 수 있다. 예를 들어 평균에 대해 알고 싶다면 표본평균을 통해 추정할 수 있고, 분산을 알고싶다면 표본분산을 통해 추정할 수 있다.

그런데 그렇게 추정한 estimate가 괜찮은(?) 값인지는 어떻게 결정할 수 있을까? 하나의 기준이 될 수 있는 것이 바로 표준오차이다. 표준오차는 통계량의 분산을 가리키는 개념이다. 우리가 뽑은 샘플은 수많은 가능한 샘플의 경우의 수 중 하나이므로, 샘플을 어떻게 뽑느냐에 따라 추정치는 달라질 것이다. 만약 샘플을 어떻게 뽑든 이 추정치가 크게 달라지지 않는다면, 그 추정치는 나쁘지 않은 추정치라고 볼 수 있을 것이다.

더 자세한 설명과 예시 문제는 다음에...

+) 깁스샘플링과 Metropolis-Hastings 알고리즘의 연관성

깁스샘플링은 Metropolis-Hastings 알고리즘의 speical case로 볼 수 있다. $$ q(y|x) = \frac{1}{d}f_{i|-i}(y|x_j,j\neq i) = \frac{1}{d}\frac{f(y)}{f(x_1,...x_{i-1},x_{i+1},...x_d)} $$

$$ \alpha = min(\frac{f(y)q(x|y)}{f(x)q(y|x)},1) = min(\frac{f(y)f(x)}{f(x)f(y)},1) = 1 $$

따라서 acceptance rate = 1 이므로 생성된 샘플을 reject할 필요가 없는 경우이다.

이 포스트는 Hogg&Craig, Introduction to Mathetmatical Statistics 7th ed.를 학습하고 정리한 것입니다.

1. Bernoulli & Binomial

베르누이 시행: random experiment with two outcomes(ex. success/failure)

let $X(success)=1, X(failure)=0, P(success)=p$, then $$X \sim Bernoulli(p)$$

let $X$=the number of success in $n$ independent Bernoulli trial, then $$X \sim Binomial(n,p)$$

베르누이 분포는 1번의 베르누이 시행에서 성공 여부를 나타내는 것으로 해석할 수 있기 때문에, 독립인 n개의 베르누이 확률변수의 합은 이항분포를 따른다.

추가. Multinomial Distribution

베르누이 분포와 이항 분포의 확장판이다. 한 번의 확률시행에서 나올 수 있는 결과의 경우의 수가 여러 개 일때 사용한다.

2. Geometric & Negative Binomial

Assume Bernoulli trial, let $X$= the number of failure / trial until 1st succes, then $$X \sim Geometric(p)$$

let $X$ = the number of failure / trial until rth success, then $$X \sim NB(r,p)$$

Geometric 분포가 Negative Binomial 분포의 special case인 것은 쉽게 유추할 수 있다. 이 분포들의 또 다른 특징은, 마치 continuous case의 Exponential 분포와 Gamma 분포의 관계와 대응된다는 것이다. 자세한 내용은 추후 공부 해야 하겠지만, 포아송 회귀(Poisson Regression)를 적용하기 위해서는 Poisson 분포의 특징과 같이 평균과 분산이 같다는 가정이 충족되어야 한다. 그러나 실제 data에서는 이 가정이 성립하지 않는 경우가 종종 발생하고, 이때 Negative Binomial 분포를 이용한다고 한다.

추가. Hyper Geometric Distribution

Geometric/Negtive Binomial 분포를 비복원 추출 시에 사용한다면, Hyper Geometric 분포는 복원 추출 시에 사용된다.

let $N$ : the number of elements in population, $n$ : the number of drawn elements without replacement, $r$ : the number of success in population, and $Y$: the number of success in $n$ randomly selected items out of N elements without replacement, then $$Y \sim HyperGeometric(N,n,r)$$

3. Poisson

(1) Poisson Process

fixed interval 에서 사건의 횟수를 생성하는 확률과정이 Poisson 분포를 따른다면, 이 확률과정을 Poisson Process 라고 부른다. Poisson process with rate $\lambda$는 다음과 같은 조건을 만족 해야 한다.

Let $g(x,w)$ : the probability of $x$ events in each interval of length $w$

1. $g(1,h) = \lambda h + o(h)$

   어떤 사건이 아주 짧은 time interval 에서 1번 발생할 확률은 해당 interval의 길이에 비례한다.

   2. $g(0,h) = 1-\lambda h + o(h)$

   3. $\sum_{i=2}^\infin g(x,h)=o(h)$

      어떤 사건이 아주 짧은 time interval 에서 2번 이상 발생할 확률은 거의 0이다. (엄밀한 설명은 아닌 것 같지만 일단 넘어가자.)

   4. "Independet Increments" : The number of events in non-overlapping intervals are independent.

   5. "Stationaty Increments" : The number of events that occur in any interval of time depends only on the length of the time interval.

(2) Poisson Distribution

위와 같은 Poisson Process 에서, let $X$ = the number of events in an interval of length t ,Then

$X\sim Poi(\lambda t)$

4. Exponential & Gamma Distribution

위와 같은 Poisson Process에서 첫번째/$k$번째 사건이 일어날 때 까지의 시간을 나타내는 분포가 exponential 분포와 gamma 분포이다.

The time until 1st event ~ $Exp(\lambda) \equiv Gamma(1,\frac{1}{\lambda})$

! 지수분포 표기는 rate 기준

The time until $k$th event ~ $Gamma(k,\frac{1}{\lambda})$

Exponential 분포가 Gamma 분포의 special case 임은 쉽게 알 수 있다.

5. Normal & Chi-square & t

정규분포는 통계학에서 가장 널리 사용되는 분포 중 하나이다. $\chi^2$와 $t$ 분포는 사실상 검정을 위한 분포라 봐도 무방한데, 기초 통계 방법론의 가설검정에서 정규분포와 chi-square 분포의 관련성, 정규분포와 t분포의 관련성을 종종 이용한다.

$Z\sim N(0,1) \Rightarrow Z^2\sim\chi^2(1)$

$T\sim t(n)$ is similar to $N(0,1)$ when $n$ is very large.

이 포스트는 Hoff, First Course in Bayesian Statistics (2009)를 학습한 내용을 정리한 것임을 밝힙니다.

1. 베이지안 통계학 이란?

전통적으로 통계학은 빈도론적 통계학(Frequentist statistics)을 중심으로 발전해왔다. 내가 배운 학부 수업의 대부분(사실상 베이지안 통계학을 뺀 전부)도 이 빈도론적 통계학을 기반으로 한 내용들이었다. 빈도론적 통계학에서 확률변수 X 는 모수(parameter) $\theta$에 의존하는 확률 분포를 갖고, 모수 $\theta$는 알려져 있지 않은 상수라고 본다.

예를 들어, 동전 하나를 던질 때 앞면이 나올 확률에는 모수에 따른 참값이 존재한다. 하지만 이 값을 정확하게 알 수 없기 때문에, 동전을 몇 번 (많이 던져볼 수록 좋다.) 던져보면서 실험(random experiment)을 하고, 이 실험에서 관찰한 데이터를 토대로 참된 확률 값을 추정하는 방식이다.

반면 베이지안 통계학에서는 확률을 주관적인 믿음이라고 본다. 확률 모수$\Theta$를 고정되어 있는 상수가 아니라 X 처럼 분포를 갖는 확률 변수로 생각하는 것이다. 따라서 모수와 관련된 확률에 참값이란 없으며, 데이터를 이용해 모수에 대해 가지고 있던 사전적인 믿음을 계속해서 업데이트 하는 방식으로 추론한다.

2. 베이지안 추론의 틀

베이지안 추론에서 가장 주요하게 쓰이는 개념은 사전분포(prior distribution), 사후분포(posterior distribution), 그리고 Likelihood 이다. 이 세 가지 개념은 다음과 같이 정의된다.

let $X|\theta \sim f(x|\theta)$ : conditional pdf of X given $\Theta$ = $\theta$ where $X$ = ($X_1, X_2, ... X_n$) : a random sample from $f(x|\theta)$

$\Theta$ ~ h($\theta$) : prior pdf of $\Theta$

$L(x|\theta)= f(x_1|\theta)f(x_2|\theta)...f(x_n|\theta)$ : likelihood

$k(\theta|x) = \frac{g(x,\theta)}{g_1(x)} = \frac{L(x|\theta)h(\theta)}{g_1(x)} \propto L(x|\theta)h(\theta)$ : posterior pdf of $\Theta$

위의 정의를 풀어보면, 어떤 확률변수의 모수 $\Theta$는 분포를 가지는 확률 변수 인데, 이때 $\Theta$의 분포를 prior distribution 이라고 한다. 이 $\Theta$가 어떠한 분포를 가질 것인지, 즉 prior distribution 이 어떠한 모양일 것인지에 대해 가지고 있는 믿음이 Bayesian approach 에서 말하는 확률에 대한 사전적인 믿음이다. 이때 $\Theta$의 분포는 직접적으로 관찰할 수 없기 때문에, 우리는 우리가 관찰한 random sample 을 토대로 추론해야한다. $\Theta$가 확률변수라는 것은 다양한 $\Theta$ 의 값이 존재함을 의미하므로, 우리가 관찰한 random sample $X$는 모수가 어떤 값 $\theta$으로 주어졌을 때 얻은 결과일 것이다. 이러한 conditional joint pdf 를 Likelihood 라고 부른다. 이제 prior distribution 과 likelihood 를 가지고 베이즈 정리를 이용하면, $\Theta|x$ 의 conditional distribution 을 얻는다. 이 분포가 관찰된 데이터($X$)를 토대로 미지의 모수에 대한 믿음을 업데이트한 결과이며, 이를 posterior distribution 이라고 부른다. Bayesian approach 에서는 이 posterior distribution 에 기반하여 추론을 진행한다.

이 포스트는 연세대학교 학부 통계전산을 수강하고 정리한 글 입니다.

통계전산은 무슨 과목일까

학부 통계전산에서는 크게 세가지 정도의 주제를 배웠다. sampling, integration, optimization인데, 학기 초반부터 중후반까지 대부분 여러가지 sampling 기법을 배운 것 같다.

(1) Sampling

R에는 통계학 연구에 사용할 수 있는 다양한 패키지들이 있다. 정규분포 등 간단한 분포로부터의 sampling은 심지어 패키지도 필요없다. 그런데 실제 사용되는 분포는 그렇게 간단할 리가 없고, 이 세상에 존재하는 수많은 분포들을 모두 패키지로 구현해놓는 것은 불가능하기 때문에, 결국 내가 원하는 분포에서 샘플링하는 법을 스스로 익혀두어야 한다. 그럼 R 등의 프로그램이 '어떻게' 원하는 분포로부터의 샘플링을 하도록 할 것인가? 가 주된 물음이다. 이를 위해서 Inverse transform method, Acceptance-Rejection method, Copula를 활용한 dependent variables 의 생성, MCMC와 Gibbs sampling 등을 배운다.

추가적으로, '좋은 ' 통계적 추정을 위해서는 표준오차(통계량의 variance)가 작은 것이 유리하다. 이와 관련하여 몇 가지 variance reduction technique 도 다룬다.

이 과목을 수강할 당시에는 베이지안 통계학에 대한 개념이 전무해서 몰랐는데, 관련 내용이 베이지안 통계학에서 아주 유용하게 사용될 수 있다고 한다. 당장 베이즈 기본서를 펼쳐보면 MCMC와 Gibbs 샘플링이 최소 10번 등장하는 것을 볼 수 있다...

(2) Optimization

학부 수리통계학을 배우면서 가장 많이 접하는 개념 중 하나가 MLE(Maximum Likelihood Estimator)이다. MLE를 구하는 방법을 익히기 위해서 온갖 문제를 풀며 손목을 혹사하지만, 사실상 손계산으로 MLE를 구할 수 있는 경우는 얼마 되지 않는다. 예를 들어, 감마분포 parameter의 mle는 어떻게 구할 수 있을까?

기본적으로 MLE는 likelihood ftn의 argmax 구하기(머신러닝에서는 이에 -를 붙여 loss function의 argmin을 구하기도 하는 것 같다), 즉 최적화(optimization) 문제와 같다. 따라서 손계산이 불가능한 경우 numerical optimization 방법을 이용해야한다. 이와 관련한 대표적인 방법으로는 Newton's method, Gradient Descent 등이 있을 것이다.

(3) Integration

베이지안 통계학에서를 예로 들어 보자. Bayes theorem을 보면, 분모 부분은 marginal density 이기 때문에 필연적으로 적분 계산을 맞닥뜨리게 된다. 그러나 학부에서 배우는 간단한 분포도 손적분이 어려운데 실제로 필요한 분포를 손적분할 수 있을리가 만무하다. 앞서 sampling에서는 필요한 분포의 sample을 많이 뽑아 적분값을 근사시켰다면, 여기서는 적분 자체를 근사하고자 한다. (ex. 적분이 불가능한 함수를 비슷한 함수로 대체) 사실 둘이 맥락이 비슷하기 때문에 완전히 분리하여 설명하는 것은 어려울 듯하다. 어찌되었든 적분 계산 방법에 대한 논의는 통계전산의 주요한 토픽 중 하나이다. 보통 convex 함수의 미분을 구할 때 일차 미분식이 0이 되도록 하는 argument를 구한다는 점을 생각해보면, optimization 문제도 여기에 포함될 수 있을 것 같다.

universal approximation 이란?

"A neural network with one hidden layer can approximate any function $f:R^N\to\R^M$ with arbitrary precision."

1. 개념을 직관적으로 이해하기

미적분 시간에 리만 적분의 idea 를 그림으로 나타낸 것을 본적이 있다. 그래프 아래의 면적을 근사하기 위하여 주어진 함수의 구간을 여러개의 아주 작은 직사각형(깍두기!) 으로 나누어, 이를 무한히 더하는 개념이다.

universal approximation의 개념도 이와 비슷한게, 여러개의 (아주 작은 구간에서 정의된) 함수의 합으로 나타내고자 하는 것이다. 다만 그 대상이 나타내고자 하는 것이 함수 아래의 면적이 아닌 함수 $g$ 이다.

이 아이디어를 기반으로, neural network로 $g$를 나타낼 수 있다는 것이 universal approximation 이다. 다음과 같은 함수 $g$를 표현하고자 한다고 가정하자.

2. Perceptron으로 만들기

기본적인 perceptron은 왼쪽과 같이 생겼는데, 이를 0.1만큼 평행이동 하여 새로운 perceptron을 만들고, 이를 원래의 perceptron 에서 (-) 하면 [0,0.1] 범위에 '깍두기' 를 만들 수 있다. 이걸 표현하고자 하는 함수에 적용해보자.

• 0과 0.1 구간의 2.7이라는 값은 다음과 같이 표현할 수 있다: $2.7(f(x)-f(x-0.1))$

• 0.2와 0.3 구간의 3.3이라는 값은 다음과 같이 표현할 수 있다: $3.3(f(x-0.3)-f(x-0.2))$

What if 아주 촘촘하게 깍두기들을 연결하면? -> 함수가 표현 된다!!!

이를 neural net으로 표현해보자.

간단한 neural net으로 g를 표현하였다! 실제로 근사의 정확도를 높이기 위해서는 엄청 나게 많은 node 가 필요할 것이다.

3. Sigmoid로 만들기

기본적인 sigmoid를 다음과 같이 변형하면, 거의 perceptron과 유사한 모양이 된다. 따라서 (2)에서와 같은 방법으로 $g$를 근사할 수 있다.

4. ReLU로 만들기

will be updated

[References]

• Youtube "혁펜하임"

• Hoff, A First Course in Bayesian Statistical Methods (2009)
• Doing Bayesian Data Analysis