문제풀이

[익명(223.131)]

# -----------------------------

# 입력값

# -----------------------------

n1 <- 20

n2 <- 23

ybar1 <- 10459

ybar2 <- 10974.6919

SSp <- 461598

# -----------------------------

# 사후분포 모수

# -----------------------------

md <- ybar1 - ybar2           # m_d'

nu <- n1 + n2 - 2            # 자유도

# -----------------------------

# Monte Carlo 반복횟수

# -----------------------------

N <- 1000000

# -----------------------------

# 역카이제곱에서 sigma^2 생성

# sigma^2 = nu*SSp / ChiSq(nu)

# -----------------------------

sigma2 <- SSp / rchisq(N, df = nu)

hist(sigma2,breaks=200,  probability = TRUE

     , main = expression("역카이제곱으로 구한 "~sigma^2~"분포"))

lines(

  density(sigma2),

  lwd = 2, col="red"

)

# -----------------------------

# 조건부 정규분포에서 mu_d 생성

# -----------------------------

mu_d <- rnorm(

  N,

  mean = md,

  sd = sqrt(sigma2 * (1/n1 + 1/n2))

)

hist(mu_d, breaks=200,  probability = TRUE

     , main = expression(sigma^2~"를 정규분포에 넣어서 구한 평균차"~mu[d]~"분포"))

lines(

  density(mu_d),

  lwd = 2, col="red"

)

# -----------------------------

# 결과 확인

# -----------------------------

cat("mu_d 평균 =", mean(mu_d), "\n")

cat("mu_d 표준편차 =", sd(mu_d), "\n")

# 비교

# t분포로 직접 생성

mu_d_t <- md +

  sqrt((SSp/nu) * (1/n1 + 1/n2)) *  rt(N, df = nu)

# 밀도 비교

plot(

  density(mu_d),

  lwd = 2,

  main = expression(mu[d]~": 역카이제곱 + 정규 = 계층적 vs 스튜던트티 차이")

)

lines(

  density(mu_d_t),

  lwd = 2,

  lty = 2, col="red"

)