Exercise 2.19 (Cuckoo birds redux)

# Load packages
library(bayesrules)
library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(janitor)

## 
## Attaching package: 'janitor'
## 
## The following objects are masked from 'package:stats':
## 
##     chisq.test, fisher.test

library(ggplot2)
library(patchwork)

# Define posibles π(tasa de éxito de sobrevivir las crías al menos una semana)
Birds <- data.frame(pi = c(0.6, 0.65, 0.7, 0.75))

# Define modelo de probabilidades a priori f(π)
prior <- c(0.3, 0.4, 0.2, 0.1)

# Simular 10000 valores de π con modelo de probabilidades a priori
set.seed(131184)
Birds_sim <- sample_n(Birds, size = 10000, weight = prior, replace = TRUE)

# Simular 10000 valores y (número de éxitos en n ensayos)  utilizando la función rbinom() con size = 15 y prob = π y agregar al data.frame
Birds_sim <- Birds_sim %>% 
  mutate(y = rbinom(10000, size = 15, prob = pi))

# Chequear Birds_sim
Birds_sim %>% 
  head(5)

##     pi  y
## 1 0.70 11
## 2 0.70 13
## 3 0.70 10
## 4 0.65 11
## 5 0.70 10

Aproximaciones del modelo de probabilidad priori f(π):

# Aproximación de probabilidades a priori f(π)

Birds_sim %>% 
  tabyl(pi) %>% 
  adorn_totals("row")

##     pi     n percent
##    0.6  2994  0.2994
##   0.65  4055  0.4055
##    0.7  1959  0.1959
##   0.75   992  0.0992
##  Total 10000  1.0000

Las aproximaciones mediante simulación para el modelo de probabilidad priori f(π) son las siguientes: f(π=0.6)=0.2994; f(π=0.65)=0.4055; f(π=0.7)=0.1959; f(π=0.75)=0.0992.

Histograma correspondiente a la simulación de f(π):

#Histograma aproximaciones de probabilidades a priori f(π)
ggplot1 <-ggplot(Birds_sim, aes(x = pi)) +ggtitle("Simulación del modelo a priori f(π) ") +
  stat_count(aes(y = ..prop..))
ggplot1

## Warning: The dot-dot notation (`..prop..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(prop)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

# Plot y por cada π (proporciones)
ggplot(Birds_sim, aes(x = y)) + 
  stat_count(aes(y = ..prop..)) + 
  facet_wrap(~ pi)

# Nos enfocamos en simulaciones en donde obtuvimos y = 10
y_ten <- Birds_sim %>% 
  filter(y == 10)

Aproximaciones del modelo de probabilidad posteriori f(π|y=10):

# Aproximación de probabilidades a posteriori f(π|y=10)
y_ten %>% 
  tabyl(pi) %>% 
  adorn_totals("row")

##     pi    n    percent
##    0.6  559 0.28755144
##   0.65  834 0.42901235
##    0.7  385 0.19804527
##   0.75  166 0.08539095
##  Total 1944 1.00000000

Las aproximaciones mediante simulación para el modelo de probabilidad posteriori f(π|y=10) son las siguientes: f(π=0.6|y=10)≈0.2876; f(π=0.65|y=10)≈0.429; f(π=0.7|y=10)≈0.198; f(π=0.75|y=10)≈0.0854.

Histograma correspondiente a la simulación de f(π|y=10):

# Plot aproximaciones de probabilidades a posteriori f(π|y=10)
# Histograma aproximaciones de probabilidades a posteriori f(π|y=10)
ggplot2 <-ggplot(y_ten, aes(x = pi)) + ggtitle("Simulación del modelo a posteriori f(π|y=10)") +
  stat_count(aes(y = ..prop..))
ggplot2

Histogramas juntos correspondientes a la simulaciones de f(π) y f(π|y=10):