[통계/R] 베이지안 회귀 기초 개념 & ADP 실기 26회 R 구현

통계학을 공부하면서 가장 어려운 개념이라고 생각하는 베이지안.  

ADP 실기가 베이즈 관련 개념들을 많이 출제하는 것 같아서

나름 공부한 내용을 포스팅하기로 했다. 더 공부하면서 지속적으로 포스팅 수정할 계획.

 

베이지안은 접할수록 쉽지 않은 개념이다. 이렇게 정리를 해도 온전히 내 것으로 만들지 못했다.. 갈길이 멀다. ㅠㅠ
어렵지만 그만큼 흥미로운 분야라고 생각한다. 복습하고 많이 접할수록 언젠간 더 가까워지지 않을까.. 

---

[ 베이지안 회귀 ]

• 확률론적인 방법을 이용하여 회귀분석을 수행하는 방법
• 베이지안 회귀는 베이즈 정리를 기반으로 사전 확률분포에서 MCMC기법에 의해 선정된 관측값을 반영하여, 사후 확률분포를 추산하는 방식
• 모든 관측값이 반영되었을 때 새로운 값 관측값(unobserved data) x에 대한 출력값 y를 확률 분포에 의해 결정할 수 있는 예측 분포(Prodictive distribution)를 만드는 모델

베이지안 회귀분석의 흐름

1. 베이지안 회귀분석에서는 모델 파라미터들의 값이 어떤 분포를 따르는지를 나타내는 사전 분포(prior distribution)를 이전 경험이나 지식에 기반하여 먼저 정의

→ 만약 ADP 실기 문제가 출제 된다면 사전 분포를 미리 정의해줄 것임 (ADP 26회 기출처럼)

2.  데이터를 관측하여 이 사전 분포를 업데이트한 사후 분포(posterior distribution)를 추정

3. 추정된 사후 분포를 통해 불확실성을 고려한 예측 분포(predictive distribution)를 계산

 

일반 회귀분석 vs. 베이지안 회귀

 

장점 : 불확실성 고려

• 베이지안 회귀분석에서는 사전분포를 설정함으로써, 회귀계수에 대한 불확실성을 추정 가능
• 미래 데이터에 대한 예측분포(predictive distribution)를 추정할 수 있어 불확실성을 고려한 예측을 수행할 수 있음

🔎 불확실성을 고려한다는 의미는 무엇일까?

베이지안 회귀분석은 추정한 파라미터(모델 계수)가 실제 값과 어느 정도의 오차 범위 내에 존재할 가능성이 있는지를 확률적으로 계산하여, 이러한 추정 결과를 확률 분포 형태로 나타내기 때문에 추정한 값이 가질 수 있는 불확실성을 표현할 수 있음.

단순 회귀분석에서는 독립 변수의 값이 주어지면 종속 변수의 값이 단 하나의 값으로 결정되는 것과 달리, 베이지안 회귀분석에서는 미래 데이터에 대한 예측 분포를 추정함으로써 불확실성을 반영.

즉, 점추정 방식의 단순 회귀와는 다르게 베이지안 방법론에서는 ‘평균이 4, 어떤 분산으로 되어있는 정규분포로 파라미터가 분포되어있다’고 말하며 확률분포를 따짐.

 

 

베이즈룰

• 베이지안 회귀분석은 베이즈룰을 기반으로 하는 통계적 모델링 방법
• 베이즈룰 : 두 확률변수의 사전확률과 사후확률 사이의 관계를 나타내는 정리

 

 

베이즈룰 Example

 

마르코프 체인 몬테카를로법 (MCMC)

• 베이지안 회귀분석에서 사후 분포를 추정하는 데에 사용되는 기법
• MCMC : 마르코프 체인 + 몬테카를로 합친 알고리즘
  • 마르코프 체인 : 어떤 상태가 바로 이전 상태의 결과에 영향을 받는 상황을 의미
  • 몬테카를로 : 무작위 샘플링을 의미

 

• 초기값(initial value)을 기반으로 사후 분포를 추정하며, 이 때 생성된 샘플은 이전 샘플 값에 의존
• 이전 샘플에서 시작해 천천히 샘플 공간을 탐색하며, 목표 분포에 수렴할 때까지 샘플링을 수행
• MCMC 알고리즘을 통해 얻은 샘플은 사후 분포의 모양을 나타내며, 이를 통해 모델 파라미터의 불확실성을 추정할 수 있음

 

ADP 기출 26회 7번 : 베이지안 회귀

• Python에서 pymc3이라는 라이브러리를 통해 구현할 수 있는 것으로 아는데, 코랩과 로컬 등 여러 환경을 조정해가면서 구현하려고 했으나 라이브러리 간 충돌이 심해서 구현에 실패했다. ADP 실기 응시자들에게 시간 효율성을 위해 R을 추천한다.. 그리고 심지어 26회 제공된 패키지에는 pymc3가 없었음. python 응시자들은 애초부터 이 문제를 맞을 확률이 낮았단 말 ,,      

❓ height(키), weight(몸무게), waist(허리둘레) 컬럼을 가진 problem7.csv파일을 가지고 다음을 분석하라.

A시의 20대 남성 411명을 임의로 추출하여 키, 몸무게, 허리둘레를 조사하여 기록한 데이터이다.
이 데이터를 이용하여 20대 남성의 키와 허리둘레가 체중에 영향을 미치는지 알아보고자 한다.

7.1 아래 조건을 참고하여 회귀계수(반올림하여 소수점 두자리)를 구하시오.
- 베이지안 회귀시드넘버 1234로 지정
- 1000번의 burn-in 이후 10,000의 MCMC를 수행
- 회귀계수의 사전분포는 부적절한 균일분포(inproper uniform prior distribution), 오차항의 분산의 사전분포는 역감마 분포로 지정.
- 이때, 형상(Shape)모수와 척도(Scale)모수는 각각 0.005로 지정.

7.2 위에서 만든 모델을 바탕으로 키 180cm, 허리둘레 85cm인 남성의 몸무게를 추정하라.

• 여러 형태의 분포에 맞게 함수 옵션을 수정해주면 된다. 

data <- read.csv("https://raw.githubusercontent.com/Datamanim/datarepo/main/adp/26/problem7.csv")

library(rjags)
library(runjags)

model_string <- "
model {
  # Define variables
  for (i in 1:n) {
    y[i] ~ dnorm(mu[i], tau)
    mu[i] <- beta[1] + beta[2]*x1[i] + beta[3]*x2[i]
  }
  beta[1] ~ dunif(-1000, 1000)  # Improper uniform prior
  beta[2] ~ dunif(-1000, 1000)  # Improper uniform prior
  beta[3] ~ dunif(-1000, 1000)  # Improper uniform prior
  tau ~ dgamma(0.005, 0.005)    # Inverse gamma prior
}
"

# Set seed
set.seed(1234)

# Prepare data
data_jags = list("x1" = data$height,
                 "x2" = data$waist,
                 "y" = data$weight,
                 "n" = nrow(data))

# Run MCMC
jags_model <- jags.model(textConnection(model_string),
                         data = data_jags,
                         n.chains = 3)
update(jags_model, 1000)
samples <- coda.samples(jags_model, variable.names = c("beta", "tau"), n.iter = 10000)

summary(samples)

Output

# Output

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 411
   Unobserved stochastic nodes: 4
   Total graph size: 1736

Initializing model

  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100%
> update(jags_model, 1000)
  |**************************************************| 100%
> samples <- coda.samples(jags_model, variable.names = c("beta", "tau"), n.iter = 10000)
  |**************************************************| 100%
> summary(samples)

Iterations = 2001:12000
Thinning interval = 1 
Number of chains = 3 
Sample size per chain = 10000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

            Mean      SD  Naive SE Time-series SE
beta[1] -26.5886 8.75340 0.0505378        3.03867
beta[2]   0.4824 0.27197 0.0015702        0.07157
beta[3]   0.2143 0.54911 0.0031703        0.12431
tau       0.2491 0.05628 0.0003249        0.01122

2. Quantiles for each variable:

            2.5%      25%      50%      75%   97.5%
beta[1] -42.7128 -32.2536 -26.5532 -19.8254 -9.5844
beta[2]  -0.1515   0.3367   0.5501   0.6677  0.8695
beta[3]  -0.5297  -0.2884   0.1777   0.4790  1.3939
tau       0.1316   0.2135   0.2567   0.2912  0.3410

# 7-2 문제 해답
키가 180cm이고 허리둘레가 85cm인 남성의 몸무게를 추정하려면 다음과 같이 계산

mu = beta[1] + beta[2]180 + beta[3]85
mu = -26.59 + 0.48180 + 0.2185
mu = 74.14

따라서 추정된 몸무게는 약 74.14kg

---

개인적으로 공부하며 요약한 자료입니다.

오류 발견 시 댓글 남겨주시면 정정하겠습니다!    

 

출처는 더보기에 ↓

Reference

• 파이썬 코드 관련 설명
  • https://www.pymc.io/projects/docs/en/stable/learn.html
  • Bayesian Additive Regression Trees: Introduction — PyMC example gallery
  • https://datascienceschool.net/03 machine learning/19.02 베이지안 회귀 분석 예제.html#id3
  • Pyro(PyTorch 기반의 확률 프로그래밍 라이브러리)를 이용한 베이지안 회귀분석
• R 코드 및 통계 설명
  • https://www.youtube.com/watch?v=z7rCogOaDh4
  • https://www.youtube.com/watch?v=_V6z8xfuk70
  • https://www.r-bloggers.com/2010/08/using-jags-in-r-with-the-rjags-package/
• 통계적인 설명
  • https://mons1220.tistory.com/212
  • https://coffeewhale.com/bayesian/linear/regression/2019/10/19/bayesian-lr/
  • https://wikidocs.net/122767
  • https://diy-multitab.tistory.com/36
  • https://ko.wikipedia.org/wiki/베이즈_정리