Aquí hay un par de algoritmos que harán una estimación robusta del parámetro de la distribución geométrica, junto con un ejemplo. Escribiré la distribución geométrica como $\text{p}(x | p) = p(1-p)^x$ .
El método 1 se basa en que $\text{p}(x+1|p)/\text{p}(x|p) = 1-p$ . Podemos ampliar esto a:
$\frac{\sum_{i=0}^a \text{p}(i|p)}{\sum_{i=a+1}^{2a}\text{p}(i|p)} = (1-p)^a$
Si especificamos $a$ lo suficientemente bajo como para que las frecuencias observadas no estén contaminadas por los valores atípicos, podemos utilizarlo para derivar un estimador robusto de $p$ :
func.rob <- function(x, a) {
1 - (sum(x>=a & x<(2*a)) / sum(x<a))^(1/a)
}
El método 2 es una media recortada. Se especifica un porcentaje de recorte $\alpha$ y el algoritmo calcula la media de todos los valores menores o iguales al $1-\alpha$ cuantil de los datos. A continuación, encontramos $p$ tal que la distribución geométrica tiene la misma media en el mismo rango.
func.trim <- function(x, a) {
require(MASS)
geom.tm <- function(p, xbar, cutoff) {(sum(dgeom(1:cutoff, p)*c(1:cutoff))-xbar)^2}
cutoff <- quantile(x, a)
xbar <- mean(x[x <= cutoff])
optim(0.5, geom.tm, lower=1.0e-05, upper=1-1.0e-05, method="L-BFGS-B",
xbar=xbar, c=cutoff)$par
}
Utilizamos la minimización del cuadrado en lugar de la búsqueda de la raíz porque no siempre es fácil poner entre paréntesis la raíz; tanto la baja $p$ y alta $p$ dará medias recortadas muy bajas para cortes que son pequeños en relación con la media real de la distribución. Además, me parece que es posible que no haya solución a este problema, especialmente para muestras pequeñas; sin embargo, no lo he observado en decenas de miles de ejecuciones con muestras de tamaño 250.
A continuación se presentan algunas comparaciones con la MLE para datos no contaminados y contaminados, con un tamaño de muestra de 250, $p=0.4$ . Primero, 10.000 ejecuciones con datos no contaminados:
> func.mle <- function(x) 1/(1+mean(x))
>
> # Uncontaminated data
> z <- matrix(rgeom(2500000, 0.4), 10000, 250)
> phat.rob <- apply(z, 1, func.rob, a=4)
> phat.trim <- apply(z, 1, func.trim, a=0.75)
> phat.mle <- apply(z, 1, func.mle)
>
> cat(" Robust estimator sample mean, std. dev.: ",mean(phat.rob)," ",sd(phat.rob),"\n",
+ "Trim estimator sample mean, std. dev.: ",mean(phat.trim)," ",sd(phat.trim),"\n",
+ "ML estimator sample mean, std. dev.: ",mean(phat.mle)," ",sd(phat.mle),"\n",
+ "Relative efficiency (robust, trim): ",var(phat.mle)/var(phat.rob)," ",
+ var(phat.mle)/var(phat.trim), "\n")
Robust estimator sample mean, std. dev.: 0.4017054 0.03044539
Trim estimator sample mean, std. dev.: 0.3874961 0.01849624
ML estimator sample mean, std. dev.: 0.4007747 0.01972546
Relative efficiency (robust, trim): 0.4197697 1.137332
>
Como podemos ver, el estimador robusto parece insesgado, mientras que el estimador recortado tiene un sesgo un poco bajo; sin embargo, el estimador robusto es bastante menos eficiente en el modelo verdadero que el MLE, mientras que el estimador recortado lo hace bastante bien.
Ahora cambiamos el 4% de los valores a 25, un claro valor atípico para esta elección de parámetros:
> # Contaminated data 4% @ 25
> z[,1:10] <- 25
> phat.rob <- apply(z, 1, func.rob, a=4)
> phat.mle <- apply(z, 1, func.mle)
> phat.trim <- apply(z, 1, func.trim, a=0.75)
> mean(phat.rob)
[1] 0.4017998
> mean(phat.trim)
[1] 0.3651195
> mean(phat.mle)
[1] 0.290966
>
El MLE ha fallado, mientras que el estimador robusto basado en la proporción lo hace bien y el estimador recortado ha caído un poco más bajo; esto se debe a que el uso de un cuantil en lugar de un número crudo da lugar a una estimación sesgada del cuantil, gracias a que todos los valores atípicos están en el lado alto. Aun así, es mucho mejor que el MLE, y es bastante fácil sustituir el cuantil por un número en bruto que esté seguro de que se encuentra por debajo del nivel "atípico".
Editar: (cosas adicionales de la respuesta)
Si decides utilizar algoritmos como los anteriores, te recomiendo que hagas experimentos de simulación como los que yo he hecho para ayudarte a entender sus propiedades.
En cuanto a la identificación de los valores atípicos, se trata de una cuestión delicada en la que interviene el juicio. Si sospechas que la frecuencia de los valores atípicos es baja, por ejemplo, $\le 1\%$ se podría utilizar la estimación robusta de $p$ calcula el $99^{\text{th}}$ percentil de la distribución, y poner todos los puntos de datos que caen por encima de eso en la categoría "sospechoso". Si tienes alguna forma de volver a la fuente de los datos, lo que parece, por un comentario anterior, que tienes, eso ayudaría a clarificar si los puntos eran realmente atípicos. También podrías descartar todo lo que esté por encima de los, por ejemplo, $99.9^{\text{th}}$ Dado que la muestra es muy inferior a 1.000, no es probable que se descarten más que unos pocos puntos de datos buenos. Yo mismo prefiero el enfoque de "filtrar y luego examinar" (cuando es factible) al de "eliminar".
