sim5.1

library(purrrfect)

Attaching package: 'purrrfect'

The following objects are masked from 'package:base':

    replicate, tabulate

library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.2.1     ✔ readr     2.2.0
✔ forcats   1.0.1     ✔ stringr   1.6.0
✔ ggplot2   4.0.3     ✔ tibble    3.3.1
✔ lubridate 1.9.5     ✔ tidyr     1.3.2
✔ purrr     1.2.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

N <- 10000

simulate_coin_draw <- \() {
  tosses <- sample(c(1, 0), size = 3, replace = TRUE) # 1 = Heads, 0 = Tails
  first_h <- match(1, tosses)
  
  y <- case_when(
    is.na(first_h) ~ -1,
    first_h == 1  ~ 1,
    first_h == 2  ~ 2,
    first_h == 3  ~ 3
  )
  return(c(sum(tosses), y))
}

sim_coins <- replicate_int(N, simulate_coin_draw()[1], .as = X) %>% 
  mutate(
    Y = map_int(X, \(x) {
      simulate_coin_draw()[2]
    })
  )

sim_coins_summary <- sim_coins %>% 
  summarize(cnt = n(), .by = c(X, Y)) %>% 
  mutate(
    joint_prob = cnt / N,
    prob_X_given_Y = cnt / sum(cnt), .by = Y
  )

sim_coins_table <- sim_coins_summary %>% 
  pivot_wider(
    names_from = Y,
    values_from = joint_prob,
    names_prefix = 'y=',
    id_cols = X,
    values_fill = 0
  )

print(sim_coins_table)

# A tibble: 4 × 5
      X  `y=2`  `y=3`  `y=1` `y=-1`
  <int>  <dbl>  <dbl>  <dbl>  <dbl>
1     1 0.0936 0.0479 0.191  0.0448
2     3 0.0325 0.0154 0.0624 0.018 
3     2 0.0912 0.0479 0.183  0.0493
4     0 0.0319 0.0146 0.0591 0.0175

p_X_given_Y1 <- sim_coins_summary %>% 
  filter(Y == 1) %>% 
  select(X, Y, cnt, prob_X_given_Y)

print(p_X_given_Y1)

# A tibble: 4 × 4
      X     Y   cnt prob_X_given_Y
  <int> <int> <int>          <dbl>
1     2     1  1833          0.370
2     1     1  1906          0.385
3     0     1   591          0.119
4     3     1   624          0.126

ggplot(sim_coins, aes(x = factor(X), y = factor(Y))) +
  geom_jitter(alpha = 0.15, color = "steelblue", width = 0.15, height = 0.15) +
  labs(
    title = "Simulated Support and Joint PMF of (X, Y)",
    x = "X (Number of Heads)",
    y = "Y (Side Bet Payout)"
  ) +
  theme_minimal()

exec_pool <- c(rep("M", 4), rep("N", 3), rep("D", 2))

draw_execs <- \() {
  selected <- sample(exec_pool, size = 3, replace = FALSE)
  x <- sum(selected == "M")
  y <- sum(selected == "N")
  return(c(x, y))
}

sim_execs <- replicate_int(N, draw_execs()[1], .as = X) %>% 
  mutate(
    Y = map_int(X, \(x) draw_execs()[2])
  )

sim_execs_summary <- sim_execs %>% 
  summarize(cnt = n(), .by = c(X, Y)) %>% 
  mutate(
    joint_prob = cnt / N,
    prob_X_given_Y = cnt / sum(cnt), .by = Y
  )

sim_execs_table <- sim_execs_summary %>% 
  pivot_wider(
    names_from = Y,
    values_from = joint_prob,
    names_prefix = 'y=',
    id_cols = X,
    values_fill = 0
  )

print(sim_execs_table)

# A tibble: 4 × 5
      X  `y=1`  `y=2`  `y=0`  `y=3`
  <int>  <dbl>  <dbl>  <dbl>  <dbl>
1     1 0.249  0.109  0.116  0.0042
2     2 0.19   0.0773 0.0812 0.0035
3     0 0.0642 0.0253 0.0284 0.0025
4     3 0.0276 0.0098 0.0111 0.0011

p_X_given_Y2 <- sim_execs_summary %>% 
  filter(Y == 2) %>% 
  select(X, Y, cnt, prob_X_given_Y)

print(p_X_given_Y2)

# A tibble: 4 × 4
      X     Y   cnt prob_X_given_Y
  <int> <int> <int>          <dbl>
1     2     2   773         0.350 
2     1     2  1086         0.491 
3     0     2   253         0.114 
4     3     2    98         0.0443

ggplot(sim_execs, aes(x = factor(X), y = factor(Y))) +
  geom_jitter(alpha = 0.15, color = "darkgreen", width = 0.15, height = 0.15) +
  labs(
    title = "Simulated Executive Selection (X = Married, Y = Never Married)",
    x = "X (Married)",
    y = "Y (Never Married)"
  ) +
  theme_minimal()