Funciones de R para modelos continuos

Siguiendo el material de la unidad 3.4, en R cada modelo tiene cuatro funciones que se distinguen por su primera letra:

Letra Función Significado
d densidad \(f(x)\)
p acumulada \(F(x)=P(X\leq x)\)
q cuantil / percentil \(x_p\) tal que \(F(x_p)=p\)
r aleatorios simulación

Para las probabilidades usaremos:

  • \(P(X\leq a)\) = p…(a)
  • \(P(X>a)\) = p…(a, lower.tail = FALSE)
  • \(P(a<X<b)\) = p…(b) - p…(a) (en variables continuas los extremos no importan)

Convenciones de parámetros en R (según el material del curso):

Modelo Parámetros en R \(E[X]\) \(V[X]\)
Uniforme \(U(a,b)\) min = a, max = b \(\frac{a+b}{2}\) \(\frac{(b-a)^2}{12}\)
Normal \(N(\mu,\sigma^2)\) mean = μ, sd = σ \(\mu\) \(\sigma^2\)
Exponencial \(Exp(\lambda)\) rate = λ \(\frac{1}{\lambda}\) \(\frac{1}{\lambda^2}\)
Lognormal meanlog = μ, sdlog = σ \(e^{\mu+\sigma^2/2}\) \((e^{\sigma^2}-1)e^{2\mu+\sigma^2}\)
Gamma \((\alpha,\lambda)\) shape = α, rate = λ \(\frac{\alpha}{\lambda}\) \(\frac{\alpha}{\lambda^2}\)
Beta \((\alpha,\beta)\) shape1 = α, shape2 = β \(\frac{\alpha}{\alpha+\beta}\) \(\frac{\alpha\beta}{(\alpha+\beta)^2(\alpha+\beta+1)}\)
Weibull \((k,\lambda)\) shape = k, scale = λ \(\lambda\,\Gamma(1+1/k)\) \(\lambda^2\left[\Gamma(1+2/k)-\Gamma^2(1+1/k)\right]\)

En el curso la gamma se parametriza con tasa \(\lambda\) (por eso rate = λ), y la Weibull con forma y escala, cuya acumulada es \(F(x)=1-e^{-(x/\lambda)^k}\).

Problema 1 — Uniforme \(U(20,50)\)

a <- 20; b <- 50

# a) Función de densidad: f(x) = 1/(b-a) para a <= x <= b
dens <- 1 / (b - a)

# b), c), d)
p1b <- punif(30, a, b)
p1c <- punif(42, a, b, lower.tail = FALSE)
p1d <- punif(40, a, b) - punif(25, a, b)

# e) x tal que P(X <= x) = 0.80
p1e <- qunif(0.80, a, b)

# f) Esperanza y varianza
EX <- (a + b) / 2
VX <- (b - a)^2 / 12

tabla(c("a) f(x) = 1/(b-a) en [20, 50]", "b) P(X <= 30)", "c) P(X > 42)",
        "d) P(25 <= X < 40)", "e) x con P(X <= x) = 0.80", "f) E[X]", "f) V[X]"),
      c(f4(dens), f4(p1b), f4(p1c), f4(p1d), f4(p1e), f4(EX), f4(VX)))
Inciso Resultado
a) f(x) = 1/(b-a) en [20, 50] 0.0333
b) P(X <= 30) 0.3333
c) P(X > 42) 0.2667
d) P(25 <= X < 40) 0.5000
e) x con P(X <= x) = 0.80 44.0000
f) E[X] 35.0000
f) V[X] 75.0000

La densidad es \(f_X(x)=\dfrac{1}{30}\) para \(20\leq x\leq 50\) y \(0\) en otro caso. Como la densidad es constante, la probabilidad de un intervalo es su longitud sobre 30: por ejemplo, \(P(X\le 30)=10/30\).

x <- c(10, 20, 20, 50, 50, 60)
y <- c(0, 0, dens, dens, 0, 0)
plot(x, y, type = "l", col = col, lwd = 2, las = 1, ylim = c(0, 0.04),
     main = "Uniforme U(20, 50)", xlab = "x", ylab = expression(f[X](x)))
xs <- seq(25, 40, length.out = 50)
polygon(c(25, xs, 40), c(0, rep(dens, 50), 0),
        col = adjustcolor(col, alpha.f = 0.35), border = NA)
abline(v = EX, lty = 2, col = "red")
legend("topright", legend = c("P(25 ≤ X < 40)", "E[X] = 35"),
       fill = c(adjustcolor(col, 0.35), NA), border = NA,
       lty = c(NA, 2), col = c(NA, "red"), bty = "n")

Problema 2 — Normal \(N(100, 15^2)\)

mu <- 100; s <- 15

p2a <- pnorm(85, mu, s)
p2b <- pnorm(120, mu, s, lower.tail = FALSE)
p2c <- pnorm(115, mu, s) - pnorm(90, mu, s)
p2d <- qnorm(0.95, mu, s)

# e) Intervalo central ~95%: media ± 1.96 sigma (regla empírica: ±2 sigma)
ic <- qnorm(c(0.025, 0.975), mu, s)

tabla(c("a) P(X <= 85)", "b) P(X > 120)", "c) P(90 <= X < 115)",
        "d) x con P(X <= x) = 0.95",
        "e) Intervalo central 95% (exacto, ±1.96σ)",
        "e) Aproximación empírica (μ ± 2σ)"),
      c(f4(p2a), f4(p2b), f4(p2c), f4(p2d),
        paste0("[", f4(ic[1]), " ; ", f4(ic[2]), "]"),
        paste0("[", mu - 2*s, " ; ", mu + 2*s, "]")))
Inciso Resultado
a) P(X <= 85) 0.1587
b) P(X > 120) 0.0912
c) P(90 <= X < 115) 0.5889
d) x con P(X <= x) = 0.95 124.6728
e) Intervalo central 95% (exacto, ±1.96σ) [70.6005 ; 129.3995]
e) Aproximación empírica (μ ± 2σ) [70 ; 130]

Por la propiedad empírica, aproximadamente el 95 % de los valores está en \(\mu\pm2\sigma=(70,\,130)\); con el cuantil exacto \(z_{0.975}=1.96\) el intervalo es \((70.60,\,129.40)\).

graficar(function(x) dnorm(x, mu, s), 40, 160, "Normal N(100, 15²)",
         sombra = ic, vline = mu, etiqueta = "μ = 100")

Problema 3 — Normal \(N(72, 8^2)\)

mu <- 72; s <- 8

p3a <- pnorm(60, mu, s)
p3b <- pnorm(80, mu, s, lower.tail = FALSE)
p3c <- pnorm(76, mu, s) - pnorm(68, mu, s)
p3d <- qnorm(0.90, mu, s)                       # percentil 90
p3e <- qnorm(0.05, mu, s, lower.tail = FALSE)   # deja 5% por encima

tabla(c("a) P(X < 60)", "b) P(X >= 80)", "c) P(68 <= X <= 76)",
        "d) Percentil 90", "e) x que deja 5% por encima"),
      c(f4(p3a), f4(p3b), f4(p3c), f4(p3d), f4(p3e)))
Inciso Resultado
a) P(X < 60) 0.0668
b) P(X >= 80) 0.1587
c) P(68 <= X <= 76) 0.3829
d) Percentil 90 82.2524
e) x que deja 5% por encima 85.1588

Los extremos del intervalo \((68,\,76)\) están a \(\pm4\) de la media, es decir, a \(\pm0.5\sigma\); la simetría de la normal permite calcular \(P(68\le X\le 76)=2\,\Phi(0.5)-1\).

graficar(function(x) dnorm(x, mu, s), 40, 104, "Normal N(72, 8²)",
         sombra = c(68, 76), vline = mu, etiqueta = "Media μ = 72")

Problema 4 — Normal \(N(500, 40^2)\)

mu <- 500; s <- 40

p4a <- pnorm(450, mu, s)
p4b <- pnorm(575, mu, s, lower.tail = FALSE)
p4c <- pnorm(540, mu, s) - pnorm(480, mu, s)
p4d <- qnorm(c(0.10, 0.90), mu, s)   # 80% central: deja 10% en cada cola
p4e <- qnorm(0.99, mu, s)

tabla(c("a) P(X <= 450)", "b) P(X > 575)", "c) P(480 < X < 540)",
        "d) Límites del 80% central", "e) Percentil 99"),
      c(f4(p4a), f4(p4b), f4(p4c),
        paste0("[", f4(p4d[1]), " ; ", f4(p4d[2]), "]"), f4(p4e)))
Inciso Resultado
a) P(X <= 450) 0.1056
b) P(X > 575) 0.0304
c) P(480 < X < 540) 0.5328
d) Límites del 80% central [448.7379 ; 551.2621]
e) Percentil 99 593.0539 
graficar(function(x) dnorm(x, mu, s), 350, 650, "Normal N(500, 40²)",
         sombra = p4d, vline = mu, etiqueta = "μ = 500 (área: 80% central)")

Problema 5 — Exponencial \(\lambda = 0.25\)

lam <- 0.25

p5a <- pexp(2, lam)
p5b <- pexp(5, lam, lower.tail = FALSE)
p5c <- pexp(6, lam) - pexp(2, lam)
p5d <- qexp(0.90, lam)
EX <- 1 / lam
VX <- 1 / lam^2

tabla(c("a) P(X <= 2)", "b) P(X > 5)", "c) P(2 < X <= 6)",
        "d) x con P(X <= x) = 0.90", "e) E[X]", "e) V[X]"),
      c(f4(p5a), f4(p5b), f4(p5c), f4(p5d), f4(EX), f4(VX)))
Inciso Resultado
a) P(X <= 2) 0.3935
b) P(X > 5) 0.2865
c) P(2 < X <= 6) 0.3834
d) x con P(X <= x) = 0.90 9.2103
e) E[X] 4.0000
e) V[X] 16.0000

El tiempo medio es \(E[X]=1/\lambda=4\) (en las mismas unidades de tiempo del enunciado) y la varianza \(V[X]=16\).

graficar(function(x) dexp(x, lam), 0, 25, "Exponencial (λ = 0.25)",
         sombra = c(2, 6), vline = EX, etiqueta = "E[X] = 4")

Problema 6 — Lognormal (\(\mu_{\ln}=3.5\), \(\sigma_{\ln}=0.40\))

Aquí \(\ln(X)\sim N(3.5,\,0.40^2)\), por lo que en R se usan meanlog = 3.5 y sdlog = 0.40.

ml <- 3.5; sl <- 0.40

p6a <- plnorm(25, ml, sl)
p6b <- plnorm(50, ml, sl, lower.tail = FALSE)
p6c <- plnorm(40, ml, sl) - plnorm(20, ml, sl)
p6d <- qlnorm(0.50, ml, sl)                 # mediana = exp(mu)
p6e <- qlnorm(0.90, ml, sl)
EX  <- exp(ml + sl^2 / 2)
VX  <- (exp(sl^2) - 1) * exp(2 * ml + sl^2)

tabla(c("a) P(X <= 25)", "b) P(X > 50)", "c) P(20 <= X < 40)",
        "d) Mediana de X (= e^3.5)", "e) Percentil 90",
        "f) E[X]", "f) V[X]"),
      c(f4(p6a), f4(p6b), f4(p6c), f4(p6d), f4(p6e), f4(EX), f4(VX)))
Inciso Resultado
a) P(X <= 25) 0.2411
b) P(X > 50) 0.1515
c) P(20 <= X < 40) 0.5779
d) Mediana de X (= e^3.5) 33.1155
e) Percentil 90 55.2916
f) E[X] 35.8735
f) V[X] 223.2930

La mediana (\(\approx 33.12\)) es menor que la media (\(\approx 35.87\)), como corresponde a una distribución con asimetría positiva.

graficar(function(x) dlnorm(x, ml, sl), 0, 120, "Lognormal (3.5, 0.40²)",
         sombra = c(20, 40), vline = p6d, etiqueta = "Mediana")

Problema 7 — Gamma (\(\alpha=3\), \(\lambda=0.5\))

al <- 3; la <- 0.5

p7a <- pgamma(4, shape = al, rate = la)
p7b <- pgamma(8, shape = al, rate = la, lower.tail = FALSE)
p7c <- pgamma(7, al, la) - pgamma(3, al, la)
p7d <- qgamma(0.95, shape = al, rate = la)
EX  <- al / la
VX  <- al / la^2

tabla(c("a) P(X <= 4)", "b) P(X > 8)", "c) P(3 <= X < 7)",
        "d) Percentil 95", "e) E[X]", "e) V[X]"),
      c(f4(p7a), f4(p7b), f4(p7c), f4(p7d), f4(EX), f4(VX)))
Inciso Resultado
a) P(X <= 4) 0.3233
b) P(X > 8) 0.2381
c) P(3 <= X < 7) 0.4880
d) Percentil 95 12.5916
e) E[X] 6.0000
e) V[X] 12.0000 
graficar(function(x) dgamma(x, al, la), 0, 25, "Gamma (α = 3, λ = 0.5)",
         sombra = c(3, 7), vline = EX, etiqueta = "E[X] = 6")

Problema 8 — Beta (\(\alpha=2.5\), \(\beta=4\))

a1 <- 2.5; b1 <- 4

p8a <- pbeta(0.30, a1, b1)
p8b <- pbeta(0.60, a1, b1, lower.tail = FALSE)
p8c <- pbeta(0.50, a1, b1) - pbeta(0.20, a1, b1)
p8d <- qbeta(0.75, a1, b1)
EX  <- a1 / (a1 + b1)
VX  <- a1 * b1 / ((a1 + b1)^2 * (a1 + b1 + 1))

tabla(c("a) P(X <= 0.30)", "b) P(X > 0.60)", "c) P(0.20 <= X < 0.50)",
        "d) Percentil 75", "e) E[X]", "e) V[X]"),
      c(f4(p8a), f4(p8b), f4(p8c), f4(p8d), f4(EX), f4(VX)))
Inciso Resultado
a) P(X <= 0.30) 0.3522
b) P(X > 0.60) 0.1300
c) P(0.20 <= X < 0.50) 0.5723
d) Percentil 75 0.5088
e) E[X] 0.3846
e) V[X] 0.0316 
graficar(function(x) dbeta(x, a1, b1), 0, 1, "Beta (α = 2.5, β = 4)",
         sombra = c(0.2, 0.5), vline = EX, etiqueta = "E[X] = 0.3846")

Problema 9 — Weibull (\(k=2\), \(\lambda=1000\))

k <- 2; sc <- 1000

p9a <- pweibull(500, shape = k, scale = sc)
p9b <- pweibull(1200, k, sc, lower.tail = FALSE)
p9c <- pweibull(1000, k, sc) - pweibull(600, k, sc)
p9d <- qweibull(0.50, shape = k, scale = sc)   # mediana
p9e <- qweibull(0.90, shape = k, scale = sc)

tabla(c("a) P(X <= 500)", "b) P(X > 1200)", "c) P(600 <= X < 1000)",
        "d) Mediana de X", "e) Percentil 90"),
      c(f4(p9a), f4(p9b), f4(p9c), f4(p9d), f4(p9e)))
Inciso Resultado
a) P(X <= 500) 0.2212
b) P(X > 1200) 0.2369
c) P(600 <= X < 1000) 0.3298
d) Mediana de X 832.5546
e) Percentil 90 1517.4271 
graficar(function(x) dweibull(x, k, sc), 0, 3000, "Weibull (k = 2, λ = 1000)",
         sombra = c(600, 1000), vline = p9d, etiqueta = "Mediana")

Minicasos

Problema 10 — Validación de pedidos (Uniforme)

a) Modelo. \(X\): tiempo de validación (segundos), \(X\sim U(1.5,\,4.5)\), con \(a=1.5\), \(b=4.5\) y \(f_X(x)=1/3\).

a <- 1.5; b <- 4.5

p10b <- punif(3, a, b)
p10c <- punif(4, a, b) - punif(2, a, b)
p10d <- qunif(0.90, a, b)     # inicio del 10% más lento
cumple <- p10b > 0.45

tabla(c("b) P(X < 3)", "c) P(2 < X < 4)",
        "d) Tiempo desde el cual está el 10% más lento (percentil 90)",
        "e) ¿P(X < 3) > 0.45?"),
      c(f4(p10b), f4(p10c), f4(p10d), ifelse(cumple, "Sí", "No")))
Inciso Resultado
b) P(X < 3) 0.5000
c) P(2 < X < 4) 0.6667
d) Tiempo desde el cual está el 10% más lento (percentil 90) 4.2000
e) ¿P(X < 3) > 0.45? Sí

Interpretación.

  • \(P(X<3)=\,0.5000\): la mitad de las validaciones tarda menos de 3 s.
  • Una operación es excepcionalmente lenta si dura más de 4.2000 s (percentil 90).
  • e) El criterio exige \(P(X<3)>0.45\); como 0.5000 > 0.45, sí se cumple el criterio de desempeño.
x <- c(0.5, a, a, b, b, 5.5)
y <- c(0, 0, 1/3, 1/3, 0, 0)
plot(x, y, type = "l", col = col, lwd = 2, las = 1, ylim = c(0, 0.4),
     main = "Tiempo de validación ~ U(1.5, 4.5)", xlab = "segundos",
     ylab = expression(f[X](x)))
polygon(c(a, a, 3, 3), c(0, 1/3, 1/3, 0), col = adjustcolor(col, 0.35), border = NA)
abline(v = p10d, lty = 2, col = "red")
legend("topright", legend = c("P(X < 3)", paste0("Percentil 90 = ", f4(p10d))),
       fill = c(adjustcolor(col, 0.35), NA), border = NA,
       lty = c(NA, 2), col = c(NA, "red"), bty = "n")

Problema 11 — Máquina empacadora (Normal)

a) Modelo. \(X\): peso de la bolsa (g), \(X\sim N(500,\,8^2)\).

mu <- 500; s <- 8

p11b <- pnorm(515, mu, s) - pnorm(485, mu, s)
p11c <- pnorm(520, mu, s, lower.tail = FALSE)
p11d <- qnorm(c(0.02, 0.98), mu, s)
p11e <- (1 - p11b) * 10000

tabla(c("b) P(485 <= X <= 515) (aceptación)", "c) P(X > 520)",
        "d) Pesos que delimitan el 2% inferior y el 2% superior",
        "e) Bolsas esperadas fuera de aceptación (de 10 000)"),
      c(f4(p11b), f4(p11c),
        paste0(f4(p11d[1]), " g  y  ", f4(p11d[2]), " g"), f4(p11e)))
Inciso Resultado
b) P(485 <= X <= 515) (aceptación) 0.9392
c) P(X > 520) 0.0062
d) Pesos que delimitan el 2% inferior y el 2% superior 483.5700 g y 516.4300 g
e) Bolsas esperadas fuera de aceptación (de 10 000) 607.9272 
graficar(function(x) dnorm(x, mu, s), 465, 535, "Peso de las bolsas ~ N(500, 8²)",
         sombra = c(485, 515), vline = c(485, 515),
         etiqueta = "Límites de aceptación")
abline(v = p11d, lty = 3, col = "darkgreen")

f) Análisis. Con \(\sigma=8\) el intervalo de aceptación equivale a \(\mu\pm1.875\sigma\), y el 6.0793 % de la producción (≈ 608 de cada 10 000 bolsas) queda fuera. Centrar el proceso en 500 g no garantiza por sí solo pocas unidades fuera de especificación: también importa la variabilidad. Con la misma media, una máquina con mayor \(\sigma\) tendría mucha más producción fuera de 485–515 g, como se ve a continuación:

sig <- c(4, 6, 8, 10, 12)
data.frame(sigma = sig,
           `% fuera de especificación` = round(100 * (1 - (pnorm(515, 500, sig) - pnorm(485, 500, sig))), 3),
           check.names = FALSE)

Problema 12 — Prueba estandarizada (Normal)

a) Modelo. \(X\): puntaje, \(X\sim N(70,\,10^2)\).

mu <- 70; s <- 10; n <- 600

p12b <- pnorm(85, mu, s) - pnorm(65, mu, s)
p12c <- qnorm(0.90, mu, s)     # mínimo del 10% superior
p12d <- qnorm(0.15, mu, s)     # máximo del 15% inferior
p12e <- pnorm(90, mu, s, lower.tail = FALSE)

tabla(c("b) P(65 <= X <= 85)", "c) Puntaje mínimo del 10% superior",
        "d) Puntaje máximo del 15% inferior", "e) P(X > 90)",
        "f) Estudiantes con reconocimiento (de 600)",
        "f) Estudiantes con acompañamiento (de 600)"),
      c(f4(p12b), f4(p12c), f4(p12d), f4(p12e), f4(0.10 * n), f4(0.15 * n)))
Inciso Resultado
b) P(65 <= X <= 85) 0.6247
c) Puntaje mínimo del 10% superior 82.8155
d) Puntaje máximo del 15% inferior 59.6357
e) P(X > 90) 0.0228
f) Estudiantes con reconocimiento (de 600) 60.0000
f) Estudiantes con acompañamiento (de 600) 90.0000

Se espera que 60 estudiantes reciban reconocimiento (puntaje ≥ 82.8155) y 90 sean convocados a acompañamiento (puntaje ≤ 59.6357). Son valores esperados: la cantidad observada puede variar de un grupo a otro.

graficar(function(x) dnorm(x, mu, s), 30, 110, "Puntajes ~ N(70, 10²)",
         sombra = c(65, 85), vline = c(p12d, p12c),
         etiqueta = "Cortes del 15% inferior y 10% superior")

Problema 13 — Interrupciones de un servidor (Exponencial)

a) Modelo. \(X\): horas hasta la interrupción, \(X\sim Exp(\lambda)\) con \(E[X]=200\), por lo tanto \(\lambda=1/200=0.005\) interrupciones por hora.

lam <- 1 / 200

p13b <- pexp(300, lam, lower.tail = FALSE)
p13c <- pexp(250, lam) - pexp(100, lam)
p13d <- qexp(0.20, lam)
# e) P(X > 250 | X > 150) = P(X > 100)
p13e_cond <- pexp(250, lam, lower.tail = FALSE) / pexp(150, lam, lower.tail = FALSE)
p13e_dir  <- pexp(100, lam, lower.tail = FALSE)

tabla(c("a) λ = 1/200", "b) P(X > 300)", "c) P(100 < X < 250)",
        "d) Tiempo con probabilidad acumulada 0.20",
        "e) P(X > 250 | X > 150)  (cálculo condicional)",
        "e) P(X > 100)  (propiedad sin memoria)"),
      c(f4(lam), f4(p13b), f4(p13c), f4(p13d), f4(p13e_cond), f4(p13e_dir)))
Inciso Resultado
a) λ = 1/200 0.0050
b) P(X > 300) 0.2231
c) P(100 < X < 250) 0.3200
d) Tiempo con probabilidad acumulada 0.20 44.6287
e) P(X > 250 | X > 150) (cálculo condicional) 0.6065
e) P(X > 100) (propiedad sin memoria) 0.6065

f) Propiedad. Ambos cálculos coinciden: es la propiedad de falta de memoria de la exponencial, \(P(X>s+t\mid X>s)=P(X>t)\). Haber operado 150 horas sin interrupción no cambia la probabilidad de operar 100 horas adicionales.

El mantenimiento preventivo debe programarse antes de 44.6287 horas, momento en que la probabilidad acumulada de interrupción llega al 20 %.

graficar(function(x) dexp(x, lam), 0, 1000, "Exponencial (λ = 1/200)",
         sombra = c(0, p13d), vline = p13d,
         etiqueta = "t con F(t) = 0.20", xlab = "horas")

Problema 14 — Costos de reparación (Lognormal)

a) Modelo. \(X\): costo (millones de pesos), \(X\sim\text{Lognormal}(\mu=1.2,\ \sigma=0.5)\), es decir, \(\ln X\sim N(1.2,\,0.5^2)\).

ml <- 1.2; sl <- 0.5

p14b <- plnorm(3, ml, sl)
p14c <- plnorm(5, ml, sl, lower.tail = FALSE)
p14d <- qlnorm(0.50, ml, sl)
p14e <- qlnorm(0.95, ml, sl)
media <- exp(ml + sl^2 / 2)

tabla(c("b) P(X < 3)", "c) P(X > 5)", "d) Costo mediano (millones)",
        "e) Valor superado por el 5% (percentil 95)", "f) Media E[X] (millones)"),
      c(f4(p14b), f4(p14c), f4(p14d), f4(p14e), f4(media)))
Inciso Resultado
b) P(X < 3) 0.4197
c) P(X > 5) 0.2064
d) Costo mediano (millones) 3.3201
e) Valor superado por el 5% (percentil 95) 7.5566
f) Media E[X] (millones) 3.7622

f) Media vs. mediana. La mediana es \(e^{\mu}=3.3201\) y la media es \(e^{\mu+\sigma^2/2}=3.7622\). No coinciden porque la lognormal es asimétrica positiva: pocas reparaciones muy costosas (cola derecha larga) “jalan” la media hacia arriba, mientras que la mediana, al depender sólo del orden de los datos, no se ve afectada por esos extremos. Por eso, para describir el costo típico conviene la mediana y, para fijar umbrales de costos altos, los percentiles (por ejemplo, el 5 % superior comienza en 7.5566 millones).

graficar(function(x) dlnorm(x, ml, sl), 0, 15, "Costo de reparación ~ Lognormal(1.2, 0.5²)",
         sombra = c(p14e, 15), vline = c(p14d, media),
         etiqueta = c("Mediana y media"), xlab = "millones de pesos")

Problema 15 — Soporte técnico (Gamma)

a) Modelo. \(X\): tiempo total de atención (min), \(X\sim Gamma(\alpha=4,\ \lambda=0.20)\): suma de 4 etapas exponenciales con tasa 0.20 por minuto.

al <- 4; la <- 0.20

EX   <- al / la
p15c <- pgamma(30, shape = al, rate = la)
p15d <- pgamma(25, al, la) - pgamma(15, al, la)
p15e <- qgamma(0.90, shape = al, rate = la)
meta <- p15c >= 0.80

tabla(c("b) E[X] (minutos)", "c) P(X < 30)", "d) P(15 < X < 25)",
        "e) Tiempo no superado por el 90% (percentil 90)",
        "f) ¿P(X < 30) >= 0.80?"),
      c(f4(EX), f4(p15c), f4(p15d), f4(p15e), ifelse(meta, "Sí", "No")))
Inciso Resultado
b) E[X] (minutos) 20.0000
c) P(X < 30) 0.8488
d) P(15 < X < 25) 0.3822
e) Tiempo no superado por el 90% (percentil 90) 33.4039
f) ¿P(X < 30) >= 0.80? Sí

f) Evaluación. La promesa es resolver en menos de 30 minutos y la meta es que al menos el 80 % la cumpla. Como \(P(X<30)=0.8488\) >= 0.80, la promesa sí cumple la meta de la administración. (Que el promedio sea 20 min no basta: lo que exige la promesa es la proporción acumulada hasta 30 min.)

graficar(function(x) dgamma(x, al, la), 0, 80, "Tiempo de atención ~ Gamma(4, 0.20)",
         sombra = c(0, 30), vline = 30, etiqueta = "Límite de 30 min",
         xlab = "minutos")

Problema 16 — Utilización de capacidad (Beta)

a) Modelo. \(X\): proporción de utilización, \(X\sim Beta(\alpha=5,\ \beta=2)\) en \([0,1]\).

a1 <- 5; b1 <- 2

EX   <- a1 / (a1 + b1)
p16c <- pbeta(0.80, a1, b1, lower.tail = FALSE)   # alta utilización
p16d <- pbeta(0.50, a1, b1)                       # baja utilización
p16e <- qbeta(0.90, a1, b1)                       # superado por el 10%
p16f <- 250 * p16c

tabla(c("b) E[X] (utilización promedio)", "c) P(X > 0.80) (alta)",
        "d) P(X < 0.50) (baja)", "e) Nivel superado sólo por el 10% de los días",
        "f) Días esperados con utilización > 80% (de 250)"),
      c(f4(EX), f4(p16c), f4(p16d), f4(p16e), f4(p16f)))
Inciso Resultado
b) E[X] (utilización promedio) 0.7143
c) P(X > 0.80) (alta) 0.3446
d) P(X < 0.50) (baja) 0.1094
e) Nivel superado sólo por el 10% de los días 0.9074
f) Días esperados con utilización > 80% (de 250) 86.1600

La utilización promedio es 71.4286 %. En promedio, en 86.1600 de 250 días (≈ 86 días) la utilización supera el 80 %. Aquí la probabilidad actúa como proporción esperada de días: número esperado = \(n\cdot P(X>0.80)\).

graficar(function(x) dbeta(x, a1, b1), 0, 1, "Utilización ~ Beta(5, 2)",
         sombra = c(0.80, 1), vline = EX, etiqueta = "E[X] = 0.7143")

Problema 17 — Vida útil de rodamientos (Weibull)

a) Modelo. \(X\): vida útil (horas), \(X\sim Weibull(k=2.5,\ \lambda=5000)\), con \(F(x)=1-e^{-(x/5000)^{2.5}}\).

k <- 2.5; sc <- 5000

p17b <- pweibull(3000, shape = k, scale = sc)
p17c <- pweibull(6000, k, sc, lower.tail = FALSE)
p17d <- qweibull(0.10, shape = k, scale = sc)
p17e <- qweibull(0.50, shape = k, scale = sc)
reemp <- floor(p17d / 100) * 100                 # tiempo redondeado hacia abajo
F_reemp <- pweibull(reemp, k, sc)

tabla(c("b) P(X < 3000)", "c) P(X > 6000)",
        "d) Tiempo con probabilidad acumulada de falla 0.10",
        "e) Mediana de la vida útil",
        "g) Tiempo de reemplazo propuesto (h)",
        "g) P(falla antes del reemplazo propuesto)"),
      c(f4(p17b), f4(p17c), f4(p17d), f4(p17e), reemp, f4(F_reemp)))
Inciso Resultado
b) P(X < 3000) 0.2434
c) P(X > 6000) 0.2065
d) Tiempo con probabilidad acumulada de falla 0.10 2032.5496
e) Mediana de la vida útil 4318.1745
g) Tiempo de reemplazo propuesto (h) 2000
g) P(falla antes del reemplazo propuesto) 0.0962

f) Parámetro de forma. La tasa de falla de la Weibull es \(h(t)=\dfrac{k}{\lambda}\left(\dfrac{t}{\lambda}\right)^{k-1}\). Como \(k=2.5>1\), la tasa de falla es creciente: el riesgo de falla aumenta a medida que el rodamiento envejece (desgaste). Con \(k=1\) sería constante (exponencial) y con \(k<1\) decreciente.

t <- seq(1, 10000, length.out = 500)
plot(t, (k / sc) * (t / sc)^(k - 1), type = "l", col = col, lwd = 2, las = 1,
     main = "Función de riesgo h(t), Weibull (k = 2.5, λ = 5000)",
     xlab = "horas", ylab = "h(t)")

g) Reemplazo preventivo. El criterio exige que no más del 10 % haya fallado, es decir \(F(t)\le0.10\), lo que da \(t\le 2032.5496\) horas. Una propuesta coherente es reemplazar a las 2000 horas: allí la probabilidad acumulada de falla es 0.0962 (≤ 0.10), cumpliendo el criterio.

graficar(function(x) dweibull(x, k, sc), 0, 10000, "Vida útil ~ Weibull(2.5, 5000)",
         sombra = c(0, p17d), vline = reemp,
         etiqueta = "Reemplazo propuesto", xlab = "horas")

Verificación por simulación

Como complemento (criterio de la guía: uso de herramientas computacionales y simulación), se comparan algunas probabilidades teóricas con la frecuencia relativa en 100 000 valores simulados con las funciones r….

set.seed(341)
N <- 100000

sim <- data.frame(
  Caso = c("P2: P(X > 120), N(100,15²)",
           "P5: P(X > 5), Exp(0.25)",
           "P7: P(X <= 4), Gamma(3,0.5)",
           "P9: P(X <= 500), Weibull(2,1000)",
           "P13: P(X > 300), Exp(1/200)",
           "P16: P(X > 0.8), Beta(5,2)"),
  Teorica = c(pnorm(120, 100, 15, lower.tail = FALSE),
              pexp(5, 0.25, lower.tail = FALSE),
              pgamma(4, 3, 0.5),
              pweibull(500, 2, 1000),
              pexp(300, 1/200, lower.tail = FALSE),
              pbeta(0.8, 5, 2, lower.tail = FALSE)),
  Simulada = c(mean(rnorm(N, 100, 15) > 120),
               mean(rexp(N, 0.25) > 5),
               mean(rgamma(N, 3, 0.5) <= 4),
               mean(rweibull(N, 2, 1000) <= 500),
               mean(rexp(N, 1/200) > 300),
               mean(rbeta(N, 5, 2) > 0.8))
)
sim$Teorica  <- round(sim$Teorica, 4)
sim$Simulada <- round(sim$Simulada, 4)
knitr::kable(sim)
Caso Teorica Simulada
P2: P(X > 120), N(100,15²) 0.0912 0.0931
P5: P(X > 5), Exp(0.25) 0.2865 0.2871
P7: P(X <= 4), Gamma(3,0.5) 0.3233 0.3224
P9: P(X <= 500), Weibull(2,1000) 0.2212 0.2213
P13: P(X > 300), Exp(1/200) 0.2231 0.2224
P16: P(X > 0.8), Beta(5,2) 0.3446 0.3419

Los valores simulados son muy cercanos a los teóricos, lo que respalda los cálculos realizados con las funciones d, p y q.