Simulation

Théo Boulakia

4 décembre 2024

Packages

library(tidyverse)

Seed

set.seed(84735)

Nombre de simulations

n = 10000

Prior

ggplot() +
  stat_function(fun = ~dbeta(.x, 9, 11))

Simulation

x = data.frame(pi = rbeta(n, 9, 11)) |> 
  mutate(y = rbinom(n, size = 50, prob = pi))

Plot

ggplot(x) +
  geom_point(aes(x = pi, y = y, color = (y == 30)))

Posterior data

posterior = x |> 
  filter(y == 30)

Plot posterior

Code 
ggplot(posterior) +
  geom_density(aes(x = pi))

Compare

Code 
ggplot(posterior) +
  geom_density(aes(x = pi)) + 
  stat_function(fun = ~dbeta(.x, 39, 31), color = "red")

Summarise

posterior |> 
  summarise(mean(pi), sqrt(var(pi)))
   mean(pi) sqrt(var(pi))
1 0.5549714    0.05779086

Compare

mean = 39 / (39 + 31)
sd = sqrt((39 * 31) / ((39 + 31) ^ 2 * (39 + 31 + 1)))
mean
[1] 0.5571429
sd
[1] 0.05895029