library(dplyr)
set.seed(123)lab04b: Contraste de hipótesis
XS3310 Teoría Estadística
Este documento ilustra de manera intuitiva, por medio de simulaciones, el concepto del contraste de hipótesis.
1 Paquetes
2 Ejemplo con la distribución exponencial
En la clase vimos el siguiente ejemplo:
Sea \(X_{1}, X_{2}, ... , X_{n}\) una muestra aleatoria tal que \(X_{j} \sim Exp(\theta)\) y suponga que se desean contrastar las hipótesis \(H_{0}: \theta = \theta_0\) contra \(H_{1}: \theta \neq \theta_0\). Encuentre el contraste de razón de verosimilitudes para un tamaño de \(\alpha_0\).
Suponga que la \(H_0: \theta_0=1\) es cierta. Generamos \(K=1000\) muestras independientes de tamaño \(10\).
K <- 1000
n <- 10
X <- list()
for(i in 1:K){
X[[i]] <- rexp(n,1)
}2.1 Estadístico de razón de verosimilitud exacto
Calculamos el estadístico de razón de verosimilitud para cada repetición \(k=1,...,K\) y luego verificamos cuántas de ellas rechaza la hipótesis nula de acuerdo con la regla de decisión \[V=\frac{2n\bar{X}}{\theta_0}< \chi^{2}_{\frac{\alpha_0}{2},2n},~~ \text{o}~~ V=\frac{2n\bar{X}}{\theta_0} > \chi^{2}_{1-\frac{\alpha_0}{2},2n}.\]
Aquí se construye una función para el cálculo del estadístico de contraste para una muestra dada X:
RV <- function(X,H0=1){
n <- length(X)
RV <- 2*n*mean(X)/H0
return(RV)
}Luego, aplicamos la misma función a cada elemento de la lista X, y verificamos que 44 de \(1000\) se rechaza la hipótesis, o sea próximo al tramaño del contraste \(5\%\).
RV_resultado <- sapply(X,RV)
sum(!between(RV_resultado,qchisq(p=0.025,2*n),qchisq(p=0.975,2*n)))/K[1] 0.0442.2 Estadístico de razón de verosimilitud (asintótico)
Por otro lado, usando el teorema de Wilks, calculamos el estadístico \[G = -2\ln(\lambda) = -2n\left[ \ln(\bar{X}) - \ln(\theta_0) + 1 - \frac{\bar{X}}{\theta_0} \right]\]
RV_wilks <- function(X,H0=1){
n <- length(X)
RV_wilks <- -2*n*(log(mean(X))-log(H0)+1-(mean(X)/H0))
return(RV_wilks)
}
RV_wilks_resultado <- sapply(X,RV_wilks)sim2 <- sum(RV_wilks_resultado>qchisq(p=0.95,1))/KEncontramos que un 0.045 \(\%\) de las \(K=1000\) muestras rechaza la hipótesis, o sea próximo al tramaño del contraste \(5\%\).
sum(RV_wilks_resultado>qchisq(p=0.95,1))/K[1] 0.045Ambos resultados son muy próximos y nos damos cuenta que con \(n=10\), el resultado asintótico aproxima bastante bien.
2.3 Estadístico de Wald
Finalmente, utilizando el estadístico de Wald: \[W=(\hat{\theta}_n - \theta_0)I(\hat{\theta})(\hat{\theta}_n - \theta_0)\]
estad_wald <- function(X,H0=1){
n <- length(X)
estad_wald <- (mean(X)-H0)^2*(n/mean(X)^2)
return(estad_wald)
}
estad_wald_resultado <- sapply(X,estad_wald)sim3 <- sum(estad_wald_resultado>qchisq(p=0.95,1))/KEncontramos que un 0.045 \(\%\) de las \(K=1000\) muestras rechaza la hipótesis, o sea próximo al tramaño del contraste \(5\%\).
sum(estad_wald_resultado>qchisq(p=0.95,1))/K[1] 0.098Los 3 resultados son muy próximos y nos damos cuenta que con \(n=10\), el resultado asintótico aproxima bastante bien.
2.4 Con \(n=50\) y \(\alpha=0.05\)
K <- 1000
n <- 50
X <- list()
for(i in 1:K){
X[[i]] <- rexp(n,1)
}
RV_resultado <- sapply(X,RV)
sum(!between(RV_resultado,qchisq(p=0.025,2*n),qchisq(p=0.975,2*n)))/K[1] 0.046RV_wilks_resultado <- sapply(X,RV_wilks)
sum(RV_wilks_resultado>qchisq(p=0.95,1))/K[1] 0.046estad_wald_resultado <- sapply(X,estad_wald)
sum(estad_wald_resultado>qchisq(p=0.95,1))/K[1] 0.0582.5 Con \(n=50\) y \(\alpha=0.1\)
K <- 1000
n <- 50
X <- list()
for(i in 1:K){
X[[i]] <- rexp(n,1)
}
RV_resultado <- sapply(X,RV)
sum(!between(RV_resultado,qchisq(p=0.05,2*n),qchisq(p=0.95,2*n)))/K[1] 0.112RV_wilks_resultado <- sapply(X,RV_wilks)
sum(RV_wilks_resultado>qchisq(p=0.9,1))/K[1] 0.114estad_wald_resultado <- sapply(X,estad_wald)
sum(estad_wald_resultado>qchisq(p=0.9,1))/K[1] 0.121