다음 영상을 보고 R코드를 만들어 달라고 chatGPT에게 요청
chatGPT에게 한 요청
흉내 낼 목표 분포 f(x)는 N(-2,1)+N(2,1) 이고,
Xt = 1.4,
Xt+1 = N(Xt, 0.4),
비율값을 저장하는 알파 = f(Xt+1) / f(Xt),
u = 균등분포(0,1),
if 알파 > u 이면 Xt를 저장 하고 아니면 버림.
위 동작을 MCMC라고 하는 것 같은데
R로 구현해 줘
단, R에서 사용하는 함수는 사용하지 말 것.
그리고, 저장된 Xt의 히스토그램을 그려줘
그랬더니 다음 코드를 줌
-----------------------------------------
# MCMC 구현
rm(list = ls())
# 두 개의 정규 분포를 합친 목표 분포 f(x)
f <- function(x) {
exp(-0.5 * (x + 2)^2) + exp(-0.5 * (x - 2)^2)
}
# Metropolis-Hastings MCMC 알고리즘 구현
set.seed(123) # 재현성을 위한 시드 설정
iterations <- 1000000 # 반복 횟수
Xt <- 1.4 # 초기값
samples <- numeric(iterations) # 샘플 저장할 공간
samples[1] <- Xt
for (t in 2:iterations) {
# Xt+1 샘플링 (N(Xt, 0.4)로부터)
Xt_candidate <- Xt + rnorm(1, 0, sqrt(0.4))
# 비율값 계산
alpha <- f(Xt_candidate) / f(Xt)
# 수락 여부 결정
u <- runif(1)
if (alpha > u) {
Xt <- Xt_candidate # 샘플을 수락
}
# 샘플 저장
samples[t] <- Xt
}
# 저장된 샘플의 히스토그램 그리기
hist(samples, breaks = 50, probability = TRUE, col = "lightblue", main = "MCMC 샘플 히스토그램")
댓글 0