Problema 2

La edad de una antigua pieza de materia orgánica se puede estimar a partir de la tasa a la que emite partículas beta como resultado del decaimiento del carbono-14. Por ejemplo , si X es el número de partículas emitidas durante diez minutos por un fragmento óseo con 10000 años de antigüedad que contiene 1 g de carbono, entonces X tiene una distribución de Poisson con media λ=45.62.Un arqueólogo descubrió un pequeño fragmento óseo que contiene 1 g de carbono. Si t es la edad desconocida del hueso, en años, el arqueólogo contar{a el número X de partículas emitidas en diez minutos y calculará una edad estimada tˆ con la fórmula:

El arqueólogo no lo sabe, pero el hueso tiene exactamente 10000 años de antigüedad, por lo que X tiene una distribución de Poisson con λ=45.62.

a)Genere una muestra simulada de 10000 valores de X y sus correspondientes valores de tˆ. b)Estime la media de tˆ. c)Estime la desviación estándar de tˆ. d)Estime la probabilidad de que tˆesté a 1000 años con una edad real de 10000 años.

set.seed(123)  # para reproducibilidad

# Parámetros
lambda <- 45.62
n <- 10000

# a) Muestra simulada
X <- rpois(n, lambda)
t_hat <- (log(15.3) - log(X / 10)) / 0.0001210

# b) Media de t_hat
media_t_hat <- mean(t_hat)

# c) Desviación estándar de t_hat
desv_t_hat <- sd(t_hat)

# d) Probabilidad de que t_hat esté a ±1000 años de 10000
prob_1000 <- mean(abs(t_hat - 10000) < 1000)

# Mostrar resultados
media_t_hat
## [1] 10114.79
desv_t_hat
## [1] 1237.73
prob_1000
## [1] 0.5886

Problema 5

Una planta de manufactura produce pistones para motores de combustión interna. El diámetro de los pistones debe ser de 10 cm, pero debido a variaciones en el proceso de producción, sigue una distribución normal con media μ=10 cm y desviación estándar σ=0.02 cm. Un pistón es rechazado si su diámetro está fuera del rango de 9.97 cm a 10.03 cm.

  1. Genere una muestra simulada de 5000 pistones.
  2. Estime la proporción de pistones rechazados por no cumplir con las especificaciones.
  3. Si el costo de reprocesamiento de un pistón defectuoso es de $15 USD, estime el costo promedio de reprocesamiento por cada 5000 pistones fabricados.
  4. Genere un histograma del diámetro de los pistones y evalúe si se ajusta a la distribución normal.
  5. ¿Qué estrategias podrían implementarse para reducir el número de pistones rechazados?
# Simulación y cálculos para el Problema 5
set.seed(456)

mu <- 10        # media (cm)
sigma <- 0.02   # desviación estándar (cm)
n <- 5000       # tamaño de muestra
lim_inf <- 9.97
lim_sup <- 10.03
costo_reproc <- 15  # USD por pistón defectuoso

# a) generar muestra
diam <- rnorm(n, mean = mu, sd = sigma)

# b) proporción rechazados
rechazado <- (diam < lim_inf) | (diam > lim_sup)
prop_rechazados <- mean(rechazado)
num_rechazados <- sum(rechazado)

# c) costo promedio de re-procesamiento por cada 5000 pistones
costo_total <- num_rechazados * costo_reproc

# d) histograma y ajuste normal (código para gráficas y test)
# Histograma:
hist(diam, breaks = 40, main = "Histograma del diámetro de pistones",
     xlab = "Diámetro (cm)")
# Agregar densidad teórica normal sobre el histograma:
xseq <- seq(min(diam), max(diam), length.out = 200)
lines(xseq, dnorm(xseq, mean = mu, sd = sigma) * (diff(hist(diam, plot = FALSE)$breaks)[1] * n),
      lwd = 2)

# Test de bondad de ajuste: Kolmogorov-Smirnov frente a N(mu,sigma)
ks_res <- ks.test(diam, "pnorm", mean = mu, sd = sigma)

# Mostrar resultados numéricos
prop_rechazados
## [1] 0.1296
num_rechazados
## [1] 648
costo_total
## [1] 9720
mean(diam)
## [1] 10.00049
sd(diam)
## [1] 0.01975138
ks_res
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  diam
## D = 0.011961, p-value = 0.4717
## alternative hypothesis: two-sided

Problema 7

Una constructora está evaluando la resistencia a la compresión de cilindros de concreto a los 28 días de curado. La resistencia sigue una distribución normal con media μ=30 MPa y desviación estándar σ=5 MPa. Se requiere que el 95% de los cilindros superen los 25 MPa.

  1. Genere una muestra simulada de 10000 cilindros de concreto.
  2. Estime el porcentaje de cilindros que cumplen con la resistencia mínima de 25 MPa. c)Estime el 10o percentil (P10) de la resistencia a la compresión.
  3. Si un cilindro que no cumple cuesta $50 USD en reprocesamiento, estime el costo promedio asociado a los rechazos en 10000 cilindros.
  4. ¿Qué medidas de control de calidad sugerirías implementar para mejorar el desempeño de la mezcla de concreto?
# Problema 7: simulación y cálculos
set.seed(789)

mu <- 30      # MPa
sigma <- 5    # MPa
n <- 10000

# a) Generar muestra
res <- rnorm(n, mean = mu, sd = sigma)

# b) porcentaje que cumple >= 25 MPa
cumplen <- res >= 25
prop_cumplen <- mean(cumplen)
prop_rechazados <- 1 - prop_cumplen
num_rechazados <- sum(!cumplen)

# c) si reproc = $50 USD por cilindro defectuoso (rechazado)
costo_total <- num_rechazados * 50

# c) (alternativa) cálculo teórico de proporción con N(30,5)
prop_cumplen_teo <- 1 - pnorm(25, mean = mu, sd = sigma)  # P(X >= 25)
prop_rechazados_teo <- 1 - prop_cumplen_teo

# d) percentil 10 teórico y simulado
p10_teo <- qnorm(0.10, mean = mu, sd = sigma)
p10_sim <- quantile(res, 0.10)

# Estadísticas adicionales
media_sim <- mean(res)
sd_sim <- sd(res)

# Mostrar resultados
list(
  prop_cumplen_sim = prop_cumplen,
  prop_rechazados_sim = prop_rechazados,
  num_rechazados = num_rechazados,
  costo_total_sim = costo_total,
  prop_cumplen_teo = prop_cumplen_teo,
  prop_rechazados_teo = prop_rechazados_teo,
  p10_teo = p10_teo,
  p10_sim = p10_sim,
  media_sim = media_sim,
  sd_sim = sd_sim
)
## $prop_cumplen_sim
## [1] 0.8448
## 
## $prop_rechazados_sim
## [1] 0.1552
## 
## $num_rechazados
## [1] 1552
## 
## $costo_total_sim
## [1] 77600
## 
## $prop_cumplen_teo
## [1] 0.8413447
## 
## $prop_rechazados_teo
## [1] 0.1586553
## 
## $p10_teo
## [1] 23.59224
## 
## $p10_sim
##      10% 
## 23.55523 
## 
## $media_sim
## [1] 29.97459
## 
## $sd_sim
## [1] 4.959753

Problema 9

El caudal máximo diario de un río se modela como una distribución Gumbel con parámetros μ=500m³/s y β=100m³/s. Un puente ha sido diseñado para resistir caudales de hasta 750 m³/s. Se desea estimar la probabilidad de que el puente sea superado en los próximos 10 años.

  1. Genere 3650 simulaciones (equivalente a 10 años) del caudal máximo diario.
  2. Estime la probabilidad de que en al menos un día en los próximos 10 años el caudal supere los 750 m³/s.
  3. Si se desea reducir el riesgo de superación al 1%, ¿qué caudal máximo debería resistir el nuevo diseño del puente?
  4. Construya el histograma de los caudales simulados y analice la distribución.
  5. ¿Qué medidas adicionales de mitigación podría implementar el ingeniero civil a cargo?
# Problema 9: Análisis de caudales con distribución Gumbel
# Parámetros de la distribución Gumbel
mu <- 500      # parámetro de ubicación (m³/s)
beta <- 100    # parámetro de escala (m³/s)
capacidad <- 750  # capacidad del puente (m³/s)
n_years <- 10     # periodo de análisis

# Función de distribución acumulada Gumbel
# F(x) = exp(-exp(-(x-mu)/beta))
pgumbel <- function(x, mu, beta) {
  exp(-exp(-(x - mu) / beta))
}

# Función de densidad Gumbel
dgumbel <- function(x, mu, beta) {
  z <- (x - mu) / beta
  (1/beta) * exp(-(z + exp(-z)))
}

# Función cuantil Gumbel (inversa)
qgumbel <- function(p, mu, beta) {
  mu - beta * log(-log(p))
}

# Función para generar aleatorios Gumbel
rgumbel <- function(n, mu, beta) {
  u <- runif(n)
  qgumbel(u, mu, beta)
}

# ============================================
# a) Generar 3650 simulaciones (10 años)
# ============================================
set.seed(123)  # para reproducibilidad
n_simulaciones <- 3650
caudales_simulados <- rgumbel(n_simulaciones, mu, beta)

# Máximo caudal por año (365 días por año)
caudales_maximos_anuales <- sapply(1:10, function(i) {
  inicio <- (i-1) * 365 + 1
  fin <- i * 365
  max(caudales_simulados[inicio:fin])
})

cat("a) Caudales máximos anuales simulados (m³/s):\n")
## a) Caudales máximos anuales simulados (m³/s):
print(round(caudales_maximos_anuales, 2))
##  [1] 1242.59 1068.74  990.47 1179.62 1264.98 1106.88 1222.74 1296.03 1260.06
## [10] 1049.60
cat("\nPromedio de máximos:", round(mean(caudales_maximos_anuales), 2), "m³/s\n\n")
## 
## Promedio de máximos: 1168.17 m³/s
# ============================================
# b) Probabilidad de superar 750 m³/s en al menos un día en 10 años
# ============================================
# P(X > 750) para un día
p_supera_un_dia <- 1 - pgumbel(capacidad, mu, beta)

# P(al menos un día supere en 10 años) = 1 - P(ningún día supere)
# P(ningún día supere en n días) = (1 - p)^n
p_supera_10_years <- 1 - (1 - p_supera_un_dia)^(365 * n_years)

cat("b) Probabilidad de que el puente sea superado:\n")
## b) Probabilidad de que el puente sea superado:
cat("   P(X > 750 en un día) =", round(p_supera_un_dia, 6), "\n")
##    P(X > 750 en un día) = 0.078806
cat("   P(al menos 1 día supere en 10 años) =", round(p_supera_10_years, 4), "\n")
##    P(al menos 1 día supere en 10 años) = 1
cat("   Porcentaje:", round(p_supera_10_years * 100, 2), "%\n\n")
##    Porcentaje: 100 %
# ============================================
# c) Reducir riesgo al 1%: ¿qué capacidad se necesita?
# ============================================
riesgo_objetivo <- 0.01

# P(al menos un día supere) = 0.01
# 1 - (1 - p)^(365*10) = 0.01
# (1 - p)^3650 = 0.99
# 1 - p = 0.99^(1/3650)
# p = 1 - 0.99^(1/3650)

p_diaria_objetivo <- 1 - (1 - riesgo_objetivo)^(1/(365 * n_years))

# Para Gumbel: P(X > x) = p, entonces P(X ≤ x) = 1 - p
# x = qgumbel(1 - p, mu, beta)
capacidad_nueva <- qgumbel(1 - p_diaria_objetivo, mu, beta)

cat("c) Para reducir el riesgo al 1%:\n")
## c) Para reducir el riesgo al 1%:
cat("   Probabilidad diaria objetivo:", format(p_diaria_objetivo, scientific = TRUE), "\n")
##    Probabilidad diaria objetivo: 2.753513e-06
cat("   Nueva capacidad requerida:", round(capacidad_nueva, 2), "m³/s\n\n")
##    Nueva capacidad requerida: 1780.26 m³/s
# ============================================
# d) Histograma de los caudales simulados
# ============================================
# Crear gráfico con histograma y curva teórica
hist(caudales_simulados, 
     breaks = 50, 
     probability = TRUE,
     col = "lightblue", 
     border = "white",
     main = "Distribución de Caudales Diarios Simulados",
     xlab = "Caudal (m³/s)",
     ylab = "Densidad")

# Agregar curva teórica
x_seq <- seq(min(caudales_simulados), max(caudales_simulados), length.out = 200)
lines(x_seq, dgumbel(x_seq, mu, beta), col = "red", lwd = 2)

# Agregar línea de capacidad del puente
abline(v = capacidad, col = "darkgreen", lwd = 2, lty = 2)
legend("topright", 
       legend = c("Densidad teórica Gumbel", "Capacidad del puente (750 m³/s)"),
       col = c("red", "darkgreen"), 
       lwd = 2, 
       lty = c(1, 2))

# Estadísticas básicas
cat("d) Estadísticas de los caudales simulados:\n")
## d) Estadísticas de los caudales simulados:
cat("   Media:", round(mean(caudales_simulados), 2), "m³/s\n")
##    Media: 557.19 m³/s
cat("   Mediana:", round(median(caudales_simulados), 2), "m³/s\n")
##    Mediana: 536.13 m³/s
cat("   Desv. estándar:", round(sd(caudales_simulados), 2), "m³/s\n")
##    Desv. estándar: 127.19 m³/s
cat("   Mínimo:", round(min(caudales_simulados), 2), "m³/s\n")
##    Mínimo: 296.23 m³/s
cat("   Máximo:", round(max(caudales_simulados), 2), "m³/s\n\n")
##    Máximo: 1296.03 m³/s
# ============================================
# e) Medidas de mitigación adicionales
# ============================================
cat("e) Medidas adicionales de mitigación:\n\n")
## e) Medidas adicionales de mitigación:
cat("1. SISTEMA DE ALERTA TEMPRANA:\n")
## 1. SISTEMA DE ALERTA TEMPRANA:
cat("   - Umbral de alerta: 600 m³/s (80% de capacidad)\n")
##    - Umbral de alerta: 600 m³/s (80% de capacidad)
cat("   - Probabilidad de exceder 600 m³/s en un día:", 
    round(1 - pgumbel(600, mu, beta), 4), "\n\n")
##    - Probabilidad de exceder 600 m³/s en un día: 0.3078
cat("2. ANÁLISIS DE PERIODOS DE RETORNO:\n")
## 2. ANÁLISIS DE PERIODOS DE RETORNO:
# Periodo de retorno T: T = 1/p donde p es la probabilidad anual
periodos <- c(2, 5, 10, 25, 50, 100)
for(T in periodos) {
  # Para periodo de retorno T años, la probabilidad anual es 1/T
  # P(X ≤ x) = 1 - 1/T
  caudal_T <- qgumbel(1 - 1/T, mu, beta)
  cat("   Caudal para T =", T, "años:", round(caudal_T, 2), "m³/s\n")
}
##    Caudal para T = 2 años: 536.65 m³/s
##    Caudal para T = 5 años: 649.99 m³/s
##    Caudal para T = 10 años: 725.04 m³/s
##    Caudal para T = 25 años: 819.85 m³/s
##    Caudal para T = 50 años: 890.19 m³/s
##    Caudal para T = 100 años: 960.01 m³/s
cat("\n3. CAPACIDAD DE ALIVIADERO:\n")
## 
## 3. CAPACIDAD DE ALIVIADERO:
percentil_95 <- qgumbel(0.95, mu, beta)
cat("   Percentil 95:", round(percentil_95, 2), "m³/s\n")
##    Percentil 95: 797.02 m³/s
cat("   Capacidad adicional recomendada:", round(capacidad_nueva - capacidad, 2), "m³/s\n")
##    Capacidad adicional recomendada: 1030.26 m³/s