Tiempos de falla

Introducción

Caso de estudio inspirado en un problema originalmente propuesto por David Draper para la asignatura Applied Bayesian Statistics (STAT206) en University of California, SC.

Un investigador del Departamento de Ingeniería Electrónica y Eléctrica de una universidad necesita analizar unos datos sobre los tiempos de falla de un determinado tipo de alambre (Tipo 1). En este problema, el tiempo de falla se define como el número de veces que una máquina podría tensionar el alambre antes de romperse. Los siguientes datos corresponden a \(n = 14\) tiempos de falla de una parte del experimento: \[ 495 \ \ \ 541 \ \ \ 1461 \ \ \ 1555 \ \ \ 1603 \ \ \ 2201 \ \ \ 2750 \ \ \ 3468 \ \ \ 3516 \ \ \ 4319 \ \ \ 6622 \ \ \ 7728 \ \ \ 13159 \ \ \ 21194 \] A partir de este contexto, Su incertidumbre acerca de estos datos antes de que fueran observados es intercambiable. Por lo tanto, resulta apropiado modelar los datos como condicionalmente independientes e idénticamente distribuidos. El modelo más simple para los datos del tiempo de falla involucra la distribución Exponencial: \[ \begin{equation} y_i \mid \lambda \stackrel{ \mbox{iid} }{ \sim } \textsf{Exp} ( \lambda )\text{,} \ \ \ \textrm{i.e.,} \ \ \ p( y_i \mid \lambda ) = \frac{ 1 }{ \lambda }\, \exp{\left( - \frac{y_i }{ \lambda }\right) } \quad\text{para $y_i > 0$ y $\lambda > 0$, con $i=1,\ldots,n$.} \end{equation} \]

Estadístico suficiente

Muestre que \(s=\sum_{i=1}^n y_i\) es un estadístico suficiente para \(\lambda\).

Solución:

La distribución muestral conjunta es \[ p(\boldsymbol{y}\mid\lambda) = \prod_{i=1}^n \frac{ 1 }{ \lambda }\, \exp{\left( - \frac{y_i }{ \lambda }\right)} = \lambda^{-n} \exp{\left( - \frac{s }{ \lambda }\right)}\,, \] donde \(\boldsymbol{y} = (y_1,\ldots,y_n)\) y \(s = \sum_{i=1}^n y_i\).

Dado que \(p(\boldsymbol{y}\mid\lambda)\) se puede escribir como \(p(\boldsymbol{y}\mid\lambda) = h(\boldsymbol{y})\cdot g(\theta, s)\), donde \[ h(\boldsymbol{y}) = 1 \qquad\text{y}\qquad g(\theta, s) = \lambda^{-n} \exp{\left( - \frac{s }{ \lambda }\right)}\,, \] por el teorema de factorización de Neyman-Fisher se tiene que \(s\) es un estadístico suficiente para \(\lambda\).

# datos
y <- c(495, 541, 1461, 1555, 1603, 2201, 2750, 3468, 3516, 4319, 6622, 7728, 13159, 21194)

# tamaño de la muestra
n <- length(y)
# suma 
s <- sum(y)

Distribución Gamma-Inversa

Se dice que la variable aleatoria \(X\) tiene distribución Gamma-Inversa con parámetros \(\alpha>0\) y \(\beta>0\), si la función de densidad de \(X\) está dada por: \[ X \sim \textsf{GI} ( \alpha, \beta )\text{,} \ \ \ \textrm{i.e.,} \ \ \ p ( x ) = \frac{ \beta^\alpha }{ \Gamma ( \alpha ) } x^{ - ( \alpha + 1 ) } \exp{ \left( - \frac{ \beta }{ x } \right) } \quad\text{para $x>0$}. \] Muestre que si \(X\sim\textsf{Gamma}(\alpha,\beta)\), entonces \(\frac{1}{X}\sim \textsf{GI} ( \alpha, \beta )\).

Solución:

Utilizando el teorema de la transformación de variables aleatorias se tiene que la función de densidad de \(Y=\frac{1}{X}\) es \[ f_Y(y) = f_X(g^{-1}(y))\,\left| \frac{\textsf{d}}{\textsf{d}y}g^{-1}(y) \right|\,, \] donde \[ f_X(x) = \frac{ \beta^\alpha }{ \Gamma ( \alpha ) }\, x^{ \alpha - 1 } \exp{ \left( - \beta\,x \right) }\,, \quad g(x) = \frac{1}{x}\,, \quad g^{-1}(y) = \frac{1}{y}\,, \quad \frac{\textsf{d}}{\textsf{d}y}g^{-1}(y) = -\frac{1}{y^2} \,, \] para \(x > 0\) y \(y > 0\) (si \(x>0\), entonces \(\frac{1}{x}>0\)).

Por lo tanto, \[ f_Y(y) = f_X(1/y)\,\left|-1/y^2 \right|\ = \frac{ \beta^\alpha }{ \Gamma ( \alpha ) }\, \left(\frac{1}{y}\right)^{ \alpha - 1 } \exp{ \left( - \beta\,\frac{1}{y} \right) }\times\frac{1}{y^2} = \frac{ \beta^\alpha }{ \Gamma ( \alpha ) }\, y^{ -(\alpha + 1) } \exp{ \left( -\frac{\beta}{y} \right) }\,,\quad y > 0, \] lo que corresponde a la función de densidad de una variable aleatoria con distribución Gamma-Inversa con parámetros \(\alpha\) y \(\beta\).

Distribución posterior

Considere la distribución previa \(\lambda \sim \textsf{GI} ( \alpha, \beta )\) junto con la distribución muestral Exponencial. Halle la distribución posterior de \(\lambda\).

Solución:

Aplicando el Teorema de Bayes se tiene que: \[\begin{align*} p(\lambda\mid\boldsymbol{y}) &\propto p(\boldsymbol{y}\mid\lambda)\,p(\lambda)\\ &= \prod_{i=1}^n p( y_i \mid \lambda ) \,p(\lambda)\\ &= \prod_{i=1}^n \frac{ 1 }{ \lambda }\, \exp{\left( - \frac{y_i }{ \lambda }\right) } \times \frac{ \beta^\alpha }{ \Gamma ( \alpha ) } \lambda^{ - ( \alpha + 1 ) } \exp{ \left( - \frac{ \beta }{ \lambda } \right) } \\ &\propto \lambda^{-n}\, \exp{\left( - \frac{s }{ \lambda }\right) } \times \lambda^{ - ( \alpha + 1 ) } \exp{ \left( - \frac{ \beta }{ \lambda } \right) } \\ &= \lambda^{-(\alpha + n - 1)}\, \exp{\left( - \frac{\beta + s}{ \lambda }\right) } \end{align*}\] donde \(\boldsymbol{y} = (y_1,\ldots,y_n)\) y \(s=\sum_{i=1}^n y_i\).

Dado que \(p(\lambda\mid\boldsymbol{y})\) es proporcional al núcleo de una distribución Gamma Inversa con parámetros \(\alpha + n\) y \(\beta + s\), se concluye que la distribución posterior de \(\lambda\) es \(\lambda\mid\boldsymbol{y}\sim\textrm{GI}(\alpha + n,\beta + s)\).

Distribución previa y posterior

Se tiene información externa de otro experimento de acuerdo con el cual la distribución previa de \(\lambda\) debería tener una media \(\mu_0 = 4500\) y una desviación estándar \(\sigma_0 = 1800\). Haga un gráfico de las distribuciones previa y posterior en el mismo gráfico.

Solución:

La información externa disponible indica que \[ \mu_0 = \textsf{E}(\lambda) = \frac{\beta}{\alpha-1} = 4500 \qquad\text{y}\qquad \sigma^2_0 = \textsf{Var}(\lambda) = \frac{\beta^2}{(\alpha-1)^2(\alpha-2)} = 1800^2\,. \] Despejando simultáneamente para \(\alpha\) y \(\beta\), se obtiene que \(\alpha = 8.25\) y \(\beta = 32625\).

# hiperparametros
a0 <- 8.25
b0 <- 32625

# media a priori
b0/(a0-1)

## [1] 4500

# varianza a priori
sqrt(b0^2/((a0-1)^2*(a0-2)))

## [1] 1800

# distribucion Gamma-Inversa
dinvgamma <- function(x, a, b, log = F) {
  # calcula la funcion de densidad de una Gamma-Inversa
  out <- a*log(b) - lgamma(a) - (a+1)*log(x) - b/x 
  if (log == FALSE) out <- exp(out)
  return(out)
}

# grafico
par(mfrow = c(1,1), mar=c(3,3,1,1), mgp=c(1.75,.75,0))
curve(expr = dinvgamma(x, a = a0 + n, b = b0 + s), from = 1000, to = 12000, main = "", xlab = expression(lambda), ylab = "Densidad", col = 4, lwd = 3, n = 1000)
curve(expr = dinvgamma(x, a = a0, b = b0), col = 1, lwd = 3, n = 1000, add = T)
legend("topright", legend = c("Previa","Posterior"), col = c(1,4), lwd = 3, bty = "n")

Estimador de máxima versoimilitud

Halle el estimador de máxima verosimilitud (MLE, por sus siglas en inglés) de \(\lambda\).

Solución:

El MLE de \(\lambda\) es \[ \hat\lambda_{\text{MLE}} = \mathop{\mathrm{arg\,max}}_{\lambda} p(\boldsymbol{y}\mid\lambda) \] donde \(p(\boldsymbol{y}\mid\lambda)= \prod_{i=1}^n p(y_i\mid\lambda) = \lambda^{-n}\exp{\left(-\frac{1}{\lambda}\sum_{i=1}^n y_i\right)}\) es la distribución muestral conjunta de \(\boldsymbol{y}\).

Dado que maximizar \(p(\boldsymbol{y}\mid\lambda)\) es equivalente a maximizar la función de log-verosimilitud \[ \log p(\boldsymbol{y}\mid\lambda) = -n\log\lambda -\frac{1}{\lambda}\sum_{i=1}^n y_i\,, \] se tiene que el punto crítico correspondiente satisface: \[ \begin{align*} \frac{\partial}{\partial\lambda}\log p(\boldsymbol{y}\mid\lambda) = -\frac{n}{\lambda} + \frac{1}{\lambda^2}\textstyle\sum_{i=1}^n y_i = 0\,, \end{align*} \] para \(\lambda >0\). Despejando \(\lambda\) de la expresión anterior se obtiene que el punto crítico es \(\lambda = \frac{1}{n}\textstyle\sum_{i=1}^n y_i = \bar{y}\).

Finalmente, para mostrar que \(\hat\lambda_{\text{MLE}} = \bar{y}\), solo se necesita mostrar que el punto crítico \(\lambda = \bar{y}\) maximiza \(\ell(\lambda)\). Para ello, se utiliza el criterio de la segunda derivada: \[ \left[\frac{\partial^2}{\partial\lambda^2}\log p(\boldsymbol{y}\mid\lambda) \right]_{\lambda = \bar{y}} = \left[ \frac{n}{\lambda^2} - \frac{2}{\lambda^3}\,n\bar{y} \right]_{\lambda = \bar{y}} = \frac{n}{\bar{y}^2} - \frac{2}{\bar{y}^3}\,n\bar{y} = -\frac{n}{\bar{y}^2} < 0\,. \] Por lo tanto, \(\lambda = \bar{y}\) maximiza \(\ell(\lambda)\), de donde \(\hat\lambda_{\text{MLE}} = \bar{y}\).

Inferencia Bayesiana y frecuentista

En la siguiente tabla se presentan los resultados asociados tanto con la inferencia Bayesiana y como frecuentista.

Método	Estimación	CV (%)	Intervalo al 95%
Bayesiano	4854.67	22.07	(3189.59 ; 7366.21)
Frec. Asintótico	5043.71	26.73	(2401.70 ; 7685.72)
Frec. Bootstrap	5036.74	29.38	(2524.85 ; 8263.22)

Bayesiano

Inferencia Bayesiana basada en la distribución posterior.

Se simulan \(B=100000\) valores de la distribución posterior \(p(\lambda\mid\boldsymbol{y})\):

• Simular \(\lambda^{(b)}\mid\boldsymbol{y} \sim \textsf{IG}(\alpha_n,\beta_n)\), con \(\alpha_n = \alpha+n =22.25\) y \(\beta_n = \beta + \sum_{i=1}^n y_i = 103237\).

# inferencia bayesiana
ap <- a0 + n
bp <- b0 + s
# muestras iid de la posterior de lambda 
B <- 100000
set.seed(1234)
lam_pos <- 1/rgamma(n = B, shape = ap, rate = bp)
# estimacion
round(mean(lam_pos), 2)

## [1] 4854.67

# cv
round(100*sd(lam_pos)/mean(lam_pos), 2)

## [1] 22.07

# intervalo al 95%
round(quantile(x = lam_pos, probs = c(0.025,0.975)), 2)

##    2.5%   97.5% 
## 3189.59 7366.21

Frecuentista asintótico

Inferencia frecuentista basada en que \(\hat\lambda_{\text{MLE}}\approx\textsf{N}(\lambda,\hat{I}^{-1})\) siempre que \(n\rightarrow\infty\), donde \(\hat\lambda_{\text{MLE}}\) es el MLE de \(\lambda\) y \(\hat{I}\) es la información observada de Fisher, que en este caso está dada por \[ \hat{I}(\hat\lambda_{\text{MLE}}) = \left[ -\frac{\partial^2}{\partial\lambda^2}\log p(\boldsymbol{y}\mid\lambda) \right]_{\lambda = \hat\lambda_{\text{MLE}}} = -\left(-\frac{n}{\bar{y}^2}\right) = \frac{n}{\bar{y}^2}\,. \]

# inferencia frecuentista, normalidad asintotica del MLE
# informacion observada de Fisher
yb <- mean(y)
Ihat <- n/yb^2
# estimacion
round(yb, 2)

## [1] 5043.71

# cv
round(100*sqrt(1/Ihat)/yb, 2)

## [1] 26.73

# intervalo al 95%
round(yb + c(-1,1)*qnorm(0.975)*sqrt(1/Ihat), 2)

## [1] 2401.70 7685.72

Frecuentista Bootstrap

Inferencia frecuentista basada en Bootstrap no paramétrico.

Se simulan \(B=100000\) re-muestras con reemplazo de tamaño \(n=14\) de la muestra original \(\boldsymbol{y}=(y_1,\ldots,y_n)\).

# inferencia frecuentista, Bootstrap no parametrico
B <- 100000
set.seed(1234)
out <- NULL
for (b in 1:B)
  out[b] <- mean(sample(x = y, size = n, replace = T))
# estimacion
round(mean(out), 2)

## [1] 5036.74

# cv
round(100*sd(out)/mean(out), 2)

## [1] 29.38

# intervalo al 95%
round(quantile(x = out, probs = c(0.025,0.975)), 2)

##    2.5%   97.5% 
## 2524.85 8263.22

Cálculo de probabilidades

Calcule e interprete \(\textsf{Pr}(\lambda < 4000\mid\boldsymbol{y})\) y \(\textsf{Pr}(y^* < 4000\mid\boldsymbol{y})\), donde \(y^*\) es un tiempo de falla futuro.

Solución:

Se simulan \(B=100000\) valores de la distribución posterior \(p(\lambda\mid\boldsymbol{y})\) y de la distribución predictiva posterior \(p(y^*\mid\boldsymbol{y})\). Para \(b=1,\ldots,B\):

• Simular \(\lambda^{(b)}\mid\boldsymbol{y} \sim \textsf{IG}(\alpha_n,\beta_n)\), con \(\alpha_n = \alpha+n =22.25\) y \(\beta_n = \beta + \sum_{i=1}^n y_i = 103237\).
• Simular \((y^*)^{(b)}\mid\lambda^{(b)}\sim\textsf{Exp}(\lambda^{(b)})\).

Usando \(\lambda^{(1)},\ldots,\lambda^{(B)}\) y \((y^*)^{(1)},\ldots,(y^*)^{(B)}\) se puede calcular: \[ \begin{align*} \textsf{Pr}(\lambda < 4000\mid\boldsymbol{y}) &= \int_{0}^{4000} \textsf{GI}(\lambda\mid\alpha_n,\beta_n)\,\, \text{d}\lambda \\ &\approx \frac{1}{B}\sum_{b=1}^B I\left(\lambda^{(b)} < 4000\right) = 0.2135\,. \end{align*} \]

La probabilidad posterior de que la media del tiempo de falla del alambre de metal (Tipo 1) sea menor que 4000 es aproximadamente 21%.

# probalidad posterior de lambda < 4000
round(mean(lam_pos < 4000), 4)

## [1] 0.2135

Además, \[ \begin{align*} \textsf{Pr}(y^* < 4000\mid\boldsymbol{y}) &= \int_0^{4000} \int_0^\infty \textsf{Exp}(y^*\mid\lambda)\,\textsf{GI}(\lambda\mid\alpha_n,\beta_n)\,\, \text{d}\lambda\,\, \text{d}y^* \\ &\approx \frac{1}{B}\sum_{b=1}^B I\left((y^*)^{(b)} < 4000\right) = 0.5705\,, \end{align*} \] donde \(I(A)\) es la función indicadora del conjunto \(A\).

La probabilidad posterior de que un tiempo de falla futuro/no observado del alambre de metal (Tipo 1) sea menor que 4000 es aproximadamente 57%..

# probabilidad posterior de y* < 4000
ynew_pos <- rexp(n = B, rate = 1/lam_pos)
round(mean(ynew_pos < 4000), 4) 

## [1] 0.5702

Prueba de hipótesis una población

Pruebe el sistema de hipótesis \(H_0:\lambda=\lambda_0\) frente a \(H_1:\lambda\neq\lambda_0\), con \(\lambda_0=4000\).

Solución:

Asumiendo que \(\textsf{Pr}(H_0)=\textsf{Pr}(H_1)=0.5\), el factor de Bayes \(B_{01}\) asociado con el sistema de hipótesis \(H_0:\lambda=\lambda_0\) frente a \(H_1:\lambda\neq\lambda_0\) es: \[ B_{01} = \frac{p(\boldsymbol{y}\mid H_0)}{p(\boldsymbol{y}\mid H_1)}\,. \] Así, \[ \begin{align*} p(\boldsymbol{y}\mid H_0) &= \int_{0}^\infty p(\boldsymbol{y}\mid\lambda,H_0)\,p(\lambda\mid H_0)\,\text{d}\lambda \\ &= \int_{0}^\infty \lambda_0^{-n}\exp{\left(-\tfrac{s}{\lambda_0}\right)}\,\delta_{\lambda_0}(\lambda)\,\text{d}\lambda \\ &= \lambda_0^{-n}\exp{\left(-\tfrac{s}{\lambda_0}\right)} \int_{0}^\infty\delta_{\lambda_0}(\lambda)\,\text{d}\lambda \\ &= \lambda_0^{-n}\exp{\left(-\tfrac{s}{\lambda_0}\right)}\,, \end{align*} \] donde \(s = \textstyle\sum_{i=1}^n y_i\) y \(\delta_a(x)\) es la función delta de Dirac, ya que \(\int_{0}^\infty\delta_{\lambda_0}(\lambda)\,\text{d}\lambda = 1\).

Además, \[ \begin{align*} p(\boldsymbol{y}\mid H_1) &= \int_{0}^\infty p(\boldsymbol{y}\mid\lambda,H_1)\,p(\lambda\mid H_1)\,\text{d}\lambda \\ &= \int_{0}^\infty \lambda^{-n}\exp{\left(-\tfrac{s}{\lambda}\right)}\,\frac{\beta^{\alpha}}{\Gamma(\alpha)}\lambda^{-(\alpha+1)}\exp{\left(-\frac{\beta}{\lambda}\right)}\,\text{d}\lambda\\ &= \frac{\beta^{\alpha}}{\Gamma(\alpha)} \int_{0}^\infty \lambda^{-(\alpha+n+1)}\exp{\left(-\tfrac{1}{\lambda}\left(\beta+s\right)\right)}\,\text{d}\lambda\\ &=\frac{\beta^{\alpha}}{\Gamma(\alpha)}\cdot\frac{\Gamma(\alpha+n)}{\left(\beta+s\right)^{\alpha+n}}\,. \end{align*} \] Por lo tanto, \[ B_{01} = \frac{ \lambda_0^{-n}\exp{\left(-\tfrac{s}{\lambda_0}\right)} }{ \frac{\beta^{\alpha}}{\Gamma(\alpha)}\cdot\frac{\Gamma(\alpha+n)}{\left(\beta+s\right)^{\alpha+n}} } = \lambda_0^{-n}\exp{\left(-\tfrac{s}{\lambda_0}\right)}\cdot \frac{\Gamma(\alpha)\, \left(\beta+s\right)^{\alpha+n} }{\Gamma(\alpha+n)\,\beta^{\alpha}} \] y \[ \log B_{01} = -n\log \lambda_0 - \frac{s}{\lambda_0} + \log \Gamma(\alpha) + (\alpha+n)\log\left(\beta+s\right) - \log\Gamma(\alpha+n) - \alpha\log\beta\,. \] Entonces, \(\log B_{01} = 0.2454\), \(B_{01} = 1.2781\), y en consecuencia, \(B_{10} = 0.7824\). Por lo tanto, no hay evidencia sustancial en contra de \(H_0\).

# factor de bayes
l0 <- 4000
logB01 <- -n*log(l0) - s/l0 + lgamma(a0) + (a0+n)*log(b0+s) - lgamma(a0+n) - a0*log(b0)
round(logB01, 4)

## [1] 0.2454

B01 <- exp(logB01)
round(B01, 4)

## [1] 1.2781

round(1/B01, 4)

## [1] 0.7824

Prueba de hipótesis dos poblaciones

Experimentación adicional bajo las mismas condiciones con otro tipo de alambre (Tipo 2) produjo los siguientes resultados: \[ 294\ \ \ 569 \ \ \ 766 \ \ \ 1576 \ \ \ 1602 \ \ \ 2015 \ \ \ 2166 \ \ \ 3885 \ \ \ 8141 \ \ \ 10285 \] Considerando modelos independientes de la forma \(y_{i,k} \mid \lambda_k \stackrel{ \mbox{iid} }{ \sim } \textsf{Exp} ( \lambda_k )\) con \(\lambda_k \sim \textsf{GI}(\alpha_0,\beta_0)\), para \(i = 1,\ldots,n_k\) y \(k=1,2\), donde \(y_{i,k}\) es el tiempo de falla del alambre \(i\) de tipo \(k\), y \(n_k\) es el número de alambres de tipo \(k\) sometidos a experimentación (la distribución previa es la misma para ambos tipos de alambre).

Pruebe el sistema de hipótesis \(H_0:\lambda_1=\lambda_2\) frente a \(H_1:\lambda_1\neq\lambda_2\). Reporte el factor de Bayes \(B_{10}\) e interprete los resultados.

Solución:

Nuevamente, asumiendo que \(\textsf{Pr}(H_0)=\textsf{Pr}(H_1)=0.5\), el factor de Bayes \(B_{01}\) asociado con el sistema de hipótesis \(H_0:\lambda_1=\lambda_2\) frente a \(H_1:\lambda_1\neq\lambda_2\) es: \[ B_{01} = \frac{p(\boldsymbol{y}\mid H_0)}{p(\boldsymbol{y}\mid H_1)}\,. \] Así, \[ \begin{align*} p(\boldsymbol{y}_1,\boldsymbol{y}_2\mid H_0) &= \int_0^{\infty} p(\boldsymbol{y}_1,\boldsymbol{y}_2\mid\lambda, H_0)\,p(\lambda\mid H_0)\,\text{d}\lambda \\ &= \int_0^\infty \lambda^{-n_1}\exp{\left(-\tfrac{s_1}{\lambda}\right)}\,\lambda^{-n_2}\exp{\left(-\tfrac{s_2}{\lambda}\right)}\, \frac{\beta^{\alpha}}{\Gamma(\alpha)}\lambda^{-(\alpha+1)}\exp{\left(-\frac{\beta}{\lambda}\right)}\,\text{d}\lambda \\ &= \frac{\beta^{\alpha}}{\Gamma(\alpha)}\int_0^\infty \lambda^{-(\alpha+n_1+n_2+1)}\exp{\left(-\frac{\beta+s_1+s_2}{\lambda}\right)}\,\text{d}\lambda \\ &= \frac{\beta^{\alpha}}{\Gamma(\alpha)}\,\frac{\Gamma(\alpha+n_1+n_2)}{(\beta+s_1+s_2)^{\alpha+n_1+n_2}}\,, \end{align*} \] donde \(\lambda=\lambda_1=\lambda_2\) y \(s_k = \textstyle\sum_{i=1}^{n_k} y_{i,k}\), para \(k=1,2\).

Además, \[ \begin{align*} p(\boldsymbol{y}_1,\boldsymbol{y}_2\mid H_1) &= \int_0^{\infty}\int_0^{\infty} p(\boldsymbol{y}_1,\boldsymbol{y}_2\mid\lambda_1,\lambda_2, H_1)\,p(\lambda_1,\lambda_2\mid H_1)\,\text{d}\lambda_1\,\text{d}\lambda_2 \\ &= \int_0^{\infty}\int_0^{\infty} \lambda_1^{-n_1}\exp{\left(-\tfrac{s_1}{\lambda_1}\right)}\,\lambda_2^{-n_2}\exp{\left(-\tfrac{s_2}{\lambda_2}\right)}\, \frac{\beta^{\alpha}}{\Gamma(\alpha)}\lambda_1^{-(\alpha+1)}\exp{\left(-\frac{\beta}{\lambda_1}\right)}\,\frac{\beta^{\alpha}}{\Gamma(\alpha)}\lambda_2^{-(\alpha+1)}\exp{\left(-\frac{\beta}{\lambda_2}\right)} \,\text{d}\lambda_1\,\text{d}\lambda_2 \\ &= \frac{\beta^{\alpha}}{\Gamma(\alpha)}\,\frac{\beta^{\alpha}}{\Gamma(\alpha)}\int_0^{\infty} \lambda_1^{-(\alpha+n_1+1)}\exp{\left(-\tfrac{\beta+s_1}{\lambda_1}\right)}\,\text{d}\lambda_1\int_0^{\infty} \lambda_2^{-(\alpha+n_2+1)}\exp{\left(-\tfrac{\beta+s_2}{\lambda_2}\right)}\,\text{d}\lambda_2\\ &= \frac{\beta^{\alpha}}{\Gamma(\alpha)}\,\frac{\beta^{\alpha}}{\Gamma(\alpha)}\,\frac{\Gamma(\alpha+n_1)}{(\beta+s_1)^{\alpha+n_1}}\,\frac{\Gamma(\alpha+n_2)}{(\beta+s_2)^{\alpha+n_2}}\,. \end{align*} \] Por lo tanto, \[ B_{01} = \frac{ \frac{\beta^{\alpha}}{\Gamma(\alpha)}\,\frac{\Gamma(\alpha+n_1+n_2)}{(\beta+s_1+s_2)^{\alpha+n_1+n_2}} }{ \frac{\beta^{\alpha}}{\Gamma(\alpha)}\,\frac{\beta^{\alpha}}{\Gamma(\alpha)}\,\frac{\Gamma(\alpha+n_1)}{(\beta+s_1)^{\alpha+n_1}}\,\frac{\Gamma(\alpha+n_2)}{(\beta+s_2)^{\alpha+n_2}} } = \frac{\Gamma(\alpha)}{\beta^{\alpha}}\,\frac{\Gamma(\alpha+n_1+n_2)}{\Gamma(\alpha+n_1)\Gamma(\alpha+n_2)}\,\frac{(\beta+s_1)^{\alpha+n_1}(\beta+s_2)^{\alpha+n_2}}{(\beta+s_1+s_2)^{\alpha+n_1+n_2}} \] y \[ \begin{align*} \log B_{01} &= \log\Gamma(\alpha)-\alpha\log\beta + \log\Gamma(\alpha+n_1+n_2) - \log\Gamma(\alpha+n_1) - \log\Gamma(\alpha+n_2) + \\ &\qquad\qquad\qquad\qquad\qquad(\alpha+n_1)\log(\beta+s_1) + (\alpha+n_2)\log(\beta+s_2) - (\alpha+n_1+n_2)\log(\beta+s_1+s_2)\,. \end{align*} \] Entonces, \(\log B_{01} = -0.1663\), \(B_{01} = 0.8468\), y en consecuencia, \(B_{10} = 1.181\). Por lo tanto, nuevamente no hay evidencia sustancial en contra de \(H_0\).

# datos
y1 <- c(495, 541, 1461, 1555, 1603, 2201, 2750, 3468, 3516, 4319, 6622, 7728, 13159, 21194)
y2 <- c(294, 569, 766, 1576, 1602, 2015, 2166, 3885, 8141, 10285)

# tamaños de la muestra
n1 <- length(y1)
n2 <- length(y2)
# suma 
s1 <- sum(y1)
s2 <- sum(y2)

# factor de bayes
logB01 <- lgamma(a0) - a0*log(b0) + lgamma(a0+n1+n2) - lgamma(a0+n1) - lgamma(a0+n2) + (a0+n1)*log(b0+s1) + (a0+n2)*log(b0+s2) - (a0+n1+n2)*log(b0+s1+s2)
round(logB01, 4)

## [1] -0.1663

B01 <- exp(logB01)
round(B01, 4)

## [1] 0.8468

round(1/B01, 4)

## [1] 1.181

Bondad de ajuste

Verifique la idoneidad del modelo para ambos tipos de alambre empleando como estadística de prueba la media del tiempo de falla. Presente sus resultados gráficamente comparando la distribución predictiva posterior con el valor observado correspondiente. Así mismo, reporte el valor \(p\) predictivo posterior en cada caso.

Solución:

En la siguiente figura se muestra la distribución predictiva posterior del promedio junto con el valor observado en la muestra para cada tipo de alambre. Se observa que los modelos capturan efectivamente la media de los datos dado que los valores \(p\) predictivos posteriores son 0.3948 y 0.5962, respectivamente.

# distribucion posterior
ap1 <- a0 + n1
bp1 <- b0 + s1
ap2 <- a0 + n2
bp2 <- b0 + s2
# muestras iid de la posterior de lambda 
B <- 100000
set.seed(1234)
lam1_pos <- 1/rgamma(n = B, shape = ap1, rate = bp1)
lam2_pos <- 1/rgamma(n = B, shape = ap2, rate = bp2)
# distribucion predictiva posterior
trep1 <- NULL
trep2 <- NULL
set.seed(1234)
for (i in 1:B) {
  trep1[i] <- mean(rexp(n = n1, rate = 1/lam1_pos[i]))
  trep2[i] <- mean(rexp(n = n2, rate = 1/lam2_pos[i]))
}

# grafico
par(mfrow = c(1,2), mar = c(3,3,1,1), mgp = c(1.75,0.75,0))
# windows(width = 5, height = 4)
hist(x = trep1, freq = F, col = "gray90", border = "gray90", xlim = c(0,20000), ylim = c(0,0.0003), main = "Tipo 1", ylab = "Densidad", xlab = "Promedio")
abline(v = mean(y1), lwd = 2, lty = 2, col = 2)
hist(x = trep2, freq = F, col = "gray90", border = "gray90", xlim = c(0,20000), ylim = c(0,0.0003), main = "Tipo 2", ylab = "Densidad", xlab = "Promedio")
abline(v = mean(y2), lwd = 2, lty = 2, col = 2)

# valor p
round(mean(trep1 > mean(y1)), 4)

## [1] 0.3948

round(mean(trep2 > mean(y2)), 4)

## [1] 0.5962