Bayesian Statistics 3판
############################################
# Example 17.1
# Student-t posterior simulation
############################################
# 데이터
y <- c(
45.6, 41.1, 44.5, 44.0, 40.6,
44.1, 39.0, 39.5, 39.5, 41.7,
42.5, 42.7, 42.1, 42.4, 44.8,
41.0, 39.9, 43.9, 41.3, 45.1,
42.0, 38.5, 42.6, 43.8, 43.0
)
############################################
# 기본 통계량
############################################
n <- length(y)
ybar <- mean(y)
SSy <- sum((y - ybar)^2)
############################################
# 1. Amber (Jeffreys)
############################################
S_amber <- SSy
kappa_amber <- n - 1
n_amber <- n
m_amber <- ybar
sigma2_amber <- (S_amber / kappa_amber)/n_amber
sqrt_amber = sqrt(sigma2_amber)
############################################
# 2. Brett (Exact)
############################################
S <- 4.094 #사전분산 9 곱하기 카이제곱표에서 자유도 1의 50%값 0.4549
kappa <- 1
m <- 40
n0 <- 1 #S의 자유도값
S_brett <- S + SSy
kappa_brett <- kappa + n
n_brett <- n0 + n
m_brett <- (n0 * m + n * ybar) / (n0 + n)
sigma2_brett <- (S_brett / kappa_brett)/n_brett
sqrt_brett = sqrt(sigma2_brett)
############################################
# 3. Chandra (Approx)
############################################
S_chandra <- S + SSy
kappa_chandra <- kappa + n - 1
n_chandra <- n0 + n
m_chandra <- (n0 * m + n * ybar) / (n0 + n)
sigma2_chandra <- (S_chandra / kappa_chandra)/n_chandra
sqrt_chandra = sqrt(sigma2_chandra)
############################################
# 결과 출력
############################################
cat("m_amber =", m_amber, "\n")
cat("sqrt_amber =", sqrt_amber, "\n\n")
cat("m_brett =", m_brett, "\n")
cat("sqrt_brett =", sqrt_brett, "\n\n")
cat("m_chandra =", m_chandra, "\n")
cat("sqrt_chandra =", sqrt_chandra, "\n\n")
############################################
# 스튜던트 t 난수 발생
############################################
N <- 1e6
cat("Generating Student-t random numbers...\n")
# Amber
mu_amber <- rt(N, df = kappa_amber) *
sqrt_amber +
m_amber
# Brett
mu_brett <- rt(N, df = kappa_brett) *
sqrt_brett +
m_brett
# Chandra
mu_chandra <- rt(N, df = kappa_chandra) *
sqrt_chandra +
m_chandra
############################################
# 밀도 계산
############################################
j<- density(mu_amber)
E<- density(mu_brett)
A<- density(mu_chandra)
############################################
# 밀도함수 plot
############################################
plot( j, lwd = 2, main="amber, brett, chandra 밀도곡선")
############################################
# 개별 밀도도 추가
############################################
lines(E, lwd = 2, lty = 2, col="blue")
lines(A, lwd = 2, lty = 3, col="red")
legend("topright",
legend = c("Amber", "Brett", "Chandra"), # 범례에 들어갈 이름 (데이터 순서에 맞게 변경)
col = c("black", "blue", "red"), # 각 선의 색상 (j, E, A 순서)
lty = c(1, 2, 3), # 각 선의 종류 (j=실선, E=lty 2, A=lty 3)
lwd = 2) # 선 두께를 범례에도 반영
############################################
# 95% credible intervals
############################################
cat("\n===== 95% Credible Intervals =====\n")
cat("\nJeffreys:\n")
print(quantile(mu_amber, c(0.025, 0.975)))
cat("\nExact:\n")
print(quantile(mu_brett, c(0.025, 0.975)))
cat("\nApprox:\n")
print(quantile(mu_chandra, c(0.025, 0.975)))
댓글 0