본문으로 바로가기



능선회귀분석 with R - 4. 베이지안 MAP 추정량으로서의 능선회귀분석

category Statistics/Linear models 2017.06.18 01:20
능선회귀분석 with R - 4. 베이지안 MAP 추정량으로서의 능선회귀분석

(이 포스팅은 The Elements of Satistical Learning의 흐름을 따라가고 있음을 밝힙니다. 책은 링크가 없어진 관계로 좀 더 쉬운 버전 링크를 공유합니다. pdf파일은 공개되어 있으니, 다운받아 참고하시면 좋을것 같습니다.)

이전 포스팅에서 우리는 실제 데이터를 가지고, 람다의 변화에 따른 능선회귀의 계수를 나타내주는 coefficient path를 그려보았다. 오늘은 조금 이론적인 내용으로 능선회귀의 계수, \[ \hat{\beta}_{Ridge}=\left(X^{T}X+\lambda I\right)^{-1}X^{T}y \] 를 베이지안적인 관점에서 어떻게 바라볼 수 있는지에 대하여 다뤄보고자 한다.

베이지안 MAP(Maximum a posteriori estimation)으로서의 능선회귀

이 포스팅은 The Elements of Satistical Learning의 Ex. 3.6.을 풀면서 작성하게 되었다. 먼저 문제는 다음과 같다.


Show that the ridge regression estimate is the mean (and mode) of the posterior distribution, under a Gaussian prior \(\beta\sim N\left(0,\tau^{2}\mathbf{I}\right)\) and Gaussian sampling model \(y\sim N\left(\mathbf{X}\beta,\sigma^{2}\mathbf{I}\right)\). Find the relationship between the regularization parameter \(\lambda\) in the ridge formula, and the variance \(\tau^2\) and \(\sigma^2\).


사실 위의 문제는 조금 불친절하게 기술되어 있다고 생각할 수 있는데, 그 이유는 문제의 서술에 있어서 어느정도 독자가 베이지안 통계에 대한 지식을 가지고 있다는 전제하에 쓴 것 이기 때문이라 생각한다.

베이지안 통계에서는 모델의 모수에 대하여 알고있는 ‘사전지식’을 ’prior’, 즉, 우리말로 바꾸면 ’이전의’라는 뜻을 가지고 있는 단어를 사용하여 나타낸다. 따라서 ’prior distribution’이라는 말은 사용자(혹은 통계학자)가 생각하는 모수의 분포에 대한 정보를 나타낸다. 따라서 위의 문제에서

  • under a Gaussian prior \(\beta\sim N\left(0,\tau^{2}\mathbf{I}\right)\)

이라는 구절이 의미하는 내용은, 각 독립변수의 모수들, \(\beta_1, ..., \beta_p\),이 \(\sim N\left(0,\tau^{2}\mathbf{I}\right)\)을 따른다는 의미이며, 이것을 바탕으로 다음의 두 가지 의미를 생각할 수 있다.

  1. 모든 베타들은 0 근처에 위치하고 있는 가능성이 가장 높다. 이는 회귀분석을 함에 있어서 기본적인 전제가 모든 독립변수들이 종속변수와 관계가 없다는 것에서 출발하겠다는 의미가 된다.
  2. 또한 베타 prior 분포의 공분산 행렬이 대각행렬이라는 의미는 각 베타끼리의 종속성이 없다는 것을 의미한다. 즉, 어느 독립변수의 계수 \(\beta_i\)가 특정값을 가질 경우, 다른 계수들이 이것에 영향을 받지 않는다고 사용자(혹은 통계학자)는 생각한다.

이러한 사전 지식을 바탕으로, 종속변수 \(y\)에 대한 분포를 가정하는데, 이를 다음과 같이 표현하고 있다.

  • Gaussian sampling model \(y\sim N\left(\mathbf{X}\beta,\sigma^{2}\mathbf{I}\right)\)

즉, 특정 베타값이 주어졌을때, 종속변수 \(y\)는 평균이 \(\mathbf{X}\beta\)이고, 분산이 \(\sigma^{2}\mathbf{I}\)을 따른다고 가정한다. 이를 바탕으로 우리는 다음의 두 가지를 유추할 수 있다.

  1. 데이터 행렬 \(\mathbf{X}\)np열이고, 베타계수가 p1열이므로, y는 n 변량 정규분포(n variate normal distribution)를 따르게 된다. 즉, 우리가 관찰한 n개의 관측치를 특정 분포를 따르는 확률변수에서 샘플을 뽑았다고 생각한다.
  2. 역시 공분산 행렬의 대각행렬을 제외한 나머지 값들이 0이라는 것은, n개의 관측치들이 독립적이라는 것을 의미한다.

이러한 측면에서 위의 가우시안 계층모형은 다음과 같이 표현하는것이 좀 더 의미가 명확해진다. |는 ’given’을 의미한다. (또한, 이 표현이 베이지안 통계에서 좀 더 보편화된 기호라고 생각한다.)

\[ \boldsymbol{\beta}\sim N\left(0,\tau^{2}\mathbf{I}\right) \\ y|\beta\sim N\left(\mathbf{X}\beta,\sigma^{2}\mathbf{I}\right) \]

베이지안 통계에서는 이러한 계층모델이 주어졌을때(결국 그냥 세팅일 뿐이다.) posterior 분포를 구하는 것에 관심이 있다. posterior는 사후분포라고 번역이 되는데, 이것은 주어진 prior와 sampling model을 바탕으로 데이터가 주어졌을때 기존의 지식(prior)과 현재 주어진 정보(데이터)를 동시에 고려한 이후의 모수에 분포를 의미하기 때문이다.

이 분포를 구하는 것은 수식적으로는 상당히 계산이 복잡하지만 아이디어는 결국 기초 통계 시간에 배운 유명한 베이즈 정리를 바탕으로 한다. Sampling model을 \(p\left(y|\beta\right)\)이라 하고, prior분포를 \(p\left(\beta\right)\)라 가정할 때, 베이즈 정리에 의하여 posterior분포는 다음과 같이 구할 수 있다.

\[ p\left(\beta|y\right)=\frac{p\left(\beta,y\right)}{p\left(y\right)}=\frac{p\left(\beta\right)p\left(y|\beta\right)}{p\left(y\right)} \]

여기서 중요한 점은, 사후 분포를 구할 때, 위의 식을 그대로 이용하는 것이 아니라는 점이다. `사후분포를 구할때는 보통 마지막 식의 분자에 위치한 \(p\left(\theta\right)p\left(y|\theta\right)\) 두 항만을 이용한다.

우리가 관심있는 사후 분포는 베타의 분포가 우리가 이전에 알고있던 모든 베타가 0 근처에 몰려있고, 각 베타의 확정성을 의미하는 분산 \(\tau\)의 값이 어떻게 변하게 되는지만 관심이 있기 때문이다. 마지막 항의 분모에 위치한 \(p\left(y\right)\)의 주된 임무(?)는 분자를 적정한 값으로 나눠주어 \(\frac{p\left(\theta\right)p\left(y|\theta\right)}{p\left(y\right)}\) 전체의 값이 확률분포의 조건, 적분해서 1이 나오는 성질을 만족시키도록 하는 것에 있다. 따라서 베이지안에서는 보통 관심있는 항만을 표시하기 위하여 다음과 같이 비례한다라는 의미를 가지고 있는 기호 \(\propto\)를 사용한다. (영어로는 ’proportional to’의 의미이다.)

\[ p\left(\beta|y\right)\propto p\left(\beta\right)p\left(y|\beta\right) \]

위 두 항을 가지고도 베타의 사후분포가 어떻게 되는지에 대한 충분히 파악할 수 있기 때문이기도 하다. 자, 이제 이것을 바탕으로 우리가 오늘 구하고 싶은 문제를 풀어보도록 하겠다. 먼저, 베타의 사전분포 pdf와 베타가 주어졌을때의 y의 conditional pdf는 다음과 같다.

\[ \begin{align*} p\left(\beta\right) & =\frac{1}{\sqrt{\left(2\pi\tau\right)^{p}}}exp\left(-\frac{1}{2}\beta^{T}\left(\tau^{2}\mathbf{I}\right)^{-1}\beta\right)\\ p\left(y|\beta\right) & =\frac{1}{\sqrt{\left(2\pi\sigma\right)^{n}}}exp\left(-\frac{1}{2}\left(y-\mathbf{X}\beta\right)^{T}\left(\sigma^{2}\mathbf{I}\right)^{-1}\left(y-\mathbf{X}\beta\right)\right) \end{align*} \]

위의 식에 대한 정보는 위키의 multivariate normal distribution 페이지를 참고하도록 하자. 이 두개의 식을 가지고 우리가 관심있어 하는건 이 두개의 정규분포 pdf의 곱이 베타를 기준으로 어떠한 분포를 따르는가에 대한 것이다. 다행인 것은 우리는 두개의 정규분포 pdf를 곱하였을 경우 이 분포는 다시 정규분포가 된다는 것을 알고 있다. 요즘 유행하는 팩트 체크를 해보면, 두 정규분포 pdf의 곱은 다시 어떤 정규분포의 pdf폼을 따른다. 하지만 두 정규분포를 따르는 확률변수의 곱은 정규분포를 따르지 않는다.

따라서, 사후 분포의 평균을 나타내는 모수와 분산을 나타내는 모수만을 알아내면 그만인 것이다. 또한, 베타의 사후분포는 여전히 p변량 정규분포를 따를 것이므로, 베타 사후분포의 pdf는 일반적인 p변량 정규분포의 pdf식인 다음의 꼴을 가지고 있을 것이다.

\[ p\left(\beta|y\right) =\frac{1}{\sqrt{\left(2\pi\right)^{p}\left|variance\right|}}exp\left(-\frac{1}{2}\left(\beta-mean\right)^{T}\left(variance\right)^{-1}\left(\beta-mean\right)\right) \]

다행스럽게도(?) 정규분포의 pdf는 지수함수 안에 평균과 분산을 나타내는 모수들이 들어있기 때문에, 지수함수 밖에 있는 항들을 신경쓰지 않아도 어떤 부분이 평균벡터인지, 어떤 부분이 공분산 행렬인지를 판단할 수 있다. 즉, 아래에서 처럼 우리가 위에 주어진 두 pdf를 곱한 것의 지수함수 안의 부분을 베타로 묶었을 경우, \(\beta^{T}\left(\right)\beta\), 괄호안에 위치한 부분이 공분산 행렬의 역행렬이 되는 것이고, \(-2\beta^{T}\left(\right)\)의 괄호부분은 공분산 행렬의 역행렬과 평균벡터의 곱한 값이 위치하게 될 것이다.

\[ \begin{align*} \left(\beta-\mu\right)^{T}\Sigma^{-1}\left(\beta-\mu\right) & =\beta^{T}\Sigma^{-1}\beta-\beta^{T}\Sigma^{-1}\mu-\mu^{T}\Sigma^{-1}\beta+\mu^{T}\Sigma^{-1}\mu\\ & =\beta^{T}\Sigma^{-1}\beta-2\beta^{T}\Sigma^{-1}\mu+\mu^{T}\Sigma^{-1}\mu \end{align*} \]

위의 사실을 이용하여 사후분포의 평균벡터와 공분산행렬을 알아내어 보자.

\[ \begin{align*} p\left(\beta|y\right) & \propto p\left(\beta\right)p\left(y|\beta\right)\\ & \propto exp\left(-\frac{1}{2}\beta^{T}\left(\mathbf{I}\right)^{-1}\beta\right)exp\left(-\frac{1}{2}\left(y-\mathbf{X}\beta\right)^{T}\left(\sigma^{2}\mathbf{I}\right)^{-1}\left(y-\mathbf{X}\beta\right)\right)\\ & \propto exp\left(-\frac{1}{2\tau^{2}}\beta^{T}\beta-\frac{1}{2\sigma^{2}}\left(y-\mathbf{X}\beta\right)^{T}\left(y-\mathbf{X}\beta\right)\right)\\ & \propto exp\left(-\frac{1}{2\tau^{2}}\beta^{T}\beta-\frac{1}{2\sigma^{2}}\left(y^{T}y-2\beta^{T}\mathbf{X}^{T}y+\beta^{T}\mathbf{X}^{T}\mathbf{X}\beta\right)\right)\\ & \propto exp\left(-\frac{1}{2\tau^{2}}\beta^{T}\beta+\frac{1}{\sigma^{2}}\beta^{T}\mathbf{X}^{T}y-\frac{1}{2\sigma^{2}}\beta^{T}\mathbf{X}^{T}\mathbf{X}\beta\right)\\ & \propto exp\left(\beta^{T}\left(-\frac{1}{2\tau^{2}}\mathbf{I}-\frac{1}{2\sigma^{2}}\mathbf{X}^{T}\mathbf{X}\right)\beta-2\beta^{T}\left(-\frac{1}{2\sigma^{2}}\mathbf{X}^{T}y\right)\right)\\ & \propto exp\left(\beta^{T}\left(-\frac{1}{2\sigma^{2}}\left(\mathbf{X}^{T}\mathbf{X}+\frac{\sigma^{2}}{\tau^{2}}\mathbf{I}\right)\right)\beta-2\beta^{T}\left(-\frac{1}{2\sigma^{2}}\mathbf{X}^{T}y\right)\right) \end{align*} \]

위의 전개에서 \(y^T y\)항은 우리가 관심있는 베타와는 관계가 없으므로 생략이 되었다. 마지막의 결과로 우리는 공분산 행렬의 역행렬과 그것과 평균 벡터와의 곱이 무엇인지 알게 되었다. 따라서 우리가 찾는 베타 사후 분포의 평균 벡터(\(\mu\))는 다음과 같다.

\[ \begin{align*} \Sigma^{-1}\mu & =-\frac{1}{2\sigma^{2}}\mathbf{X}^{T}y\\ \Rightarrow\mu & =\Sigma\left(-\frac{1}{2\sigma^{2}}\mathbf{X}^{T}y\right)\\ & =\left(-2\sigma^{2}\left(\mathbf{X}^{T}\mathbf{X}+\frac{\sigma^{2}}{\tau^{2}}\mathbf{I}\right)^{-1}\right)\left(-\frac{1}{2\sigma^{2}}\mathbf{X}^{T}y\right)\\ & =\left(\mathbf{X}^{T}\mathbf{X}+\frac{\sigma^{2}}{\tau^{2}}\mathbf{I}\right)^{-1}\mathbf{X}^{T}y \end{align*} \]

어디서 많이 봤던 꼴이다. 그렇다. 주어진 베이지안 계층 모형에서의 베타 모수의 사후분포의 평균은 우리가 유도했던 능선회귀의 해, \[ \hat{\beta}_{Ridge}=\left(X^{T}X+\lambda I\right)^{-1}X^{T}y \] 에서 \(\lambda\)부분을 \(\frac{\sigma^{2}}{\tau^{2}}\)로 바꾼것과 동일하다. 또한 정규분포는 대칭인 분포이기 때문에 pdf의 최고점을 의미하는 mode와 평균이 일치한다. 따라서 사후분포가 정규분포를 따를 경우, 사후분포의 pdf의 mode 사용하는 베이지안의 MAP(Maximum a posteriori) 추정량은 정규분포의 평균과 동일하게 된다. 이러한 관점에서 능선회귀의 베타 추정은 (~여러 베이지안들의 관점에서는 지극히~)베이지안 취향의 추정방식이라고 생각할 수 있는 것이다.

Reference

[1] Friedman, Jerome, Trevor Hastie, and Robert Tibshirani. The elements of statistical learning. Vol. 1. Springer, Berlin: Springer series in statistics, 2001.

[2] MIT OpenCourseWare - Probabilistic Modeling and Bayesian Analysis (https://ocw.mit.edu/courses/sloan-school-of-management/15-097-prediction-machine-learning-and-statistics-spring-2012/lecture-notes/MIT15_097S12_lec15.pdf)

[3] Lyndon White (https://math.stackexchange.com/users/1505/lyndon-white), Product of two Gaussian PDFs is a Gaussain PDF, but Produt of two Gaussan Variables is not Gaussian, URL (version: 2017-04-13): https://math.stackexchange.com/q/1112866

SHARE TO