Tiempo de falla de un alambre

Preguntas del caso de estudio:

Primera pregunta:

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

Desarrollo:

Para demostrar esta propiedad sobre la estadística \(s\) haremos uso del Criterio de factorización de Fisher-Neyman (Bickel and Doksum, 2016).

Este criterio establece que en un modelo regular, una estadística \(T(X)\) con rango \(\tau\) es suficiente para \(\theta\) si, y solo si, existe una función \(g(t, \theta)\) definida para \(t\) en \(\tau\) y \(\theta \in \Theta\) y una función \(h\) definida sobre \(\mathbb R^n\), tal que:

\[p(x; \theta) = g (T (x) ; \theta) h(x) \hspace{0.5cm} ∀x \in \mathbb R^n, \hspace{0.1cm} \theta\in\Theta\]

Para nuestro caso de estudio, \(y_i|\lambda \overset{iid}{\sim}Exp(\lambda)\) con \(y_i>0\) y \(\lambda>0\), \(i=1,...,n\) de aquí tenemos que:

\[\begin{aligned}p(\bf{y}|\lambda)&=p(y_1,...,y_n|\lambda)\\&\overset{1}{=}p(y_1|\lambda),...,p(y_n|\lambda)\\&\overset{2}{=}\prod_{i=1}^np(y_i|\lambda)\\&=\prod_{i=1}^n\frac{1}{\lambda}e^{- y_i/\lambda}*\bf{1}_{y\in(0,\infty)}\\&=\frac{1}{\lambda^n}e^{-\sum_{i=1}^n y_i/\lambda}*\bf{1}_{y\in(0,\infty)}\\&=g(T(y),\lambda)h(y)\end{aligned} \] Teniendo en cuenta:

Hipótesis de independencia condicional dado \(\lambda\).
Hipótesis de variables aleatorias idénticamente distribuidas.

Por la última igualdad se cumple la factorización de Fisher-Neyman, siendo así \(T(y)=\sum_{i=1}^n y_i\) un estadístico suficiente para\(\lambda\), \(g(T(y),\lambda)=\frac{e^{-\sum_{i=1}^n y_i/\lambda}}{\lambda^n}\) es la función que depende de la estadística suficiente y del parámetro \(\lambda\), además \(h(y)=\bf{1}_{y\in(0,\infty)}\).

Segunda pregunta:

Se dice que la variable aleatoria X tiene distribución Gamma-Inversa con parámetros \(\alpha\) y \(\beta\), si la función de densidad de \(\textbf{X}\) está dada por:

\[X \sim \mathrm{Gl}(\alpha, \beta), \quad \text { i.e., } \quad 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 Gamma(\alpha,\beta)\), entonces \(\frac{1}{X}\sim GI(\alpha,\beta)\).

Desarrollo:

Sea \(X\sim Gamma(\alpha,\beta)\) con \(x\in (0,\infty)\), \(\alpha\) y \(\beta\), considere la variable aleatoria \(Y=\frac{1}{X}\) por el teorema de la transformación (L.Blanco,2004) su función de densidad viene dada por

\[f_Y(y)=\left\lbrace \begin{aligned}&f_X(h^{-1}(y))|\frac{d}{dy}h^{-1}(y)|\hspace{1cm} si\ y=h(x)\ para\ algún\ x>0\\&0 \hspace{4.8cm} en\ otro\ caso \end{aligned} \right.\] con \(h^{-1}(.)\) la inversa de \(h(.)\).

En este caso como \(y=\frac{1}{x}=h(x)\), entonces \(h^{-1}(y)=x=\frac{1}{y}\). Note que la función de densidad dada por \(h(x)=\frac{1}{x}\) es una función estrictamente monótona y diferenciable, por lo tanto:

\[\begin{aligned}f_Y(y)&=\frac{\beta^\alpha}{\varGamma(\alpha)}\left(\frac{1}{y}\right)^{\alpha-1}e^{-\beta/y}|\frac{d}{dy}\frac{1}{y}|\\&=\frac{\beta^\alpha}{\varGamma(\alpha)}\left(\frac{1}{y}\right)^{\alpha-1}e^{-\beta/y}\frac{1}{y^2}\\&=\frac{\beta^\alpha}{\varGamma(\alpha)}\left(\frac{1}{y}\right)^{\alpha+1}e^{-\beta/y}\\&=\frac{\beta^\alpha}{\varGamma(\alpha)}\left(y\right)^{-(\alpha+1)}e^{-\beta/y}\end{aligned}\] Dado que \(x\in (0,\infty)\), luego \(y\in (0,\infty)\), Con lo que concluimos que \(Y\sim GI(\alpha,\beta)\) con \(y>0\).

Tercera pregunta:

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

Desarrollo:

Sea \(\lambda\sim GI(\alpha,\beta)\) y \(y_i|\lambda \overset{iid}{\sim}Exp(\lambda)\) para \(y_i>0\) y \(\lambda>0\) con \(i=1,...,n\). Hipótesis 1. Sabemos que las variables aleatorias del vector y son independientes e identicamente distribuidas, entonces

\[\begin{aligned}p(\lambda|\bf y)&\propto p(\bf y|\lambda)p(\lambda)\\&=p(y_1,y_2,...,y_n|\lambda)p(\lambda)\\&\overset{1}{=}p(y_1|\lambda)p(y_2|\lambda),...,p(y_n|\lambda)p(\lambda)\\&\overset{1}{=}\prod_{i=1}^np(y_i|\lambda)p(\lambda)\\&=\left(\prod_{i=1}^n\frac{1}{\lambda}e^{-y_i/\lambda}\right)\left(\frac{\beta^\alpha}{\varGamma(\alpha)}\lambda^{-(\alpha+1)}e^{-\beta/\lambda}\right)\\&=\left(\frac{1}{\lambda^n}e^{-\sum_{i=1}^ny_i/\lambda}\right)\left(\frac{\beta^\alpha}{\varGamma(\alpha)}\lambda^{-(\alpha+1)}e^{-\beta/\lambda}\right)\\&=\frac{\beta^\alpha}{\varGamma(\alpha)}\lambda^{-(\alpha+n+1)}e^{-(s+\beta)/\lambda}\end{aligned}\] Con \(s=\sum_{i=1}^ny_i\), vemos que esta última igualdad es el núcleo de una gamma inversa por lo que se tiene una familia conjugada, así

\[\lambda|\bf{y}\sim GI(\alpha+n,s+\beta)\]

Cuarta pregunta:

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.

Desarrollo:

La media de la previa es \(\frac{\beta}{\alpha-1}\) y la desviación estandar es \(\frac{\beta}{(\alpha-1)\sqrt{\alpha-2}}\) de acuerdo con la información de otro experimento se tiene el siguiente sistema de ecuaciones.

\[1.)\frac{\beta}{\alpha-1}=4500\hspace{2cm}2.)\frac{\beta}{(\alpha-1)\sqrt{\alpha-2}}=1800 \] despejando 1 en 2

\[\begin{aligned}\frac{4500}{\sqrt{\alpha-2}}&=1800\\\frac{5}{2}&=\sqrt{\alpha-2}\\\frac{25}{4}&=\alpha-2\\\frac{33}{4}&=\alpha\end{aligned}\] luego \(\beta=4500(\frac{33}{4}-1)=4500(\frac{29}{4})=32625\), además \(s=\sum_{i=1}^ny_i\) dadas las mediciones del experimento.

La previa tiene distribución gamma inversa y gracias a la información de un experimento anterior tiene hiperparámetros \(\alpha=\frac{33}{4}\) y \(\beta=32625\), por lo tanto, la posterior al ser una familia conjugada tiene distribución gamma inversa de parámetros \(\alpha_p=\alpha+n=\frac{89}{4}\) y \(\beta_p=\beta+s=32625+70612=103237\). A continuación se muestra la gráfica de la distribución previa y la distribución posterior.

Datos tomados del experimento

##    [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
## y1  495  541 1461 1555 1603 2201 2750 3468 3516  4319  6622  7728 13159 21194

cantidad de datos tomados

## [1] 14

Estadístico Suficiente \(s=\sum_{i=1}^ny_i=70612\)

s=sum(c(495,541,1461,1555,1603,2201,2750,3468,3516,4319,6622,7728,13159,21194));s

## [1] 70612

Gráfica de las distribuciones previa y posterior

Quinta pregunta:

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

Desarrollo:

Se debe maximizar la función de la muestra \(p(\bf{y}|\lambda)\) en función de \(\lambda\) y para facilidad en los cálculos se usa la función log-verosimilitud.

\[\begin{aligned}L(\lambda)&=log(\frac{e^{-s/\lambda}}{\lambda^n})\\&=-\frac{s}{\lambda}-nlog(\lambda)\end{aligned}\]

Hallamos sus puntos crítivos con el criterio de la primera derivada respecto a \(\lambda\) e igualamos a cero.

\[\frac{\partial}{\partial\lambda}L(\lambda)=\frac{s}{\lambda^2}-\frac{n}{\lambda}=0\] de donde

\[\begin{aligned}s&=\frac{n\lambda^2}{\lambda}\\\frac{s}{n}&=\lambda\end{aligned}\]

Luego \(\hat{\lambda}=\frac{s}{n}=\bar{y}\) es un punto crítico. Para que \(\hat{\lambda}\) sea un máximo, se´gun el criterio de la segunda derivada, debe cumplirse que:

\[\frac{\partial^2}{\partial\lambda^2}L(\lambda)=\frac{-2s\lambda}{\lambda^4}+\frac{n}{\lambda^2}=\frac{-2s\lambda^3+\lambda^4n}{\lambda^6}<0\]

Tenemos que:

\[\begin{aligned}\frac{-2s\lambda^3+\lambda^4n}{\lambda^6}&<0\\\frac{\lambda^4n}{\lambda^6}&<\frac{2s\lambda^3}{\lambda^6} \\ \frac{n}{\lambda^2}&<\frac{2s}{\lambda^3}\end{aligned}\] Y reemplazando \(\lambda=\hat{\lambda}=\bar{y}\)

\[\begin{aligned} \frac{n}{\bar{y}^2}&<\frac{2s}{\bar{y}^3}\\ \bar{y}&<\frac{2s}{n}\\ \bar{y}&<2\bar{y}\end{aligned}\] Como esto siempre se tiene en este caso, la media muestral \(\bar{y}\) es el estimador de máxima verosimilitud de \(\lambda\).

Sexta pregunta:

Completar la tabla:

Desarrollo:

Método Bayesiano

La estimación bayesiana es el valor esperado de la distribución posterior gamma inversa.

Parámetros de la distribución previa.

a=33/4
b=32625

Parámetros de la distribución posterior.

ap=a+n
bp=b+s

Estimación Bayesiana: \(E(\lambda|\bf{y})=\frac{\beta+s}{\alpha+n-1}\).

mediaP=round(bp/(ap-1),3);mediaP              #Media posterior

## [1] 4858.212

Desviación estándar Bayesiana

desvP=round(bp/((ap-1)*sqrt(ap-2)),3);desvP   #Desviación posterior

## [1] 1079.603

Método Frecuentista

Media frecuentista \(\bar{y}\).

mediaF=round(mean(yi),3);mediaF               #Media muestral

## [1] 5043.714

Desviación estándar frecuentista.

desvF=round(sd(yi),3);desvF                        #Desviación estandar muestral

## [1] 5770.487

Intervalo de confianza frecuentista.

intcF=round(c(mediaF-(qnorm(0.975)*desvF/sqrt(n)),mediaF+(qnorm(0.975)*desvF/sqrt(n))),3);intcF

## [1] 2021.004 8066.424

Método de Bootstrap

Se calcula un vector de 1000 medias de muestras con reemplazo de los datos muestrales dados en el ejercicio para calcular su media, desviación estandar y su intervalos de confianza.

## [1] 4521.714 5324.286 4509.429 3788.214 5513.286 3917.286

Media de la distribución Boostrap.

mediaB=round(mean(medias),3);mediaB

## [1] 5089.995

Desviación estándar de la distribución Boostrap.

desvB=round(sd(medias),3);desvB

## [1] 1488.127

Intervalo de confianza

El intervalo de confianza es basado en los percentiles de la distribución Boostrap dado a que el sesgo es muy bajo. El sesgo se calcula restando la media de la muestra y la media de la distribución Boostrap.

## [1] 46.281

La tabla pedida por el ejercicio es la siguiente:

Método	Estimación	CV	L. Inf.	L. Sup.
Bayesiano	4858.212	0.222	3186.027	7381.964
Frec. Asintótico	5043.714	1.144	2021.004	8066.424
Frec. Boostrap	5089.995	3.42	2465.225	8402.573

Séptima pregunta:

Calcule e interprete \(Pr(\lambda<4000|\bf{y})\) y \(Pr(y^{*}<4000|\bf{y})\)

Desarrollo:

Para hallar \(Pr(\lambda<4000|\bf{y})\) se simulan de la distribución posterior 1 millón de datos, cuentan aquellos que son menores que 4000 se divide entre los datos generados para obtener la probabilidad.

set.seed(1)
prob=rinvgamma(n = 1000000,ap,bp)
round(sum(prob<4000)/1000000,3)

## [1] 0.215

Por lo tanto \(Pr(\lambda<4000|\bf{y})=0.215\).

Para hallar la \(Pr(y^*<4000|\bf{y})\) debemos calcular la función de distribución como sigue.

\[\begin{aligned}P(y^*|\bf{y})&=\int_\Lambda P(y^*|\bf{y})P(\lambda|y)dy\\&=\int_0^\infty \frac{e^{-y^*/\lambda}}{\lambda}\frac{(s+\beta)^{\alpha+n}}{\varGamma(\alpha+n)}\lambda^{-(\alpha+n+1)}e^{-(s+\beta)/\lambda}dy\\&=\frac{(s+\beta)^{\alpha+n}}{\varGamma(\alpha +n)}\int_0^\infty e^{-(s+\beta+y^*)/\lambda}\lambda^{-(\alpha+n+1+1)}dy\\&=\frac{(s+\beta)^{\alpha+n}}{\varGamma(\alpha +n)}\frac{\varGamma(\alpha+n+1)}{\varGamma(\alpha+n+1)}\frac{(s+\beta+y^*)^{\alpha+n+1}}{(s+\beta+y^*)^{\alpha+n+1}}\int_0^\infty e^{-(s+\beta+y^*)/\lambda}\lambda^{-(\alpha+n+1+1)}dy\\&=\frac{(s+\beta)^{\alpha+n}}{\varGamma(\alpha +n)}\frac{\varGamma(\alpha+n+1)}{(s+\beta+y^*)^{\alpha+n+1}}\\&=\frac{(s+\beta)^{\alpha+n}}{(s+\beta+y^*)^{\alpha+n+1}}(\alpha+n)\end{aligned}\] Ahora como de \(P(y^*|\bf{y})\) conocemos todos los valores menos \(y^*\) nos queda \[\frac{(s+\beta)^{\alpha+n}}{(s+\beta+y^*)^{\alpha+n+1}}(\alpha+n)=\frac{103237^{22,5}*22,5}{(103237+y^*)^{23.5}}\]

es una función de distribución, por lo que \[Pr(y^*<4000|y)=\int_0^{4000}\frac{103237^{22,5}*22,5}{(103237+y^*)^{23.5}}dy^*=0.5748\]

Octava pregunta:

Pruebe el sistema de hipótesis \(H_0: \lambda=\lambda_0\) frente a \(H_1=\lambda\neq\lambda_0\), con \(\lambda_0=4000\). Para ellos tenga en cuenta que

\[p(\textbf{y}|H_0)=\int_{0}^{\infty}\lambda_0^{-n}\exp\left(-\frac{1}{\lambda_0}\sum_{i=1}^{n}y_i\right)\delta_{\lambda_0}(\lambda)d\lambda\] y

\[p(\textbf{y}|H_1)=\int_{0}^{\infty}\lambda_0^{-n}\exp\left(-\frac{1}{\lambda_0}\sum_{i=1}^{n}y_i\right)\frac{\beta_0^{\alpha_{0}}}{\Gamma(\alpha_0)}\lambda^{-(\alpha_0 +1)}\exp\left(-\frac{\beta_0}{\lambda}\right) d\lambda\]

donde \(\delta_{\lambda_0}(\lambda)\) es la función delta de Dirac. Reporte el factor de Bayes \(B_{10}\) e interprete los resultados.

Desarrollo:

Para nuestro investigador es de interés conocer si hay evidencia estadística de que el parámetro \(\lambda\) tome el valor \(\lambda_0 = 4000\). Para dar respuesta a su interrogante, planteamos el sistema de hipótesis:

\[\left\{ \begin{array}{ll} H_{0}: & \lambda=4000\\ H_{1}: & \lambda\neq 4000 \end{array} \right.\]

Desde el paradigma bayesiano es importante notar que para juzgar este sistema de hipótesis basta con calcular las probabilidades \(P(H_0|\bf{y})\) y \(P(H_1|\bf{y})\), y evaluar cuál de estas dos tiene la probabilidad más alta de ocurrencia, en caso de que \(P(H_0|\bf{y})<P(H_1|\bf{y})\) diremos que hay evidencia estadística para rechazar \(H_0\).

Teniendo en cuenta que no tenemos información o conocimiento a favor de algunas de nuestras hipótesis, asumiremos que \[P(H_0)=P(H_1)=0.5\]

Por definición, \(Pr(H_k|y)\) para \(k=1,2\) es:

\[Pr(H_k|y)=\frac{Pr(y|H_k)*Pr(H_k)}{Pr(y|H_0)*Pr(H_0)+Pr(y|H_1)*Pr(H_1)}\] Del enunciado, tenemos que:

\[\begin{aligned} Pr(y|H_0) &=\int_{0}^{\infty}\lambda_0^{-n}\exp\left(-\frac{1}{\lambda_0}\sum_{i=1}^{n}y_i\right)\delta_{\lambda_0}(\lambda)d\lambda\\ &=\lambda_0^{-n}\exp\left(-\frac{1}{\lambda_0}\sum_{i=1}^{n}y_i\right)\int_{0}^{\infty}\delta_{\lambda_0}(\lambda)d\lambda\\ &\overset{1}=\lambda_0^{-n}\exp\left(-\frac{1}{\lambda_0}\sum_{i=1}^{n}y_i\right) \hspace{0.3cm} (2) \end{aligned}\] Teniendo en cuenta que: 1. \(\int_{0}^{\infty}\delta_{\lambda_0}(\lambda)=1\) (P.A.M.Dirac, 1934)

De manera análoga, tenemos que:

\[\begin{aligned}p(\textbf{y}|H_1)&=\int_{0}^{\infty}\lambda_0^{-n}\exp\left(-\frac{1}{\lambda_0}\sum_{i=1}^{n}y_i\right)\frac{\beta_0^{\alpha_{0}}}{\Gamma(\alpha_0)}\lambda^{-(\alpha_0 +1)}\exp\left(-\frac{\beta_0}{\lambda}\right) d\lambda \\ &=\frac{\beta_0^{\alpha_0}}{\Gamma(\alpha_0)}\int_{0}^{\infty}\lambda^{-(\alpha_0+n+1)}\exp\left(-\frac{1}{\lambda}(s+\beta_0)\right)d\lambda \\ &= \frac{\beta_0^{\alpha_0}}{\Gamma(\alpha_0)}*\frac{\Gamma(\alpha_0+n)}{(\beta_0+s)^{\alpha_0+n}} \hspace{0.5cm} (3) \end{aligned}\] Usando las ecuaciones (2) y (3) tenemos que:

v1<-c(495,541,1461,1555,1603,2201,2750,3468,3516,4319,6622,7728,13159,21194) #Vector de observaciones
s1=sum(v1) #Estadístico suficiente
n1=14 

lo=4000 #lambda_0
ao=8.25 #alpha_0
bo=32625 #beta_0

pyho=((lo)^(-n1))*exp(-s1/lo) #P(y|Ho)
pyh1=((bo^ao)/gamma(ao))*(gamma(ao+n1)/(s1+bo)^(ao+n1)) #P(y|H1)

#Usando la definición de P(H_k|y)

phoy=(pyho*0.5)/(pyho*0.5+pyh1*0.5) # P(H_0|y)
ph1y=(pyh1*0.5)/(pyho*0.5+pyh1*0.5) # P(H_1|y)

phoy

## [1] 0.5610441

ph1y

## [1] 0.4389559

#Sea el factor de Bayes B_10
pyh1/pyho

## [1] 0.782391

Luego, no habría evidencia estádistica para rechazar la hipótesis nula. Y el factor de Bayes \(B_{10}\) nos indica que es insignificativa la diferencia.

Novena pregunta:

Experimentación adicional bajo las mismas condiciones con otro tipo de alambre (Tipo 2) produce los siguientes resultados

\[294\hspace{0.5cm} 569 \hspace{0.5cm} 766 \hspace{0.5cm} 1576\hspace{0.5cm} 1602\hspace{0.5cm} 2015\hspace{0.5cm} 2166\hspace{0.5cm} 3885\hspace{0.5cm} 8141\hspace{0.5cm} 10285\] Considerando modelos independientes de la forma \(y_{i,k}|\lambda_k\stackrel{\text { iid }}{\sim}\text{Exp}(\lambda_k)\) con \(\lambda_k\sim \text{GI}(\alpha_0,\beta_0)\), para \(i=1,...,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.

Desarrollo:

El sistema de hipótesis que queremos contrastar es: \[\left\{ \begin{array}{ll} H_{0}: & \lambda_1=\lambda_2=\lambda\\ H_{1}: & \lambda_1\neq \lambda_2 \end{array} \right.\]

Por definición, tenemos que: \[\begin{aligned}Pr(y_1,y_2|H_0)&=\int_0^{\infty} p(y_1,y_2|\lambda)p(\lambda)d\lambda\\ &\overset{1}= \int_0^{\infty} p(y_1|\lambda)p(y_2|\lambda)p(\lambda)d\lambda \\ &= \int_0^{\infty} \frac{1}{\lambda^{n_1}}\exp\left\lbrace-\frac{s_1}{\lambda}\right\rbrace*\frac{1}{\lambda^{n_2}}\exp\left\lbrace-\frac{s_2}{\lambda}\right\rbrace*\frac{\beta_0^{\alpha_0}}{\Gamma(\alpha_0)}\lambda^{-(\alpha_0+1)}\exp\left\lbrace-\frac{\beta_0}{\lambda}\right\rbrace d\lambda\\ &= \frac{\beta_0^{\alpha_0}}{\Gamma(\alpha_0)}\int_0^{\infty}\lambda^{-(\alpha_0+n_1+n_2+1)}\exp\left\lbrace-\frac{\beta_0+s1+s2}{\lambda}\right\rbrace d\lambda\\ &\overset{2}= \frac{\beta_0^{\alpha_0}}{\Gamma(\alpha_0)}* \frac{\Gamma(\alpha_0+n_1+n_2)}{(\beta_0+s1+s2)^{\alpha_0+n1+n2}} \hspace{0.5cm} (4) \end{aligned} \] Teniendo presente que: 1. Se justifica en el supuesto de independencia condicional dado \(\lambda\) de las variables \(y_i,k\), para \(i=1,...,n\) y \(k=1,2\). 2. La expresión de la integral es el kernel(núcleo) de una distribución Gamma Inversa con parámetros \(\alpha_0+n_1+n_2\) y \(\beta_0+s_1+s_2\)

Por otro lado, tenemos que:

\[\begin{aligned} \Pr(y_1,y_2|H_1)&= \int_0^{\infty}\int_0^{\infty}p(y_1,y_2|\lambda_1,\lambda_2)p(\lambda_1)p(\lambda_2)d\lambda_1 d\lambda_2 \\&\overset{1}=\int_0^{\infty}\int_{0}^{\infty}p(y_1|\lambda_1)p(y_2|\lambda_2)p(\lambda_1)p(\lambda_2)d\lambda_1d\lambda_2 \\&= \int_0^{\infty}p(y_1|\lambda_1)p(\lambda_1)d\lambda_1\int_0^{\infty}p(y_2|\lambda_2)p(\lambda_2)d\lambda_2\\ &\overset{2}= \frac{\Gamma(\alpha_0+n1)}{(\beta_0+s1)^{\alpha_0+n1}}*\frac{\Gamma(\alpha_0+n2)}{(\beta_0+s2)^{\alpha_0+n2}} \hspace{0.5cm} (5) \end{aligned}\] De (4), (5), la definición de \(Pr(H_k|y)\) para \(k=0,1\) y usando los datos que requiere el investigador encontramos que:

v2<-c(294,569,766,1576,1602,2015,2166,3885,8141,10285) #Vector de valores
s2=sum(v2) #Estadístico suficiente del alambre tipo 2
n2=10 #Tamaño de la muestra

py1y2ho=((bo^ao)/(gamma(ao)))*((gamma(ao+n1+n2))/((bo+s1+s2)^(ao+n1+n2))) #P(y1,y2|H_0)
py1y2h1=((gamma(ao+n1))/((bo+s1)^(ao+n1)))*((gamma(ao+n2))/((bo+s2)^(ao+n2))) #P(y1,y2|H_1)

#Calculamos las probabilidades de interés
ph0y1y2=py1y2ho*0.5/(py1y2ho*0.5+py1y2h1*0.5) #P(H_0|y_1,y_2) 
ph1y1y2=py1y2h1*0.5/(py1y2ho*0.5+py1y2h1*0.5) #P(H_1|y_1,y_2)

ph0y1y2

## [1] 1

ph1y1y2

## [1] 2.784724e-67

py1y2h1/py1y2ho #Factor de Bayes B_10

## [1] 2.784724e-67

py1y2ho/py1y2h1 #Factor de Bayes B_01

## [1] 3.59102e+66

Dado que \(P(H_0|y_1,y_2)>P(H_1|y_1,y_2)\) no hay evidencia para rechazar la diferencia entre los parámetros de las distribuciones de las cuáles surgen los datos de \(y_1\) y \(y_2\). Asimismo, el factor de Bayes \(B_{10}\) también nos indica que no hay diferencias significativas. Sin embargo, el valor del factor de Bayes \(B_{01}\) nos indica que hay evidencia decisiva a favor de \(H_0\)

Décima pregunta:

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.

Desarrollo:

Distribución Posterior de alambre tipo 2

Datos tomados del segundo experimento

y2=c(294,569,766,1576,1602,2015,2166,3885,8141,10285)

Cantidad de datos tomados \(n_2=10\)

n2=length(y2);n2

## [1] 10

Estadístico \(s_2=\sum_{k=1}^ny_k= 31299\)

s2=sum(y2);s2

## [1] 31299

La distribución posterior es \(\lambda_2|y_2\sim GI(\alpha+n_2,\beta+s_2)\) con los mismos parámetros encontrados en el punto 4.

set.seed(1)
ap2=a+n2
bp2=b+s2
sim1=rinvgamma(n=1000000,ap,bp)
sim2=rinvgamma(n=1000000,ap2,bp2)

Distribución predictiva posterior

Encontramos la predictiva posterior dados por la distribución \(y_{i,k}^*|\lambda_k\overset{iid}{\sim}Exp(\lambda_k)\) donde \(\lambda_{k}\) son simulados por el método de montecarlo por los valores de las distribuciones posteriores dadas por \(\lambda_k|y_{i,k}\sim GI(\alpha+n_k,\beta+s_k)\).

B=100000 #Número de simulaciones
set.seed(1)
y1p=rpois(n=B,lambda = sim1);
y2p=rpois(n=B,lambda = sim2);

Idoneidad del modelo alambre tipo 1

La media de los datos observada para el primer experimento es

ob1=round(s/n,3);ob1

## [1] 5043.714

Valor p predictivo posterior del primer experimento es

mean(y1p<ob1)

## [1] 0.62497

Por lo tanto el modelo es adecuado para caracterizar la media de la población.

Idoneidad del modelo alambre tipo 2

La media de los datos observada para el segundo experimento es

ob2=round(s2/n2,3);ob2

## [1] 3129.9

Valor p predictivo posterior del segundo experimento es

mean(y2p<ob2)

## [1] 0.28483

Por lo tanto el modelo es adecuado para caracterizar la media de la población

Tiempo de falla de un alambre

Por: Sebastián Franco y Rodolfo Alejandro García

Caso de estudio:

Preguntas del caso de estudio:

Primera pregunta:

Desarrollo:

Segunda pregunta:

Desarrollo:

Tercera pregunta:

Desarrollo:

Cuarta pregunta:

Desarrollo:

Quinta pregunta:

Desarrollo:

Sexta pregunta:

Desarrollo:

Método Bayesiano

Método Frecuentista

Método de Bootstrap

Séptima pregunta:

Desarrollo:

Octava pregunta:

Desarrollo:

Novena pregunta:

Desarrollo:

Décima pregunta:

Desarrollo:

Distribución Posterior de alambre tipo 2

Distribución predictiva posterior

Idoneidad del modelo alambre tipo 1

Idoneidad del modelo alambre tipo 2

Referencias