Here are a couple of algorithms that will do robust estimation of the parameter of the geometric distribution, along with an example. I'll write the geometric distribution as $\text{p}(x | p) = p(1-p)^x$.
Method 1 is based on the fact that $\text{p}(x+1|p)/\text{p}(x|p) = 1-p$. We can extend this to:
$\frac{\sum_{i=0}^a \text{p}(i|p)}{\sum_{i=a+1}^{2a}\text{p}(i|p)} = (1-p)^a$
If we specify $a$ low enough so that the observed frequencies are not contaminated by the outliers, we can use this to derive a robust estimator of $p$:
func.rob <- function(x, a) {
1 - (sum(x>=a & x<(2*a)) / sum(x<a))^(1/a)
}
Method 2 is a trimmed mean. We specify a trimming percentage $\alpha$ and the algorithm calculates the mean of all the values less than or equal to the $1-\alpha$ quantile of the data. We then find $p$ such that the geometric distribution has the same mean over the same range.
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
}
We use minimization of the square rather than rootfinding because it's not always easy to bracket the root; both low $p$ and high $p$ will give very low trimmed means for cutoffs that are small relative to the true mean of the distribution. (Generally this isn't a good thing to do, however.) Additionally, it seems to me that it is possible for there to be no solution to this problem, especially for smaller samples; I haven't observed it in tens of thousands of runs with samples of size 250, though.
Here are some comparisons against the MLE for uncontaminated and contaminated data, sample size 250, $p=0.4$. First, 10,000 runs with uncontaminated data:
> 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
>
As we can see, the robust estimator appears unbiased, while the trimmed estimator is biased a little low; however, the robust estimator is quite a bit less efficient at the true model than the MLE, while the trimmed estimator does quite well.
Now we change 4% of the values to 25, a clear outlier for this choice of parameters:
> # 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
>
The MLE has failed, while the ratio-based robust estimator is doing well and the trimmed estimator has fallen a little lower; this is because the use of a quantile instead of a raw number results in a biased estimate of the quantile, thanks to all the outliers being on the high side. Still, it's much better than the MLE, and it's easy enough to replace the quantile with a raw number that you are confident lies below the "outlier" level.
Edit: (additional answer stuff)
If you decide to use algorithms such as the above, I recommend doing simulation experiments such as I've done to help you understand their properties.
As for identifying outliers - that's a chancy thing, where judgement plays a role. If you suspect the frequency of outliers is low, say, $\le 1\%$, you might use the robust estimate of $p$, calculate the $99^{\text{th}}$ percentile of the distribution, and put all the data points that fall above that in the "suspect" category. If you have some way of going back to the source of the data, which it appears from a comment above that you do, that would help clarify whether the points were indeed outliers. Alternatively, you could throw out everything above the, for example, $99.9^{\text{th}}$ percentile, given that your sample is well under 1000 you're not likely to be throwing out more than a very few good data points. I prefer the "filter, then examine" approach (when it's feasible) to the "delete" approach, myself.
