Robust estimation of a geometric random variable

I have a bunch of data which is assumed to be instances of a geometric random variable with outliers. How can I do a robust estimation of the parameter $p$ so that the effect of outliers is minimized?

As a part of the estimation process, I also need to know which are the outliers in the data.

I am looking for some solution using R, matlab or C/C++.

EDIT 1

Is any algorithm available which can be used for robust estimation of geometric random variable, though not readily available as a function in R or matlab ?

EDIT 2

The maximum likelihood estimate of p of geometric random variable is the mean of the instance values. So if we do robust estimate of mean of the instance values, can we say that we are doing robust estimation of the underlying geometric random variable? If so, which method is most suitable for doing the robust mean estimation. (I am a newbie in statistics and R)

• How are you identifying outliers? Typically an outlier is a point that has a different distribution from the rest of the data set (perhaps due to an error in data collection/recording). In many cases, it is problematic to identify outliers purely based on their numeric value (except in cases like where you see an observed human height of 30 feet). How do you know large values aren't part of the true data generating process? – Macro May 1 '12 at 14:56
• Well, @macro, one excellent reason to find a robust estimator $\hat{p}$ is precisely to respond to your question: the robust estimator will be relatively uninfluenced by outliers and therefore is expected to be a good summary of the mass of data. As such, the geometric distribution with parameter $\hat{p}$ becomes a reference against which the data can be compared in order to address questions about their distribution, including checking the expectation that they are geometric with possibly a small proportion of extreme (outlying) values. – whuber May 1 '12 at 15:59
• @Macro my data comes from face/fingerprint recognition systems and by looking at the data having excessively large values, I am able to see that they typically correspond to poor registration/noisy image etc. And I completely agree with the reasoning given by whuber. But the question is, how to do robust estimation? – suresh May 1 '12 at 19:36
• How large a sample do you have? Do you have a substantial fraction of zeroes, or is the distribution really spread out? – jbowman May 2 '12 at 1:30
• @jbowman the data has substantial fraction of zeros. The datasize may be a few hundreds. Btw, why did you ask this question? – suresh May 2 '12 at 21:02

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)  0.4017998 > mean(phat.trim)  0.3651195 > mean(phat.mle)  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.

• Thank you very much for the answer. But how do I find out which values are the outliers? – suresh May 3 '12 at 2:32