Summary

욕조물을 적정 온도 42도로 데우려고 한다. 데워졌을 것이라 판단되어 온도계로 온도를 쟀다. 단 사용한 온도계는 정밀도가 낮으며 계측된 온도 x는 실제 온도 theta 가 평균, 표준편차 2도인 정규분포를 따르는 확률분포이다. 지금 온도계가 가리키는 온도는 첫 번째가 40도, 두 번째가 41도였다.

실제 온도는 몇 도일까? ([1] p272 ~ 280)

library(rjags)
library(coda)

set.seed(2022) # seed 설정

# 모델 문자열 설정
mod_string <- " model {

  # 우도 설정: 실제 온도 theta, 표준편차 2를 따른다고 가정
  for (i in 1:n) {
    y[i] ~ dnorm(theta, 1.0/sig2) # 정밀도 = 1/분산
  }
  
  # 사전분포 설정: 
  # 적정 온도인 42도로 데우려고 하므로 평균이 42, 표준편차 3인 정규분포로 가정
  theta ~ dnorm(42, 1/3^2)
  sig2 <- 2^2
} "

# 측정된 온도값
y <- c(40, 41)
n <- length(y)

# 데이터 설정
data_jags <- list(y = y, n = n)

# 모니터링 변수 설정
# 물 온도를 추정하려는거지.
params <- c('theta')

# 초기값 함수는 옵션인데, 여기서는 측정된 값으로 설정한다.
inits <- function () {
  inits <- list('theta' = c(40))
}

# 체인은 3개를 실행
mod <- jags.model(textConnection(mod_string), data = data_jags, inits = inits, n.chains = 3)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 2
##    Unobserved stochastic nodes: 1
##    Total graph size: 12
## 
## Initializing model
update(mod, 5e3) # 번인 처리

# 사후분포 샘플 생성
mod_sim <- coda.samples(mod, variable.names = params, n.iter = 5e3)
mod_csim <- do.call(rbind, mod_sim) # 3개의 체인의 샘플을 하나로 묶음

Model Checking

plot(mod_sim) # 3개의 체인은 적절하게 수렴하는지?

gelman.diag(mod_sim) # Potential scale reduction factors는 1에 가까운지?
## Potential scale reduction factors:
## 
##       Point est. Upper C.I.
## theta          1          1
gelman.plot(mod_sim) # shrink factor 가 4000 이후로 1에 수렴함을 볼 수 있다.

autocorr.diag(mod_sim) # 초기 lag에 대한 자기상관 값을 검토
##               theta
## Lag 0   1.000000000
## Lag 1   0.002584077
## Lag 5   0.003080737
## Lag 10 -0.004432652
## Lag 50 -0.001792495
autocorr.plot(mod_sim)

effectiveSize(mod_sim) # 효과표본 크기는 샘플 수에 가까운지  검토
## theta 
## 15000

사후 분포 추정

사후분포 평균을 계산한다.

# estimation
pm_coeff <- colMeans(mod_csim) # 사후분포 평균값 계산
pm_coeff
##    theta 
## 40.78435
summary(mod_sim)
## 
## Iterations = 5001:10000
## Thinning interval = 1 
## Number of chains = 3 
## Sample size per chain = 5000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##           Mean             SD       Naive SE Time-series SE 
##       40.78435        1.28265        0.01047        0.01047 
## 
## 2. Quantiles for each variable:
## 
##  2.5%   25%   50%   75% 97.5% 
## 38.25 39.93 40.79 41.64 43.32

40.78 정도로 나온다. 책에서는 이 값이 40.77 로 계산된다.

HPD 구간을 구한다.

HPDinterval(as.mcmc(mod_csim))
##          lower    upper
## theta 38.06591 43.11818
## attr(,"Probability")
## [1] 0.95

추가적으로 실제 온도가 41 ~ 42도일 확률은?

mean(mod_csim >= 41 & mod_csim <= 42)
## [1] 0.2591333

References

[1] 고지마 히로유키, 2015, 세상에서 가장 쉬운 베이즈통계학 입문, 지상사