Pruebas de Bondad de Ajuste en R

Las pruebas de bondad de ajuste consisten en comparar nuestra función de distribución empírica $F_n$ generada por nuestros datos, y una función de distribución cualesquiera $F_0$, y contrastaremos hipótesis del estilo: $$H_0:F=F_0,\quad H_1:F\neq F_0$$ Donde $F$ es la distribución real de nuestros datos.
Consideremos una m.a. $X_1,X_2,...,X_n$, de la cual definimos nuestra  función de distribución empírica como sigue $$F_n(x)=\frac{1}{n}\sum_{i=1}^{n}\mathbb{1}_{\{X_i\le x\}}$$ Es decir a cada dato se le asigna un probabilidad de que ocurra de $1/n$, algunas propiedades de utilidad son las siguientes:
1. $\mathbb{E}(F_n(x))=F(x)=\mathbb{P}(X_1\le x)$
2. $\sqrt{n} (F_n(x)-F(x))\overset { d }{ \rightarrow  } N(0,F(x)(1-F(x)))$
3. $\underset {  x\in \mathbb{R} }{\sup  }|F_n(x)-F(x)|\overset { \mathbb{P} }{ \rightarrow  } 0$

Consideraremos 3 tipos de pruebas que nos darán criterios para aceptar o rechazar nuestras hipótesis.

1.- Kolgomorov-Smirnov 

Para esta prueba consideraremos el estadístico $$D_n=\sqrt{n}\underset {  x\in \mathbb{R} }{\sup  }|F_n(x)-F_0(x)|$$ El cual esta motivado por las propiedades 1 y 2 de arriba, lo que mide el estadístico es el máximo de la distancia entre la distribución de la hipótesis y la empírica, por lo que nuestra búsqueda sera en $D_n$'s pequeñas, podemos fijar un nivel de confianza, es decir aceptar la hipotesis $H_0$ si $D_n\le c$ y en otro caso rechazar. 
Se puede probar que el estadistico efectivamente no depende de $F$ y que $$\mathbb{P}(D_n\le t)\rightarrow K(t)=1-2\sum_{i=1}^n (-1)^{i-1}e^{-2i^2 t}$$ Es decir converge en distribución a K -la distribución de Kolgomorov-Smirnov- Gracias a esto podemos calcular de que tamaño son nuestros intervalos $\{ D_n\le c \}$, si $n$ es lo suficientemente grande $$\alpha=\mathbb{P}(D_n\ge c | H_0)\approx 1- K(c)$$
Así es como podremos calcular tamaños computacionalmente, adicionalmete se pueden tabular los valores de $K(.)$ si fuera necesario. 
Uso de R

En R haremos uso de la función ks.test() que se alimenta del vector de datos y la función de distribución con la que se quiere comparar. Como en el siguiente ejemplo :
> y<-c(0.3902,0.508,0.9787,0.888,0.1024,
0.6723,0.313605,0.7254,0.9409,0.05277)
> ks.test(y,runif,0,1)

        One-sample Kolmogorov-Smirnov test

data:  y
D = 0.75171, p-value = 2.304e-06
alternative hypothesis: two-sided

En este caso se compara la muestra y con una uniforme $(0,1)$, podemos establecer un p-value y poder contrastar es decir rechazar $H_0$ si $p(X)\le \alpha$ 

2.  Cramer-von Mises 

Podemos darle otro enfoque en la manera en la que comparamos nuestra distancia entre la distribución teórica y la empírica, en esta prueba se comparará todos los puntos de la curva, nuestro estadístico de prueba es el siguiente; $$W_n^2=n\int_{-\infty}^{\infty}[F_n(x)-F(x)]^2\varphi (F(x)) dF(x)$$ donde $\varphi (x)\ge 0$. 
La función $\varphi$ nos ayudara a ponderar ciertas regiones de la función de acumulación donde tenemos mas interés en estimar de mejor manera, el caso de Cramer-von Mises se dará igual importancia a cualquier región por lo que $\varphi(x)=1$, así nuestro estadistico de prueba es: $$\omega_n^2=n\int_{-\infty}^{\infty}[F_n(x)-F(x)]^2 dF(x)$$
La expresión no se ve nada agradable, pero puede ser aproximada de la siguiente manera: $$\omega_n^2=\sum_{i=1}^{n}[F(x)-\frac{2i-1}{2n}]^2+\frac{1}{12n}$$ Obviamente se rechazarán valores grandes de $\omega_n^2$ 
Uso de R
Se usa el paquete goftest y de ahi la función cvm.test, su uso es analogo al de ks.test. 
> y<-c(0.3902,0.508,0.9787,0.888,0.1024,
0.6723,0.313605,0.7254,0.9409,0.05277)
> cvm.test(y,runif,0,1)

        Cramer-von Mises test of goodness-of-fit
        Null hypothesis: distribution ‘runif’

data:  y

omega2 = 0.1048, p-value = 0.5706

Como vemos esta prueba es mas rigurosa que k-s. 

3. Anderson-Darling

En este caso $\varphi (x)=\frac{1}{x(1-x)}$, es de esta manera pues en los valores cercanos a $0$ y $1$ se tiene una mayor ponderación, así se le da mayor importancia a los valores de las colas.
Nuestro estadístico de prueba será el siguiente $$A_n^2=n\int_{-\infty}^{\infty}\frac{[F_n(x)-F(x)]^2}{F(x)(1-F(x))}dF(x)$$ afortunadamente igual podrá se aproximado de la siguiente manera: $$A_n^2=-n-\frac{1}{n}\sum_{i=1}^n (2i-1)[\log u_{(i)} -\log (1-u_{(n-j+1)}]$$ Donde $F(x_{(i)})=u_{(i)}$. Igualmente nos interesan valores pequeños de $A_n^2$ por lo que se rechazarán los grandes. 
Uso de R
Del paquete goftest se usa la función ad.test y se usa de igual forma que las anteriores:
> y<-c(0.3902,0.508,0.9787,0.888,0.1024,
0.6723,0.313605,0.7254,0.9409,0.05277)
> ad.test(y,runif,0,1)

        Anderson-Darling test of goodness-of-fit
        Null hypothesis: distribution ‘runif’

data:  y

An = 0.99798, p-value = 0.3561

Gracias por su visita 

Comentarios