Aquí hay un par de algoritmos que va a hacer sólida la estimación del parámetro de la distribución geométrica, junto con un ejemplo. Voy a escribir la distribución geométrica como $\text{p}(x | p) = p(1-p)^x$.
Método 1 se basa en el hecho de que $\text{p}(x+1|p)/\text{p}(x|p) = 1-p$. Podemos extender este:
$\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 baja para que las frecuencias observadas no son contaminados por los valores extremos, se puede utilizar para obtener 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 limitada. Podemos especificar un recorte porcentaje $\alpha$ y el algoritmo calcula la media de todos los valores inferiores o iguales a la $1-\alpha$ los cuantiles de los datos. Podemos encontrar $p$ de manera tal que la distribución geométrica tiene la misma media que a través de la misma gama.
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 de la plaza en lugar de rootfinding porque no siempre es fácil de soporte de la raíz; ambos bajo $p$ y alto $p$ le dará muy baja tapizados significa para los puntos de corte que son pequeños en relación a la verdadera media de la distribución. (En general, esto no es una buena cosa que hacer, sin embargo.) Además, a mí me parece que es posible que no haya ninguna solución a este problema, especialmente para las pequeñas muestras; no he observado en decenas de miles de carreras con muestras de tamaño 250, aunque.
Aquí están algunas comparaciones contra el MLE para no contaminada y contaminada de datos, el tamaño de la muestra 250, $p=0.4$. En primer lugar, de 10.000 pistas con incontaminada de datos:
> 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 aparece imparcial, mientras que el tapizado estimador es sesgado un poco baja; sin embargo, el estimador robusto es un poco menos eficiente en el verdadero modelo de la MLE, mientras que el tapizado estimador hace bastante bien.
Ahora cambiamos el 4% de los valores de 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 la relación basada en sólidas estimador está haciendo bien y los tapizados estimador ha caído un poco menor; esto es debido a que el uso de un cuantil en lugar de un número sin procesar los resultados de una estimación sesgada de la cuantil, gracias a todos los valores atípicos de estar en el lado de alta. Aún así, es mucho mejor que el MLE, y es bastante fácil de reemplazar el cuantil con una cruda número que usted está seguro se encuentra por debajo del "valor atípico".
Edit: (adicional contesto)
Si usted decide utilizar algoritmos como el de arriba, yo recomiendo hacer experimentos de simulación como yo lo he hecho para ayudarle a entender sus propiedades.
Como para la identificación de valores atípicos - que es una cosa arriesgada, donde el juicio juega un papel importante. Si usted sospecha que la frecuencia de los valores extremos es baja, digamos, $\le 1\%$, se puede usar la estimación robusta de $p$, calcular el $99^{\text{th}}$ percentil de la distribución, y poner todos los puntos de datos que caen de arriba que en la "sospecha" de la categoría. Si usted tiene alguna manera de volver a la fuente de los datos, de la cual se desprende de un comentario anterior que hacer, que ayudaría a aclarar si los puntos de hecho fueron los valores atípicos. Alternativamente, usted podría tirar todo por encima de el, por ejemplo, $99.9^{\text{th}}$ percentil, dado que la muestra está bien por debajo de los 1000 que no te va a tirar más de un muy buen par de puntos de datos. Yo prefiero el "filtro", luego de examinar el enfoque de" (cuando es posible) a la "eliminar" en el enfoque, a mí mismo.
