Inferencia Bayesiana: modelos conjugados

La inferencia bayesiana combina información previa (prior) y evidencia de los datos (likelihood) para obtener una distribución posterior sobre el parámetro de interés: \[ p(\theta\mid y)\propto p(y\mid \theta)\,p(\theta). \] Una ventaja clave es el uso de familias conjugadas, donde prior y posterior pertenecen a la misma familia, facilitando el cálculo analítico.

1) Modelo Beta–Binomial (Bernoulli/Binomial)

Modelo

Sea \(y_1,\dots,y_n\) una muestra Bernoulli con probabilidad de éxito \(\theta\): \[ y_i\mid \theta \sim \text{Bernoulli}(\theta), \qquad \text{o equivalentemente}\quad t=\sum_{i=1}^n y_i \mid \theta \sim \text{Binomial}(n,\theta). \] La verosimilitud para \(t\) es: \[ p(t\mid \theta)=\binom{n}{t}\theta^{t}(1-\theta)^{n-t}. \]

Prior conjugada

\[ \theta \sim \text{Beta}(a,b), \qquad p(\theta)\propto \theta^{a-1}(1-\theta)^{b-1}. \]

Posterior

Multiplicando verosimilitud y prior: \[ p(\theta\mid t)\propto \theta^{t+a-1}(1-\theta)^{(n-t)+b-1}, \] por lo tanto: \[ \theta\mid t \sim \text{Beta}(a+t,\; b+n-t). \] Definiendo \(a' = a+t\) y \(b' = b+n-t\), se obtienen momentos: \[ \mathbb{E}[\theta\mid t]=\frac{a'}{a'+b'}, \qquad \mathrm{Var}(\theta\mid t)=\frac{a'b'}{(a'+b')^2(a'+b'+1)}. \]

Intervalo de credibilidad

Un intervalo creíble del \(95\%\) se obtiene con cuantiles de la Beta posterior: \[ \left[q_{\;0.025},q_{\;0.975}\right] = \left[ F^{-1}_{\text{Beta}(a',b')}(0.025), F^{-1}_{\text{Beta}(a',b')}(0.975) \right]. \]

Predictiva posterior

Para una nueva observación \(y_{\text{new}}\): \[ \mathbb{P}(y_{\text{new}}=1\mid t)=\mathbb{E}[\theta\mid t]=\frac{a'}{a'+b'}. \] Para \(m\) ensayos futuros, \(T_{\text{new}}\mid t\) sigue una Beta–Binomial: \[ T_{\text{new}}\mid t \sim \text{BetaBinomial}(m,a',b'). \]

2) Modelo Gamma–Poisson (Poisson con tasa desconocida)

Modelo

Sea \(y_1,\dots,y_n\) una muestra Poisson con tasa \(\lambda\): \[ y_i\mid \lambda \sim \text{Poisson}(\lambda), \qquad p(y_i\mid \lambda)=\frac{e^{-\lambda}\lambda^{y_i}}{y_i!}. \] La verosimilitud conjunta es: \[ p(\mathbf{y}\mid \lambda)\propto e^{-n\lambda}\lambda^{\sum_{i=1}^n y_i}. \] Defina \(S=\sum_{i=1}^n y_i\).

Prior conjugada

Usando parametrización Gamma por shape–rate: \[ \lambda \sim \text{Gamma}(\alpha,\beta), \qquad p(\lambda)\propto \lambda^{\alpha-1}e^{-\beta\lambda}. \]

Posterior

\[ p(\lambda\mid \mathbf{y})\propto \left(e^{-n\lambda}\lambda^{S}\right)\left(\lambda^{\alpha-1}e^{-\beta\lambda}\right) = \lambda^{(\alpha+S)-1}e^{-(\beta+n)\lambda}, \] por lo tanto: \[ \lambda\mid \mathbf{y}\sim \text{Gamma}(\alpha+S,\;\beta+n). \] Momentos: \[ \mathbb{E}[\lambda\mid \mathbf{y}] = \frac{\alpha+S}{\beta+n}, \qquad \mathrm{Var}(\lambda\mid \mathbf{y}) = \frac{\alpha+S}{(\beta+n)^2}. \]

Predictiva posterior (Negativo Binomial)

Para una nueva cuenta \(y_{\text{new}}\), la distribución predictiva marginal es Negativo Binomial: \[ p(y_{\text{new}}\mid \mathbf{y})=\int p(y_{\text{new}}\mid \lambda)p(\lambda\mid \mathbf{y})\,d\lambda, \] y se obtiene: \[ y_{\text{new}}\mid \mathbf{y}\sim \text{NegBin}\left(r=\alpha+S,\; p=\frac{\beta+n}{\beta+n+1}\right), \] bajo la convención donde: \[ \mathbb{P}(Y=k)=\binom{k+r-1}{k}(1-p)^k p^{r},\quad k=0,1,2,\dots \] Además, \[ \mathbb{E}[y_{\text{new}}\mid \mathbf{y}] = \mathbb{E}[\lambda\mid \mathbf{y}] = \frac{\alpha+S}{\beta+n}. \]

Resumen (conjugación)

library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.4.3
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.4.3
## 
## Adjuntando el paquete: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
set.seed(123)

# Muestra
y <- rbinom(10, 1, 0.4)
n <- length(y)
t <- sum(y)

# Prior
a <- 3
b <- 1

# Posterior
ap <- a + t
bp <- b + (n - t)

# Soporte
x <- seq(0.001, 0.999, by = 0.001)
df_plot <- bind_rows(

  # Previas
  data.frame(
    x = x,
    densidad = dbeta(x, 1, 1),
    grupo = "Previa",
    etiqueta = "Beta(1,1)"
  ),
  data.frame(
    x = x,
    densidad = dbeta(x, 2, 2),
    grupo = "Previa",
    etiqueta = "Beta(2,2)"
  ),
  data.frame(
    x = x,
    densidad = dbeta(x, 4, 4),
    grupo = "Previa",
    etiqueta = "Beta(4,4)"
  ),

  # Prior
  data.frame(
    x = x,
    densidad = dbeta(x, a, b),
    grupo = "Prior vs Posterior",
    etiqueta = "Prior"
  ),

  # Posterior
  data.frame(
    x = x,
    densidad = dbeta(x, ap, bp),
    grupo = "Prior vs Posterior",
    etiqueta = "Posterior"
  )
)
ggplot(df_plot, aes(x = x, y = densidad, color = etiqueta)) +
  geom_line(linewidth = 1.1) +
  facet_wrap(~ grupo, scales = "free_y") +
  labs(
    title = "Modelo Beta-Binomial",
    x = expression(theta),
    y = "Densidad"
  ) +
  theme_minimal() +
  theme(
    legend.position = "bottom",
    strip.text = element_text(face = "bold")
  )

# Muestra
y<-rbinom(10,1,0.4)
y
##  [1] 1 0 1 0 0 1 0 0 0 1
# Previa
a<-1
b<-1
x<-seq(0.01,0.99,by=0.01)
betax<-dbeta(x,a,b)
plot(x,betax,type="l",col=3)

a<-2
b<-2
x<-seq(0.01,0.99,by=0.01)
betax<-dbeta(x,a,b)
plot(x,betax,type="l",col=3)

a<-4
b<-4
x<-seq(0.01,0.99,by=0.01)
betax<-dbeta(x,a,b)
plot(x,betax,type="l",col=3)

#Previa
# Historicamente 3 exitos por cada un fracaso
# Hiperparámetros
a<-3
b<-1
n<-length(y)
n
## [1] 10
t<-sum(y)
t
## [1] 4
# Métdo de Monte Carlo
theta_p<-rbeta(100000,a+t,b+n-t)
mean(theta_p)
## [1] 0.5003136
(a+t)/(a+b+n)
## [1] 0.5
ap<-a+t
bp<-b+n-t
var(theta_p)
## [1] 0.01665372
(ap*bp)/((ap+bp)**2*(ap+bp+1))
## [1] 0.01666667
#Intervalo de credibilidad
c(qbeta(0.025,ap,bp),qbeta(0.975,ap,bp))
## [1] 0.2513455 0.7486545
# Probailidad de cofnianza del 95%

Modelo Poisson

# =========================================
# =========================================
# MODELO GAMMA - EJEMPLO COMPLETO (CORREGIDO)
# =========================================

# Librerias
library(ggplot2)
library(dplyr)
library(MASS)   # fitdistr()
## Warning: package 'MASS' was built under R version 4.4.3
## 
## Adjuntando el paquete: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
set.seed(123)

# -------------------------
# 1) Datos simulados (verdadero)
# -------------------------
shape_true <- 3.0
rate_true  <- 0.8
scale_true <- 1 / rate_true

n_sim <- 800
x_sim <- rgamma(n_sim, shape = shape_true, rate = rate_true)

# -------------------------
# 2) Datos "historicos" (ejemplo)
# -------------------------
n_hist <- 500
x_hist <- rgamma(n_hist, shape = 2.4, rate = 0.9)

df_raw <- bind_rows(
  data.frame(x = x_sim,  fuente = "Simulado"),
  data.frame(x = x_hist, fuente = "Historico")
)

# -------------------------
# 3) Ajuste Gamma por MLE (historico)
# -------------------------
fit <- fitdistr(x_hist, densfun = "gamma")
## Warning in densfun(x, parm[1], parm[2], ...): Se han producido NaNs
shape_hat <- fit$estimate["shape"]
rate_hat  <- fit$estimate["rate"]
scale_hat <- 1 / rate_hat

# -------------------------
# 4) Curvas teoricas
# -------------------------
x_grid <- seq(0, quantile(df_raw$x, 0.995), length.out = 800)

df_curve <- bind_rows(
  data.frame(
    x = x_grid,
    densidad = dgamma(x_grid, shape = shape_true, rate = rate_true),
    curva = "Densidad verdadera (simulada)"
  ),
  data.frame(
    x = x_grid,
    densidad = dgamma(x_grid, shape = shape_hat, rate = rate_hat),
    curva = "Densidad ajustada (MLE)"
  )
)

# -------------------------
# 5) Grafico 1: Histograma + densidades
# -------------------------
p1 <- ggplot(df_raw, aes(x = x)) +
  geom_histogram(aes(y = after_stat(density)),
                 bins = 40,
                 fill = "gray85",
                 color = "gray40") +
  geom_line(data = df_curve,
            aes(x = x, y = densidad, linetype = curva),
            linewidth = 1.1) +
  facet_wrap(~ fuente, scales = "free_y") +
  labs(
    title = "Modelo Gamma: Histograma y densidad teorica",
    x = "x",
    y = "Densidad"
  ) +
  theme_minimal() +
  theme(legend.position = "bottom")

print(p1)

# -------------------------
# 6) Grafico 2: ECDF vs CDF Gamma (historico)
# -------------------------
df_ecdf <- data.frame(x = sort(x_hist)) %>%
  mutate(
    ecdf = ecdf(x_hist)(x),
    cdf_fit = pgamma(x, shape = shape_hat, rate = rate_hat)
  )

p2 <- ggplot(df_ecdf, aes(x = x)) +
  geom_line(aes(y = ecdf), linewidth = 1.1) +
  geom_line(aes(y = cdf_fit),
            linetype = "dashed",
            linewidth = 1.1) +
  labs(
    title = "Historico: ECDF vs CDF Gamma ajustada",
    x = "x",
    y = "Probabilidad acumulada"
  ) +
  theme_minimal()

print(p2)

# -------------------------
# 7) Grafico 3: Q-Q plot Gamma (historico)
# -------------------------
p <- ppoints(length(x_hist))
q_teo <- qgamma(p, shape = shape_hat, rate = rate_hat)
q_emp <- sort(x_hist)

df_qq <- data.frame(
  q_teo = q_teo,
  q_emp = q_emp
)

p3 <- ggplot(df_qq, aes(x = q_teo, y = q_emp)) +
  geom_point(alpha = 0.7) +
  geom_abline(linetype = "dashed", linewidth = 1.1) +
  labs(
    title = "Q-Q plot: Historico vs Gamma ajustada",
    x = "Cuantiles teoricos (Gamma ajustada)",
    y = "Cuantiles empiricos (historico)"
  ) +
  theme_minimal()

print(p3)

# -------------------------
# 8) Simulacion predictiva desde Gamma ajustada
# -------------------------
n_pred <- 500
x_pred <- rgamma(n_pred, shape = shape_hat, rate = rate_hat)

df_pred <- data.frame(x = x_pred)

p4 <- ggplot(df_pred, aes(x = x)) +
  geom_histogram(aes(y = after_stat(density)),
                 bins = 40,
                 fill = "gray85",
                 color = "gray40") +
  geom_line(
    data = data.frame(
      x = x_grid,
      densidad = dgamma(x_grid, shape = shape_hat, rate = rate_hat)
    ),
    aes(x = x, y = densidad),
    linewidth = 1.1
  ) +
  labs(
    title = "Predictivo: nuevas observaciones (Gamma ajustada)",
    x = "x",
    y = "Densidad"
  ) +
  theme_minimal()

print(p4)

# -------------------------
# 9) Resumen teorico Gamma
# -------------------------
cat("\n--- Ley Gamma ---\n")
## 
## --- Ley Gamma ---
cat("PDF: f(x) = (rate^shape / Gamma(shape)) * x^(shape-1) * exp(-rate*x)\n")
## PDF: f(x) = (rate^shape / Gamma(shape)) * x^(shape-1) * exp(-rate*x)
cat("Soporte: x > 0\n")
## Soporte: x > 0
cat("Media: shape / rate\n")
## Media: shape / rate
cat("Varianza: shape / rate^2\n")
## Varianza: shape / rate^2
cat("\n--- Parametros ajustados (MLE, historico) ---\n")
## 
## --- Parametros ajustados (MLE, historico) ---
cat("shape =", shape_hat, "\n")
## shape = 2.743421
cat("rate  =", rate_hat, "\n")
## rate  = 1.004718
cat("media =", shape_hat / rate_hat, "\n")
## media = 2.730538
cat("var   =", shape_hat / (rate_hat^2), "\n")
## var   = 2.717715
rm(list = ls())
graphics.off()

# ---- Paquetes ----
# install.packages(c("ggplot2","dplyr","tidyr","scales"))
library(ggplot2)
library(dplyr)
library(tidyr)
library(scales)
## Warning: package 'scales' was built under R version 4.4.3
set.seed(20260118)

# ---- 1) Datos (simulados) ----
nA <- 200
nB <- 100

# Parametros en escala log
muA <- log(7.5) - 0.5*(0.55^2)  # A mean ~ 7.5
sdA <- 0.55

muB <- log(9.5) - 0.5*(0.70^2)  # B mean ~ 9.5
sdB <- 0.70

tA <- rlnorm(nA, meanlog = muA, sdlog = sdA)
tB <- rlnorm(nB, meanlog = muB, sdlog = sdB)

df <- tibble(
  grupo = c(rep("A_NetCo", nA), rep("B_Partner", nB)),
  horas = c(tA, tB)
)

# ---- 2) Descriptivo ----
desc <- df %>%
  group_by(grupo) %>%
  summarise(
    n = n(),
    mean = mean(horas),
    median = median(horas),
    sd = sd(horas),
    p90 = quantile(horas, 0.90),
    p95 = quantile(horas, 0.95),
    .groups = "drop"
  )

print(desc)
## # A tibble: 2 × 7
##   grupo         n  mean median    sd   p90   p95
##   <chr>     <int> <dbl>  <dbl> <dbl> <dbl> <dbl>
## 1 A_NetCo     200  7.29   6.44  4.25  11.7  15.6
## 2 B_Partner   100  9.09   7.56  6.37  16.0  18.7
# ---- 3) Modelo Bayesiano Lognormal por grupo ----
# y = log(horas) ~ Normal(mu, sigma^2)
# Priors conjugados: Normal-Inverse-Gamma
posterior_normal_NIG <- function(y, m0=0, k0=0.01, a0=2, b0=2, S=50000) {
  n <- length(y)
  ybar <- mean(y)
  ss <- sum((y - ybar)^2)

  kn <- k0 + n
  mn <- (k0*m0 + n*ybar)/kn
  an <- a0 + n/2
  bn <- b0 + 0.5*ss + (k0*n*(ybar - m0)^2)/(2*kn)

  sig2 <- 1 / rgamma(S, shape = an, rate = bn)   # Inv-Gamma
  mu <- rnorm(S, mean = mn, sd = sqrt(sig2/kn))

  list(mu=mu, sig2=sig2)
}

yA <- log(df$horas[df$grupo == "A_NetCo"])
yB <- log(df$horas[df$grupo == "B_Partner"])

m0 <- log(8)     # prior debil centrado en 8h
k0 <- 0.01
a0 <- 2
b0 <- 2
S  <- 60000

postA <- posterior_normal_NIG(yA, m0=m0, k0=k0, a0=a0, b0=b0, S=S)
postB <- posterior_normal_NIG(yB, m0=m0, k0=k0, a0=a0, b0=b0, S=S)

muA_s <- postA$mu; s2A_s <- postA$sig2
muB_s <- postB$mu; s2B_s <- postB$sig2

# ---- 4) Metricas en horas (escala original) ----
# E[T]=exp(mu+0.5*sigma^2), median=exp(mu)
meanA_s <- exp(muA_s + 0.5*s2A_s)
meanB_s <- exp(muB_s + 0.5*s2B_s)

medianA_s <- exp(muA_s)
medianB_s <- exp(muB_s)

diff_mean_s  <- meanB_s - meanA_s
ratio_mean_s <- meanB_s / meanA_s

summ <- function(x) c(
  mean = mean(x),
  median = median(x),
  sd = sd(x),
  q025 = unname(quantile(x, 0.025)),
  q975 = unname(quantile(x, 0.975))
)

# Probabilidades (decision-friendly)
p_B_gt_A <- mean(diff_mean_s > 0)
p_diff_gt_1 <- mean(diff_mean_s > 1)
p_diff_gt_2 <- mean(diff_mean_s > 2)
p_ratio_gt_1_2 <- mean(ratio_mean_s > 1.2)

rope <- 0.5
p_in_rope <- mean(abs(diff_mean_s) <= rope)

cat("\n--- Posterior mean(T) in hours ---\n")
## 
## --- Posterior mean(T) in hours ---
print(rbind(A = summ(meanA_s), B = summ(meanB_s)))
##       mean   median        sd     q025      q975
## A 7.372684 7.361959 0.3109553 6.797382  8.016328
## B 9.219177 9.178980 0.6479222 8.065194 10.613867
cat("\n--- Posterior median(T) in hours ---\n")
## 
## --- Posterior median(T) in hours ---
print(rbind(A = summ(medianA_s), B = summ(medianB_s)))
##       mean   median        sd     q025     q975
## A 6.316006 6.312134 0.2473972 5.841867 6.816716
## B 7.526575 7.513791 0.4776428 6.630175 8.508300
cat("\n--- Posterior diff of means (B - A) in hours ---\n")
## 
## --- Posterior diff of means (B - A) in hours ---
print(summ(diff_mean_s))
##      mean    median        sd      q025      q975 
## 1.8464925 1.8095296 0.7201659 0.5189098 3.3690365
cat("\n--- Posterior ratio of means (B / A) ---\n")
## 
## --- Posterior ratio of means (B / A) ---
print(summ(ratio_mean_s))
##      mean    median        sd      q025      q975 
## 1.2526861 1.2460794 0.1029577 1.0675199 1.4734879
cat("\n--- Key probabilities ---\n")
## 
## --- Key probabilities ---
cat(sprintf("P(B worse than A in mean)       = %.4f\n", p_B_gt_A))
## P(B worse than A in mean)       = 0.9975
cat(sprintf("P(diff > 1 hour)                = %.4f\n", p_diff_gt_1))
## P(diff > 1 hour)                = 0.8875
cat(sprintf("P(diff > 2 hours)               = %.4f\n", p_diff_gt_2))
## P(diff > 2 hours)               = 0.3961
cat(sprintf("P(ratio > 1.2) (>=20%% worse)     = %.4f\n", p_ratio_gt_1_2))
## P(ratio > 1.2) (>=20% worse)     = 0.6845
cat(sprintf("P(|diff| <= %.1f h) (ROPE)       = %.4f\n", rope, p_in_rope))
## P(|diff| <= 0.5 h) (ROPE)       = 0.0231
# ---- 5) Graficas de datos (SIN tildes para evitar utf8towcs) ----
p1 <- ggplot(df, aes(x=horas, fill=grupo)) +
  geom_histogram(alpha=0.5, bins=30, position="identity") +
  labs(title="Repair times (hours) - Histogram", x="Hours", y="Count") +
  theme_minimal()

p2 <- ggplot(df, aes(x=horas, color=grupo)) +
  geom_density(linewidth=1) +
  labs(title="Repair times (hours) - Density", x="Hours", y="Density") +
  theme_minimal()

p3 <- ggplot(df, aes(x=grupo, y=horas, fill=grupo)) +
  geom_violin(alpha=0.5, trim=FALSE) +
  geom_boxplot(width=0.15, outlier.alpha=0.3) +
  labs(title="Distribution by group (violin + box)", x="", y="Hours") +
  theme_minimal() +
  theme(legend.position="none")

p4 <- ggplot(df, aes(x=horas, color=grupo)) +
  stat_ecdf(linewidth=1) +
  labs(title="ECDF (empirical CDF) by group", x="Hours", y="F(h)") +
  theme_minimal()

print(p1); print(p2); print(p3); print(p4)

# ---- 6) Graficas del posterior ----
post_df <- tibble(diff_mean = diff_mean_s, ratio_mean = ratio_mean_s)

p5 <- ggplot(post_df, aes(x=diff_mean)) +
  geom_density(linewidth=1) +
  geom_vline(xintercept=0, linetype="dashed") +
  geom_vline(xintercept=c(-rope, rope), linetype="dotted") +
  labs(
    title="Posterior: difference of means (B - A)",
    subtitle=sprintf("P(B>A)=%.3f | P(|diff|<=%.1fh)=%.3f", p_B_gt_A, rope, p_in_rope),
    x="Difference (hours)", y="Density"
  ) +
  theme_minimal()

p6 <- ggplot(post_df, aes(x=ratio_mean)) +
  geom_density(linewidth=1) +
  geom_vline(xintercept=1, linetype="dashed") +
  geom_vline(xintercept=1.2, linetype="dotted") +
  labs(
    title="Posterior: ratio of means (B / A)",
    subtitle=sprintf("P(ratio>1)=%.3f | P(ratio>1.2)=%.3f", mean(ratio_mean_s>1), p_ratio_gt_1_2),
    x="Ratio", y="Density"
  ) +
  theme_minimal()

print(p5); print(p6)

# ---- 7) Posterior Predictive Check (PPC) on mean ----
S_ppc <- 8000
idx <- sample(seq_along(muA_s), S_ppc)

sigmaA <- sqrt(s2A_s[idx])
sigmaB <- sqrt(s2B_s[idx])

TrepA_mean <- numeric(S_ppc)
TrepB_mean <- numeric(S_ppc)

for(i in 1:S_ppc){
  TrepA <- rlnorm(nA, meanlog = muA_s[idx[i]], sdlog = sigmaA[i])
  TrepB <- rlnorm(nB, meanlog = muB_s[idx[i]], sdlog = sigmaB[i])
  TrepA_mean[i] <- mean(TrepA)
  TrepB_mean[i] <- mean(TrepB)
}

ppc <- tibble(
  grupo = rep(c("A_NetCo","B_Partner"), each=S_ppc),
  mean_rep = c(TrepA_mean, TrepB_mean)
)

obs_means <- df %>% group_by(grupo) %>% summarise(mean_obs = mean(horas), .groups="drop")

p7 <- ggplot(ppc, aes(x=mean_rep, fill=grupo)) +
  geom_density(alpha=0.5) +
  geom_vline(data=obs_means, aes(xintercept=mean_obs), linetype="dashed") +
  labs(title="PPC: distribution of replicated mean", x="Replicated mean (hours)", y="Density") +
  theme_minimal()

print(p7)

bayes_p_A <- mean(TrepA_mean >= obs_means$mean_obs[obs_means$grupo=="A_NetCo"])
bayes_p_B <- mean(TrepB_mean >= obs_means$mean_obs[obs_means$grupo=="B_Partner"])

cat("\n--- PPC Bayesian p-values for mean ---\n")
## 
## --- PPC Bayesian p-values for mean ---
cat(sprintf("Group A: %.3f\n", bayes_p_A))
## Group A: 0.551
cat(sprintf("Group B: %.3f\n", bayes_p_B))
## Group B: 0.520
# ---- 8) Ejemplo regla tipo regulador ----
# "Sancionar si P(diff > 2h) > 0.95"
decision <- ifelse(p_diff_gt_2 > 0.95, "STRONG evidence (intervene)", "Not enough evidence")
cat("\nDecision rule: if P(diff>2h)>0.95 => ", decision, "\n")
## 
## Decision rule: if P(diff>2h)>0.95 =>  Not enough evidence