Consideramos el test de Kolmogorov-Smirnov con hipótesis nula compuesta: \[ H_0: F \in \{F_{\boldsymbol\theta}:\boldsymbol\theta\in\Theta\subset{\mathbb R}^r\}, \] es decir, queremos contrastar si \(F\) pertenece a una familia paramétrica de distribuciones de probabilidad.

Por ejemplo, queremos contrastar si los datos del fichero kevlar.txt siguen una distribución exponencial de tasa \(\lambda>0\), cuya función de distribución es \[ F_{\lambda}(x) = 1 - e^{-\lambda x}, \quad x>0. \] Por tanto, el contraste que queremos hacer es: \[ H_0: F \in \{F_{\lambda}:\lambda\in(0,\infty)\subset{\mathbb R}\}. \] Paso 1: Calculamos la estimación de máxima verosimilitud del parámetro desconocido.

kev <- scan("http://verso.mat.uam.es/~amparo.baillo/MatEstII/kevlar.txt")
m <- mean(kev)
emvlambda <- 1/m
emvlambda

## [1] 0.9765455

Paso 2: Calculamos el estadístico de Kolmogorov-Smirnov con el parámetro estimado, \(\hat{D}_{n} = \|F_{n}-F_{\hat{\lambda}}\|_{\infty}\):

KStest <- ks.test(kev,pexp,emvlambda)
Dn <- KStest$statistic
Dn

##          D 
## 0.08969841

Paso 3: La región de rechazo del contraste al nivel de significación \(\alpha\) es \(R=\{\hat{D}_{n}>c\}\), donde el nivel crítico \(c\) es el percentil 95 de la distribución de \(\hat D_n\) bajo \(H_0\) y se tiene que aproximar mediante métodos Monte Carlo (simulación). La idea es que, si \(H_0\) fuera cierta, entonces \(F\) sería una distribución exponencial con valor aproximado \(\hat{\lambda}\). Vamos a generar muchas (\(N\)) muestras Monte Carlo (del mismo tamaño \(n\) que la original) de la distribución exponencial de parámetro \(\hat\lambda\). Para cada muestra Monte Carlo calculamos la e.m.v. de \(\lambda\), que la denotaremos \(\hat{\hat\lambda}\), y luego el valor del estadístico de Kolmogorov-Smirnov: \[ \begin{array}{cccccc} \mbox{Muestra Monte Carlo 1} & x_1^{(1)},\ldots,x_n^{(1)} & \rightarrow & \hat{\hat\lambda}_1 & \rightarrow & \hat{\hat D}_n^{(1)} = \|F_n^{(1)}-F_{\hat{\hat\lambda}_1}\|_\infty \\ \mbox{Muestra Monte Carlo 2} & x_1^{(2)},\ldots,x_n^{(2)} & \rightarrow & \hat{\hat\lambda}_2 & \rightarrow & \hat{\hat D}_n^{(2)} = \|F_n^{(2)}-F_{\hat{\hat\lambda}_2}\|_\infty \\ \vdots \\ \mbox{Muestra Monte Carlo N} & x_1^{(N)},\ldots,x_n^{(N)} & \rightarrow & \hat{\hat\lambda}_N & \rightarrow & \hat{\hat D}_n^{(N)} = \|F_n^{(N)}-F_{\hat{\hat\lambda}_N}\|_\infty \end{array} \] Aproximamos el valor de \(c\) mediante el percentil 95 de \(\hat{\hat D}_n^{(1)}, \hat{\hat D}_n^{(2)},\ldots,\hat{\hat D}_n^{(N)}\).

n <- length(kev)
N <- 1000
DnMC <- rep(0,N)
for (i in 1:N){
  MuestraMC <- rexp(n,rate=emvlambda)
  emvlambdaMC <- 1/mean(MuestraMC)
  KStestMC <- ks.test(MuestraMC,pexp,emvlambdaMC)
  DnMC[i] <- KStestMC$statistic
}
c <- quantile(DnMC,probs=0.95)
c

##       95% 
## 0.1074688