Solution to exercise 3M6 of Statistical Rethinking book

3M6.R

library(rethinking)

set.seed(42)
p_grid <- seq(0, 1, length.out=1000)
prior <- rep(1, 1000)
real_water <- 0.7

posteriors <- function(N) {
  likelihod <- dbinom(round(N*real_water), size=N, prob=p_grid)
  posterior <- likelihod * prior
  posterior <- posterior/sum(posterior)
  return(posterior)
}

for (i in 1:3000) {
  posterior <- posteriors(i)
  samples <- sample(p_grid, prob=posterior, size=1e4, replace=T)
  interval <- PI(samples, prob=0.99)
  names(interval) <- NULL
  diff <-interval[2] - interval[1]
  print(diff)
  if (diff < 0.05) {
    print(i)
    break
  }
}