TAREA 4
Mariana Urrutia Barragán
2025-11-26
Una variable aleatoria \(Z\) se dice que tiene una distribución Power Muth (PM) con parámetros \(\beta\) y \(\gamma\) si tiene función de densidad de probabilidad dada por:
\[ f_Z(z; \beta, \gamma) = \frac{\gamma}{\beta^\gamma} z^{\gamma - 1} \left( e^{\left(\frac{z}{\beta}\right)^\gamma} - 1 \right) \exp\left[\left(\frac{z}{\beta}\right)^\gamma - e^{\left(\frac{z}{\beta}\right)^\gamma} + 1 \right], \qquad z > 0, \]
donde \(\gamma > 0\) es el
parámetro de forma y \(\beta > 0\)
el parámetro de escala.
Lo denotamos por: \[
Z \sim \text{PM}(\beta, \gamma).
\]
Aplicación
Consideramos que un conjunto de datos corresponde a la vida útil de ruptura por tensión del 49/epoxy hebras de Kevlar que fueron sometidas a una presión sostenida constante en el 70% nivel de estrés hasta que todos hubieran fallado. Los tiempos de falla se dan en horas.
Los \(n = 49\) valores son:
1051, 1337, 1389, 1921, 1942, 2322, 3629, 4006, 4012, 4063, 4921, 5445, 5620, 5817, 5905, 5956, 6068, 6121, 6473, 7501, 7886, 8108, 8546, 8666, 8831, 9106, 9711, 9806, 10205, 10396, 10861, 11026, 11214, 11362, 11604, 11608, 11745, 11762, 11895, 12044, 13520, 13670, 14110, 14496, 15395, 16179, 17092, 17568, 17568.
En el Cuadro 1 se muestra la estadística descriptiva para la vida útil de la rotura por fatiga donde \(b_1\) y \(b_2\) son los coeficientes de asimetría y curtosis de la muestra, respectivamente.
Estadísticas Descriptivas
rm(list=ls()); set.seed(123)
tiempos <- c(1051, 1337, 1389, 1921, 1942, 2322, 3629, 4006, 4012, 4063,
4921, 5445, 5620, 5817, 5905, 5956, 6068, 6121, 6473, 7501,
7886, 8108, 8546, 8666, 8831, 9106, 9711, 9806, 10205, 10396,
10861, 11026, 11214, 11362, 11604, 11608, 11745, 11762, 11895,
12044, 13520, 13670, 14110, 14496, 15395, 16179, 17092, 17568, 17568)
n <- length(tiempos)
n## [1] 49# Estadistica descriptiva
summary_s <- c(n = n,
mean = mean(tiempos),
median = median(tiempos),
sd = sd(tiempos),
var = var(tiempos),
skewness = { # simple skewness (moment)
m3 <- mean((tiempos-mean(tiempos))^3)
m2 <- mean((tiempos-mean(tiempos))^2)
m3 / m2^(3/2)
},
kurtosis = { # excess kurtosis
m4 <- mean((tiempos-mean(tiempos))^4)
m2 <- mean((tiempos-mean(tiempos))^2)
m4 / m2^2 - 3
})Cuadro 1:
print(round(summary_s,4))## n mean median sd var
## 49.0000 8805.6939 8831.0000 4553.9153 20738145.0085
## skewness kurtosis
## 0.0967 -0.8276Comparamos la distribución PM con otros modelos como Gamma (G), Weibull (W), Log-Normal (LN) y Birnbaum-Saunders (BS). Calculamos el criterio de información de Akaike AIC y el criterio de información Bayesiano BIC.
Se sabe que:
\[ AIC = 2k - 2 \, \text{loglik} \]
\[ BIC = k \log(n) - 2 \, \text{loglik} \]
donde \(k\) es el número de parámetros en el modelo, \(n\) es el tamaño de la muestra y loglik es el valor maximizado de la función de verosimilitud.
Definir Densidades:
Distribución Power–Muth (PM)
con parámetros parámetros \(β>0\) y \(γ>0\)
\[ f_{\text{PM}}(x;\beta,\gamma) = \frac{\gamma}{\beta^\gamma}\, x^{\gamma-1} \left( e^{(x/\beta)^\gamma} - 1 \right) \exp\left( (x/\beta)^\gamma - \left( e^{(x/\beta)^\gamma}-1 \right) \right), \qquad x>0. \] \[ f_{\text{G}}(x;k,\theta) = \frac{1}{\Gamma(k)\,\theta^{k}} x^{k-1} \exp\!\left(-\frac{x}{\theta}\right), \qquad x>0. \]
###Distribución Gamma (G)
Parametros: \(k>0\) y \(θ>0\)
\[ f_{\text{G}}(x;k,\theta) = \frac{1}{\Gamma(k)\,\theta^{k}} x^{k-1} \exp\!\left(-\frac{x}{\theta}\right), \qquad x>0. \]
###Distribución Weibull (W) parametros \(a>0\) y \(b>0\)
\[ f_{\text{W}}(x;a,b) = \frac{a}{b} \left( \frac{x}{b} \right)^{a-1} \exp\!\left[ -\left( \frac{x}{b} \right)^{a} \right], \qquad x>0. \]
###Distribución Log-Normal (LN)
parametros \(μ y𝜎>0\) \[ f_{\text{LN}}(x;\mu,\sigma) = \frac{1}{x\,\sigma\sqrt{2\pi}} \exp\!\left( -\frac{(\ln x - \mu)^{2}}{2\sigma^{2}} \right), \qquad x>0. \]
###Distribución Birnbaum–Saunders (BS)
paremetros: forma \(α>0\) y escala \(β>0\)
\[ f_{\text{BS}}(x;\alpha,\beta) = \frac{\sqrt{\frac{x}{\beta}}+\sqrt{\frac{\beta}{x}}} {2\alpha x} \, \phi\!\left( \frac{1}{\alpha} \left( \sqrt{\frac{x}{\beta}} - \sqrt{\frac{\beta}{x}} \right) \right), \qquad x>0, \]
#Power-Muth (PM)
dpm <- function(x, beta, gamma){
z <- (x / beta)^gamma
out <- (gamma / (beta^gamma)) * x^(gamma - 1) * (exp(z) - 1) * exp( z - (exp(z) - 1) )
return(out)
}
# log
logdpm <- function(x, beta, gamma){
z <- (x / beta)^gamma
log_expz_minus1 <- ifelse(z > 50, z + log1p(-exp(-z)), log(exp(z)-1)) # approximate for large z
res <- log(gamma) - gamma*log(beta) + (gamma - 1)*log(x) + log_expz_minus1 + ( z - (exp(z)-1) )
return(res)
}
# Gamma
logdgamma_manual <- function(x, shape, scale){
res <- (shape-1)*log(x) - x/scale - lgamma(shape) - shape*log(scale)
return(res)
}
# Weibull
logdweibull_manual <- function(x, shape, scale){
res <- log(shape) - log(scale) + (shape-1)*(log(x) - log(scale)) - (x/scale)^shape
return(res)
}
# Log-Normal
logdln_manual <- function(x, mu, sigma){
res <- -log(x) - log(sigma) - 0.5*log(2*pi) - ( (log(x)-mu)^2 ) / (2*sigma^2)
return(res)
}
# Birnbaum-Saunders
logdbs_manual <- function(x, alpha, beta){
s1 <- sqrt(x/beta)
s2 <- sqrt(beta/x)
arg <- (s1 - s2)/alpha
res <- log(s1 + s2) - log(2) - log(alpha) - log(x) + dnorm(arg, log=TRUE)
return(res)
}Log-verosimilitudes y optim (maximización)
nll_pm <- function(par){
tb <- par[1]
tg <- par[2]
beta <- exp(tb)
gamma <- 0.5 + exp(tg)
v <- logdpm(tiempos, beta, gamma)
val <- -sum(v)
if(!is.finite(val)) val <- 1e20
return(val)
}
init_pm <- c(log(median(tiempos)), log(1-0.5 + 1e-8))
init_pm <- c(log(median(tiempos)), log(0.5))
opt_pm <- optim(init_pm, nll_pm, method="Nelder-Mead", control=list(maxit=20000))
beta_hat_pm <- exp(opt_pm$par[1])
gamma_hat_pm <- 0.5 + exp(opt_pm$par[2])
loglik_pm <- -opt_pm$value
#Gamma MLE
nll_gamma <- function(par){
# par: log(shape), log(scale)
shape <- exp(par[1])
scale <- exp(par[2])
v <- logdgamma_manual(tiempos, shape, scale)
val <- -sum(v)
if(!is.finite(val)) val <- 1e20
return(val)
}
#MM
mom_shape <- (mean(tiempos)^2) / var(tiempos)
mom_scale <- var(tiempos) / mean(tiempos)
init_gamma <- log(c(mom_shape, mom_scale))
opt_gamma <- optim(init_gamma, nll_gamma, method="Nelder-Mead", control=list(maxit=20000))
shape_hat_gamma <- exp(opt_gamma$par[1])
scale_hat_gamma <- exp(opt_gamma$par[2])
loglik_gamma <- -opt_gamma$value
# Weibull MLE
nll_weibull <- function(par){
# par: log(shape), log(scale)
shape <- exp(par[1])
scale <- exp(par[2])
v <- logdweibull_manual(tiempos, shape, scale)
val <- -sum(v)
if(!is.finite(val)) val <- 1e20
return(val)
}
init_weibull <- log(c(1, mean(tiempos)))
opt_weibull <- optim(init_weibull, nll_weibull, method="Nelder-Mead", control=list(maxit=20000))
shape_hat_weibull <- exp(opt_weibull$par[1])
scale_hat_weibull <- exp(opt_weibull$par[2])
loglik_weibull <- -opt_weibull$value
#Lognormal MLE
nll_lognorm <- function(par){
# par: mu, log(sigma)
mu <- par[1]
sigma <- exp(par[2])
v <- logdln_manual(tiempos, mu, sigma)
val <- -sum(v)
if(!is.finite(val)) val <- 1e20
return(val)
}
init_logn <- c(mean(log(tiempos)), log(sd(log(tiempos))))
opt_logn <- optim(init_logn, nll_lognorm, method="Nelder-Mead", control=list(maxit=20000))
mu_hat_logn <- opt_logn$par[1]
sigma_hat_logn <- exp(opt_logn$par[2])
loglik_logn <- -opt_logn$value
# Birnbaum-Saunders MLE
nll_bs <- function(par){
# par: log(alpha), log(beta)
alpha <- exp(par[1])
beta <- exp(par[2])
v <- logdbs_manual(tiempos, alpha, beta)
val <- -sum(v)
if(!is.finite(val)) val <- 1e20
return(val)
}
init_bs <- log(c(0.5, median(tiempos)))
opt_bs <- optim(init_bs, nll_bs, method="Nelder-Mead", control=list(maxit=20000))
alpha_hat_bs <- exp(opt_bs$par[1])
beta_hat_bs <- exp(opt_bs$par[2])
loglik_bs <- -opt_bs$valueAIC, BIC (k = params)
k_pm <- 2; k_gamma <- 2; k_weibull <- 2; k_logn <- 2; k_bs <- 2
AIC_pm <- 2*k_pm - 2*loglik_pm
BIC_pm <- k_pm * log(n) - 2*loglik_pm
AIC_gamma <- 2*k_gamma - 2*loglik_gamma
BIC_gamma <- k_gamma * log(n) - 2*loglik_gamma
AIC_weibull <- 2*k_weibull - 2*loglik_weibull
BIC_weibull <- k_weibull * log(n) - 2*loglik_weibull
AIC_logn <- 2*k_logn - 2*loglik_logn
BIC_logn <- k_logn * log(n) - 2*loglik_logn
AIC_bs <- 2*k_bs - 2*loglik_bs
BIC_bs <- k_bs * log(n) - 2*loglik_bsCuadro Comparativo
res <- data.frame(
Modelo = c("PM(β,γ)","Gamma(shape,scale)","Weibull(shape,scale)","LogNormal(mu,sigma)","Birnbaum-Saunders(alpha,beta)"),
Param1 = c(round(beta_hat_pm,4), round(shape_hat_gamma,4), round(shape_hat_weibull,4), round(mu_hat_logn,4), round(alpha_hat_bs,4)),
Param2 = c(round(gamma_hat_pm,4), round(scale_hat_gamma,4), round(scale_hat_weibull,4), round(sigma_hat_logn,4), round(beta_hat_bs,4)),
logLik = c(round(loglik_pm,4), round(loglik_gamma,4), round(loglik_weibull,4), round(loglik_logn,4), round(loglik_bs,4)),
AIC = c(round(AIC_pm,4), round(AIC_gamma,4), round(AIC_weibull,4), round(AIC_logn,4), round(AIC_bs,4)),
BIC = c(round(BIC_pm,4), round(BIC_gamma,4), round(BIC_weibull,4), round(BIC_logn,4), round(BIC_bs,4))
)
print(res)## Modelo Param1 Param2 logLik AIC BIC
## 1 PM(β,γ) 8604.2578 0.8506 -479.9281 963.8561 967.6398
## 2 Gamma(shape,scale) 2.7792 3168.6881 -483.1364 970.2727 974.0564
## 3 Weibull(shape,scale) 2.0149 9907.6461 -480.8479 965.6959 969.4795
## 4 LogNormal(mu,sigma) 8.8928 0.7012 -487.8733 979.7466 983.5302
## 5 Birnbaum-Saunders(alpha,beta) 0.7520 6800.4018 -488.4345 980.8689 984.6526Histograma
grid_x <- seq(min(tiempos)*0.8, max(tiempos)*1.05, length.out = 400)
pdf_pm <- sapply(grid_x, function(x) dpm(x, beta_hat_pm, gamma_hat_pm))
pdf_gamma <- sapply(grid_x, function(x) exp(logdgamma_manual(x, shape_hat_gamma, scale_hat_gamma)))
pdf_weib <- sapply(grid_x, function(x) exp(logdweibull_manual(x, shape_hat_weibull, scale_hat_weibull)))
pdf_logn <- sapply(grid_x, function(x) exp(logdln_manual(x, mu_hat_logn, sigma_hat_logn)))
pdf_bs <- sapply(grid_x, function(x) exp(logdbs_manual(x, alpha_hat_bs, beta_hat_bs)))
# Reducir márgenes
par(mar = c(4, 4, 2, 1)) # bottom, left, top, right
hist(tiempos, breaks=15, prob=TRUE, xlab="Tiempos (horas)", main="Histograma y densidades ajustadas")
lines(grid_x, pdf_pm, lwd=2, lty=1)
lines(grid_x, pdf_gamma, lwd=2, lty=2)
lines(grid_x, pdf_weib, lwd=2, lty=3)
lines(grid_x, pdf_logn, lwd=2, lty=4)
lines(grid_x, pdf_bs, lwd=2, lty=5)
legend("topright", legend=c("PM","Gamma","Weibull","LogNormal","Birnbaum-Saunders"),
lty=1:5, lwd=2, bty="n")
Conclusiones
- El modelo con menor AIC Y BIC es \(PM(β,γ)\). Por tanto se concluye que el modelo \(PM(β,γ)\) se ajusta mejor a los datos.