Statistical Rethinking 제3장 연습문제 풀이

3E1~3E3 풀이

문제에서 p는 "proportion of water"이다.

p_grid = seq(from=0, to=1, length.out=1000)
prior = rep(1, 1000)
likelihood = dbinom(x=6, size=9, prob=p_grid)
posterior = likelihood * prior
posterior = posterior / sum(posterior)

set.seed(100)
samples = sample(p_grid, prob=posterior, replace=TRUE, size=10000)

# 3E1
mean(samples < 0.2)
# 4e-04

# 3E2
mean(samples > 0.8)
#0.1116

# 3E3
mean(samples>0.2 & samples<0.8)
# 0.888

# 3E4
quantile(samples, probs=0.2)
#       20% 
# 0.5185185 

# 3E5
quantile(samples, probs=0.8)
#        80% 
#  0.7557558 

# 3E6
library(rethinking)
HPDI(samples, prob=0.66)
#      |0.66     0.66| 
#  0.5085085 0.7737738 

# 3E7
PI(samples, prob=0.66)
#        17%       83% 
#  0.5025025 0.7697698 

library(rethinking)은 책 서론 installing the rethkinking R package를 참조하여 rethinking package를 설치한 후에 실행하여야 한다.

 

3M1~3E5 풀이

rm(list = ls())

# 3M1
p_grid = seq(from=0, to=1, length.out=1000)
prior = rep(1, 1000)
likelihood = dbinom(x=8, size=15, prob=p_grid)
posterior = likelihood * prior
posterior = posterior / sum(posterior)

# 3M2
set.seed(100)
samples = sample(p_grid, prob=posterior, replace=TRUE, size=10000)

library(rethinking)
HPDI(samples, prob=0.90)

# 3M3
w = rbinom(10000, size=15, prob=samples)
mean(w == 8)


# 3M4
w = rbinom(10000, size=9, prob=samples)
mean(w == 6)


# 3M5 - 1
p_grid = seq(from=0, to=1, length.out=1000)
prior = ifelse(p_grid < 0.5, 0, 1)
likelihood = dbinom(x=8, size=15, prob=p_grid)
posterior = likelihood * prior
posterior = posterior / sum(posterior)

# 3M5 - 2
set.seed(100)
samples = sample(p_grid, prob=posterior, replace=TRUE, size=10000)
HPDI(samples, prob=0.9)

# 3M5 - 3
w = rbinom(10000, size=15, prob=samples)
mean(w == 8)

# 3M5 - 4
w = rbinom(10000, size=9, prob=samples)
mean(w == 6)

 

3M6은 문제가 좀 이상함...

지구에 물의 비율을 정확히 추정하는것이 아니라 지구의 실제 물의 비율 0.7을 정확히 추정한다고 해야 문제가 성립.

지구본을 던지는 횟수를 N이라 놓고 N에 따라 PI가 얼마나 달라지는지 눈으로 확인해가면서 N의 근사값을 구할수 있다.

# 3M6
f = function(N) {
  p_true= 0.7
  w = rbinom(1, size=N, prob=p_true)
  p_grid = seq(from=0, to=1, length.out=1000)
  prior = rep(1, 1000)
  likelihood = dbinom(w, size=N, prob=p_grid)
  posterior = likelihood * prior
  posterior = posterior / sum(posterior)
  samples = sample(p_grid, prob=posterior, size=10000, replace=TRUE)
  PI99 = PI(samples, 0.99)
  as.numeric(PI99[2] - PI99[1])
}

위와 같이 스크립트를 작성하고 m6.R이라고 저장한다.

> rm(list = ls())
> setwd("d:/tmp/rcode/")
> getwd()
[1] "d:/tmp/rcode"
> source("m6.R")
> f(10)
[1] 0.5705806
> f(100)
[1] 0.2272272
> f(1000)
[1] 0.07507508
> f(10000)
[1] 0.02302803
> f(5000)
[1] 0.03403403
> f(2500)
[1] 0.04904905

대략 2500번 정도 실행하면 된다.