exestudiary

개요

Population(모집단)은 Population의 parameter들을 통해 최소한은 설명할 수 있다. 

Population의 parameter들은 mean, proportion, variance 등을 말하는데, 여기서 문제가 있다.

Population은 매우 많고 크기 때문에 조사하기에는 어려움이 있다. 

즉, Population의 실제 parameter들을 알기에는 다소 어려운게 있다는 뜻이다. 

방법

여기서 Statistics 또는 point estimator들은 Population parameter들을 estimate(추정)하는데 사용된다. 

Estimator는 sample의 함수, 즉, sample을 통해 parameter의 estimate(추정치)를 계산하는 함수라고 생각하면 된다. 

그럼 여기서 말하는 estimate는 뭘까?

Estimate는 sample을 통해 계산된 estimator의 값을 의미한다. 

예시

회귀분석에서 구하는 $\hat{\beta}_{OLS} = (X^TX)^{-1}X^TY$는 그럼 estimator일까? estimate일까?

Population의 parameter를 추정하기 위해 데이터로 정의된 수식이다.

데이터(sample)를 통해 Population의 parameter의 estimater를 계산하기 위한 함수라고 볼 수 있기 때문에 Estimator.

그럼 estimate는 이러한 $\hat{\beta}_{OLS}$를 통해 얻은 값이기에 만약에 관측 데이터 $X,Y$가 주어졌을 때, $\hat{\beta} = 2.5$라고 하면 2.5가 estimate가 된다. 

Frequentist vs Bayesian

이렇게 estimator와 estimate의 차이를 알았다면 Frequentist 관점에서 볼때 뭐가 Fixed value이고 뭐가 random variable인지, Bayesian 관점에서 볼때 뭐가 fixed value이고 뭐가 random variable인지 알아보자.

Frequentist 

Population의 parameter는 항상 고정되어 있지만 우리는 모르는 값이다. 

회귀에서 그렇다면 parameter들($\beta_0, \beta_1, \cdots$)은 항상 fixed되어 있지만 우리는 모르는 값이다. 

Estimator는 데이터에 따라 달라진다. 그렇기 때문에 Random variable로 취급된다. 

회귀에서 그렇다면 $\hat{\beta} = (X^TX)^{-1}X^TY$에서 $Y$가 random이기에 $\hat{\beta}$도 random이다. 

Bayesian

Bayesian에서는 parameter를 random variable로 취급을 한다.

그렇기에 prior를 부여해서 $\beta$가 가질 수 있는 값의 uncertainty(불확실성)을 모델링한다. 

그렇기 때문에 point estimate보다는 posterior distribution 자제가 관심대상이 되는 것이다. 

pdf

Gibbs sampler

• The Gibbs sampler is another special case of the Metropolis-Hastings sampler.
• The Gibbs sampler is often applied when the target is a multivariate distribution, i.e. we are interested in generating sample from $$f(x)  = f(x_1, x_2, ..., x_d)$$
• Let $X  = (X_1, ... , X_d)$ be a d-demensional random vector.
• Define $$X_{(-j)} = (X_1, .... X_{(j-1)}, X_{(j+1)}, ..., X_d)$$
• Denote the conditional distribution of $X_j$ given $X_{(-j)}$ by $f(X_j|X_{(-j)})$
• The Gibbs sampler can be performed by using the conditional distributions $f(X_j|X_{(-j)}), j= 1, 2. ..,d$

Algorithm

• The Gibbs sampler is performed by the following algorithm.
  1. Set the initial value $X^{*} = (X_1^{*}, X_2^{*}, ..., X_d^{*})$ and t = 1.
  2. Repeat for $j = 1, 2....,d$;
     • Update $X_j^* \sim f(X_j|X_{(-j)}^*)$
  3. Set $X^{(t)} = X^{*}$
  4. Set $t \leftarrow t+1$ and go to Step 2.

Example: Gibbs sampler

Bivariate Normal Distribution

• Our goal is to generate samples from $$X = (X_1, X_2) \sim N_2\{\begin{pmatrix} \mu_1 \\ \mu_2\end{pmatrix},\begin{pmatrix} \sum_{11} &\sum_{12} \\ \sum_{21} & \sum_{22} \end{pmatrix}\}$$.
• It is well known that $$X_1 | X_2 \sim N(\mu_1 + \sum_{12} {\sum_{22}}^{-1}(X_2 - \mu_2), \sum_{11} - \sum_{12}{\sum_{22}}^{-1}\sum_{21})$$. $$ X_2 | X_1 \sim N(\mu_2 + \sum_{21} {\sum_{11}}^{-1}(X_1 - \mu_1), \sum_{22} - \sum_{21}{\sum_{11}}^{-1}\sum_{12}) $$.
• Hence, the Gibbs sampler is to generate samples from these conditional distributions iteratively.

R code 

# Gibbs sampler 
mu = c(1, 2) 
Sigma = matrix(c(1, 0.5, 0.5, 1), 2, 2) #target
x = c(0, 0) #initial value 
m = 10000 #length of chain
X = matrix(0, m, 2)

for(i in 1:m){
	mu1 = mu[1] + Sigma[1,2] / Sigma[2,2] * (x[2] - mu[2])
	var1 = Sigma[1,1] - Sigma[1,2]/Sigma[2,2]*Sigma[2,1]
	x[1] = rnorm(1, mu1, sqrt(var1))

	mu2 = mu[2] + Sigma[2,1] / Sigma[1,1] * (x[1] - mu[1])
	var1 = Sigma[1,1] - Sigma[1,2]/Sigma[2,2]*Sigma[2,1]
	x[2] = rnorm(1, mu2, sqrt(var2))
	
	X[i,] = x
}

var(X)
apply(X,2,mean)

var(X[-(1:3000),]) 
# 충분히 iteration이 끝난 뒤에 Bivariate normal을 따른다고 할 수 있음. 
# 그래서 Burn-in period를 잘라줘야함. 
apply(X[-(1:3000), ],2, mean)

pdf

이번에 Markov Chain Monte Carlo(MCMC)를 배우면서 계속 생각했는데

도대체 왜 Bayesian에서는 MCMC를 그렇게 강조하는 것인가?

Bayesian Computation with R의 5장을 공부하면서도 그냥 막연하게 Rejection Sampling 배우네.. 엇 Importance Sampling 배우네.. 엇 SIR(Sampling Importance Resampling) 배우네... 하고 넘겼는데 

6장 Markov Chain Monte Carlo Methods를 배우면서 왜 MCMC를 그렇게 강조했는지 알겠다!

Bayesian Statistics에서는 $$p(\theta | y) = \frac{f(y|\theta)p(\theta)}{\int f(y|\theta)p(\theta) d\theta}$$라는 수식을 통해 $y$를 update하는 개념을 사용해 데이터가 반영된 parameter distribution을 추정하는 것을 목표로 한다. 

그런데 $\int f(y|\theta)p(\theta) d\theta$ 부분의 경우 추정 대상 parameter가 늘어남에 따라 막대한 계산량이 동반된다. 여기서 MCMC 기법이 사용되는 것이다. 

그래서 MCMC(Markov Chain Monte Carlo)가 뭐냐면 Markov Chain property를 사용하여 simulation을 통해 $\int f(y|\theta)p(\theta) d\theta$ 부분을 계산하지 않아도 $p(\theta|y)$를 찾아낼 수 있는 방법을 말하는 것인데, 

Markov Chain property에 따르면 Transition Probability Matrix가 Detailed Balance Equation을 만족한다면 n이 충분히 클 때 특정한 한 분포로 수렴하는데 Transition Probability Matrix를 Bayes' Thm의 Likelihood * Prior 부분(Density)으로 설정하면 Detaild Balance Equation(A Markov Chain is said to be reversible if there exists a probability distribution on $\pi$ on $\Omega$ such that$P(x | y)\pi(y) = P(y | x)\pi(x)$)을 만족시키면서 Stationary Distribution이 $P(\theta|y)$이 확률 과정 문제로 만들 수 있다. 이 때 확률과정의 샘플링을 충분히 한다면 샘플들의 분포는 모수 분포를 따르게 될 것이다. 

그래서!!

직접 샘플링은 어렵지만, Markov Chain을 잘 설계해서 결국 우리가 원하는 분포에 도달하도록 만들자.

라는 것인데, 

즉, 어떤 초기 state에서 출발하여 

Transition Matrix와 같은 rule에 입각하여 따라가면

Posterior과 비슷한 Stationary distribution을 따르는 Markov Chain을 만들 수 있다.

그래서 이 Stationary distribution에서 Sampling을 하면 이게 곧!! Posterior Sample이 되는 것이다. 

아하..

이제 왜 배우는지 알고 나니깐 더 열심히 해야겠다는 생각이 든다. 

방학 때 꼭 Bayesian Computation with R 끝내보자...!!!!

Let X be a random variable with pdf (or pmf) f(x). 

Our goal is to generate a sample of X, but it is impossible to draw a sample from $f(x)$ directly.

Let Y be a random variable with pdf(or pmf) $g(y)$ such that $f(t) / g(t) \leq c$ for some constant c.  

The AR method proceeds as follows:

1. Generate y from the distribution with g.
2. Generate u from Uniform(0,1). 
3. If $u < f(y)/{cg(y)}$, then accept y. i.e., x =y;
   otherwise reject y go to step1. 

Note that g and c are determined by us. 

(We call g as proposal distribution.)

Ex

Let X ~ Beta(2,2). Then, the pdf is $f(x) = 6x(1-x), 0 < x <1$.

Let g(x) be the Uniform(0,1) density, i.e., $g(x) = 1, 0 < x < 1$.

(We choose the uniform distribution on (0,1) as the proposal distribution.)

Note that $f(x) / g(x) = \frac{6x(1-x)}{1} \leq c \equiv c$

Hence, we accept if $u < f(y)/cg(y) = y(1-y)$.  

n = 1000
x = rep(0,n) # storage for generated x-values
k = 0 # count of accepted samples

while(k < n){
	u = runif(1)
    y = runif(1)
    
    if(y * (1-y) > u){
    	k = k + 1
        x[k] = y
    }
}

Remark on AR method

Let T be the number of iterations.

It can be shown that $P(Accept) = 1/c$.

Hence, the average number of the accepted samples is $T/c$.

This implies that choosing c as small as possible improves the efficiency of AR method. 

요즘에 공부하고 있는 분야가 3개 있는데,

다음과 같다.

1. 인과추론(What if + Causal Inference for Statistics, Social etc),

2. 강화학습 (Model Free control 부분)

3. Bayesian Computation with R

여기서 3번째인 Bayesian Computation 부분을 중점적으로 공부하고 있는데,

Bayesian Computation with R 교재를 사용하고 있다.

교재는 다음과 같이 생겼다.

여기서 2, 3, 4장까지는 꽤 즐겁게 공부를 했던 것 같다.

어쩌면 나는 Bayesian computation의 천재일지도 모르겠다라는 생각까지 했다. 

그러다 5장 Rejection Sampling 부분을 만나게 됐는데, 

3학년 1학기 Statistical Computation 과목에서 배운 Acceptance-Rejection method와 매우 유사하게 생겨서

'혹시 이 부분이 그 부분인가?' 하고 확인하였을 때, 수식이 똑같이 써져 있던 것을 확인할 수 있었다. 

ㅋㅋㅋㅋㅋㅋㅋㅋㅋ 그 당시에 이 부분이 제대로 이해가 안되서 난리난 흔적을 볼 수 있다. 

어찌됐든 이렇게 5장을 들어가니 Statistical Computation 과목이 얼마나 중요한 과목이었는지를 알 수 있었다. 

그리하여 이 책의 5장, Introduction to Bayesian Computation부분을 공부할 때, 

Statistical Computation 과목에서 배웠던 내용을 연계해서 이해도를 충분히 높이는 작업이 필요할 듯하다.

그래서 결심을 했는데, 

절대로 내가 학기 끝나고 안 하는 행동인 '학기 끝나고 학기동안 배웠던 과목 정리하기' 작업을 진행해볼까한다.

꽤 귀찮은 작업이 될 것 같기는하다.

아 근데 물론 한 과목만 정리할거다. Statistical Computation과목만.. ㅎㅎ

그리고 뭔가 Statistical Computation과목을 잘 정리하면 앞으로의 생활이 편해질 것 같은 느낌도 들긴 한다 음음

또 이상해보일수도 있지만 요새 이쪽으로 공부하는게 조금 재미있다.. 이해하는데 힘들때도 있는데 풀리는 쾌감이 많이 좋아서 즐길 수(?) 있는 것 같다.

오케이 그럼 앞으로 Bayesian Computation과 Statistical Computation 한번 잘 정리해봐야겠다.

파이팅~!!

+

밑에 제 사진은 내용과는 관계없이 오랜만에 티스토리 쓰는거라 넣고싶어서 넣은겁니다답

폼 한번 잡아봤습니다 네네

새 창에서 열기

새 창에서 열기