# Robust mean estimation with O(1) update efficiency

I am looking for a robust estimation of the mean that has a specific property. I have a set of elements for which I want to calculate this statistic. Then, I add new elements one at a time, and for each additional element I would like to recalculate the statistic (also known as an online algorithm). I would like this update calculation to be fast, preferably O(1), i.e. not dependent on the size of the list.

The usual mean has this property that it can be updated efficiently, but is not robust to outliers. Typical robust estimators of the mean, like inter-quartile mean and trimmed mean cannot be updated efficiently (since they require maintaining a sorted list).

I would appreciate any suggestions for robust statistics that can be calculated/updated efficiently.

• Why not just use an initial segment of the data--such as the first 100 or first 1000 or whatever--to erect "fences" for screening outliers? You don't have to update them again, so there's no need to maintain additional data structures. – whuber Apr 18 '13 at 16:50
• @whuber I cannot guarantee that the initial sample will represent the rest of the data. For example, the order in which I am given the data is not random (imagine a scenario where I am first given higher values and then lower values). – Bitwise Apr 18 '13 at 17:12
• That is a crucial observation. It implies you need to take more care than usual, because initially you will be obtaining a "robust" estimate of the mean high outliers. By continuing to update that estimate, you could wind up throwing out all the lower values. Thus you will need a data structure in which key parts of the entire distribution of data are recorded and periodically updated. Check out our threads with keywords "online" and "quantile" for ideas. Two such promising ones are at stats.stackexchange.com/questions/3372 and stats.stackexchange.com/q/3377. – whuber Apr 18 '13 at 17:18
• I would offer a bounty but i don't have enough reputation – Jason S Apr 4 '14 at 15:24
• To continue with the idea in @whuber's first comment, you can maintain a uniformly sampled random subset of size $100$ or $1000$ from all the data seen thus far. This set and the associated "fences" can be updated in O(1) time. – Innuo Apr 22 '14 at 16:15

You might think of relating your problem to that of the recursive control chart. Such a control chart will evaluate whether a new observation is in control. If it is, this observation is included in the new estimate of the mean and variance (necessary to determine control limits).

Some background on robust, recursive, univariate control charts can be found here. One of the classic texts on quality control and control charts appears to be available online here.

Intuitively, using the a mean, $\mu_{t-1}$ and a variance $\sigma^2_{t-1}$ as inputs, you can determine whether a new observation at time $t$ is an outlier by a number of approaches. One would be to declare $x_t$ an outlier if it is outside of a certain number of standard deviations of $\mu_{t-1}$ (given $\sigma^2_{t-1})$, but this may run into problems if the data does not conform to certain distributional assumptions. If you want to go this road, then supposing you have determined if a new point is not an outlier, and would like to include it in your mean estimate with no special rate of forgetting. Then you can't do better than:

$\mu_t = \frac{t-1}{t}\mu_{t-1}+\frac{1}{t}x_t$

Similarly, you will need to update the variance recursively:

$\sigma^2_t = \frac{t-1}{t}\sigma^2_{t-1}+\frac{1}{t-1}(x_t-\mu_t)^2$

However, you might want to try some more conventional control charts. Other control charts which are more robust to the distribution of the data and can still handle non-stationarity (like the $\mu$ of your process slowly going higher) are the EWMA or CUSUM are recommended (see the textbook linked to above for more details on the charts and their control limits). These methods will typically be less computationally intensive than a robust because they have the advantage of simply needing to compare a single new observation to information derived from non-outlier observations. You can refine your estimates of the long term process $\mu$ and $\sigma^2$ used in the control limit calculations of these methods with the updating formulas given above if you like.

Regarding a chart like the EWMA, which forgets old observations and gives more weight to new ones, if you think that your data is stationary (meaning the parameters of the generating distribution do not change) then there is no need to forget older observations exponentially. You can set the forgetting factor accordingly. However, if you think that it is non-stationarity you will need to select a good value for the forgetting factor (again see the textbook for a way to do this).

I should also mention that before you begin monitoring and adding new observations online, you will need to obtain estimates of $\mu_0$ and $\sigma^2_0$ (the initial parameter values based on a training dataset) that are not influenced by outliers. If you suspect there are outliers in your training data, you can pay the one-time cost of using a robust method to estimate them.

I think an approach along these lines will lead to the fastest updating for your problem.

• Using control charts is an interesting idea. It seems that it might be difficult to overcome the challenges outlined in comments to the question, however. In the non-stationary case if you are "forgetting" older values it seems possible the estimates could be highly biased. For instance, how would your suggestions perform on a stream of data given by $x_t=\cos(\pi t/10^6)+2\lfloor t/10^6\rfloor$? (This drops very gradually, suddenly jumps and rises very gradually, suddenly jumps again, and so on.) – whuber Apr 22 '14 at 16:26
• @Bitwise says initial sample may not represent future data. Without info about how different the rest of the data will be you essentially can't do anything. However if the initial data has info on the process non-stationarity (say a downward trend) then new observations can be allowed in accounting for the fact that we expect them to be lower. However, some info about the non-stationarity is needed. You propose one pathological type of non-stationarity. Some methods, eg the EWMA are optimal for a specific process but are generally pretty good. Your process would require a more custom job. – Deathkill14 Apr 22 '14 at 17:26
• (I detect a mathematician in you, because it's a very mathematical move to dismiss as "pathological" something you can't handle :-). But I beg to differ with your prognosis: methods like those suggested by @Innuo can indeed protect against such "pathologies" and everything else the real world might throw at you, especially when randomization is incorporated into the sampling. – whuber Apr 22 '14 at 18:32
• Actually, I agree that one should not dismiss a problem that one is faced by. Could you please link me to the methods @Innuo discussed (I can't find them from this post - were they in the links you provided above and I missed them?). Thank you. – Deathkill14 Apr 22 '14 at 18:41
• @Innuo posted a brief comment at stats.stackexchange.com/questions/56494/… suggesting a uniform random sample of all previously observed data could be maintained in $O(1)$ time. Although it's not quite clear exactly how that would be done, coupling it with a robust estimator of the mean would constitute a universal solution, applicable to any stream of data whatsoever. – whuber Apr 22 '14 at 18:46

This solution implements a suggestion made by @Innuo in a comment to the question:

You can maintain a uniformly sampled random subset of size 100 or 1000 from all the data seen thus far. This set and the associated "fences" can be updated in $O(1)$ time.

Once we know how to maintain this subset, we may select any method we like to estimate the mean of a population from such a sample. This is a universal method, making no assumptions whatsoever, that will work with any input stream to within an accuracy that can be predicted using standard statistical sampling formulas. (The accuracy is inversely proportional to the square root of the sample size.)

This algorithm accepts as input a stream of data $x(t),$ $t=1, 2, \ldots,$ a sample size $m$, and outputs a stream of samples $s(t)$ each of which represent the population $X(t) = (x(1),x(2), \ldots, x(t))$. Specifically, for $1\le i \le t$, $s(i)$ is a simple random sample of size $m$ from $X(t)$ (without replacement).

For this to happen, it suffices that every $m$-element subset of $\{1,2,\ldots,t\}$ have equal chances of being the indexes of $x$ in $s(t)$. This implies the chance that $x(i),$ $1\le i\lt t,$ is in $s(t)$ equals $m/t$ provided $t \ge m$.

At the beginning we just collect the stream until $m$ elements have been stored. At that point there is only one possible sample, so the probability condition is trivially satisfied.

The algorithm takes over when $t=m+1$. Inductively suppose that $s(t)$ is a simple random sample of $X(t)$ for $t\gt m$. Provisionally set $s(t+1) = s(t)$. Let $U(t+1)$ be a uniform random variable (independent of any previous variables used to construct $s(t)$). If $U(t+1) \le m/(t+1)$ then replace a randomly chosen element of $s$ by $x(t+1)$. That's the entire procedure!

Clearly $x(t+1)$ has probability $m/(t+1)$ of being in $s(t+1)$. Moreover, by the induction hypothesis, $x(i)$ had probability $m/t$ of being in $s(t)$ when $i\le t$. With probability $m/(t+1) \times 1/m$ = $1/(t+1)$ it will have been removed from $s(t+1)$, whence its probability of remaining equals

$$\frac{m}{t}\left(1 - \frac{1}{t+1}\right) = \frac{m}{t+1},$$

exactly as needed. By induction, then, all the inclusion probabilities of the $x(i)$ in the $s(t)$ are correct and it's clear there is no special correlation among those inclusions. That proves the algorithm is correct.

The algorithm efficiency is $O(1)$ because at each stage at most two random numbers are computed and at most one element of an array of $m$ values is replaced. The storage requirement is $O(m)$.

The data structure for this algorithm consists of the sample $s$ together with the index $t$ of the population $X(t)$ that it samples. Initially we take $s = X(m)$ and proceed with the algorithm for $t=m+1, m+2, \ldots.$ Here is an R implementation to update $(s,t)$ with a value $x$ to produce $(s,t+1)$. (The argument n plays the role of $t$ and sample.size is $m$. The index $t$ will be maintained by the caller.)

update <- function(s, x, n, sample.size) {
if (length(s) < sample.size) {
s <- c(s, x)
} else if (runif(1) <= sample.size / n) {
i <- sample.int(length(s), 1)
s[i] <- x
}
return (s)
}


To illustrate and test this, I will use the usual (non-robust) estimator of the mean and compare the mean as estimated from $s(t)$ to the actual mean of $X(t)$ (the cumulative set of data seen at each step). I chose a somewhat difficult input stream that changes quite smoothly but periodically undergoes dramatic jumps. The sample size of $m=50$ is fairly small, allowing us to see sampling fluctuations in these plots.

n <- 10^3
x <- sapply(1:(7*n), function(t) cos(pi*t/n) + 2*floor((1+t)/n))
n.sample <- 50
s <- x[1:(n.sample-1)]
online <- sapply(n.sample:length(x), function(i) {
s <<- update(s, x[i], i, n.sample)
summary(s)})
actual <- sapply(n.sample:length(x), function(i) summary(x[1:i]))


At this point online is the sequence of mean estimates produced by maintaining this running sample of $50$ values while actual is the sequence of mean estimates produced from all the data available at each moment. The plot shows the data (in gray), actual (in black), and two independent applications of this sampling procedure (in colors). The agreement is within expected sampling error:

plot(x, pch=".", col="Gray")
lines(1:dim(actual)[2], actual["Mean", ])
lines(1:dim(online)[2], online["Mean", ], col="Red")


For robust estimators of the mean, please search our site for and related terms. Among the possibilities worth considering are Winsorized means and M-estimators.

• it's not clear to me how the rejection threshold look like in this approach (e.g. the threshold beyond which observations are rejected as outliers). Can you add them to the plot? – user603 Apr 23 '14 at 6:46
• @user603 The "rejection threshold," or whatever robust method is used to estimate the mean, is irrelevant: choose whatever method you wish to estimate the mean. (Not all robust methods work by erecting thresholds and rejecting data, BTW.) This would be done in the code of my answer by replacing summary by a robust variant. – whuber Apr 23 '14 at 13:56
• Something is not clear to me in this example. Are the grey data "good" or "outliers." If the prior, it seems the fit is biased down (it should fit them better since the situation would be similar to @Bitwise's downward trend we would like to follow). If the grey data at higher index values is outliers, then it seems the fit is biased upwards. What is the target you want to fit here? The current fit seems torn between those two scenarios. – Deathkill14 Apr 23 '14 at 20:42
• @Death As explained in the text immediately preceding the figure, the gray data are the original stream of data. Its running mean is the black curve. The colored curves are based on the algorithm. Vertical deviations of the colored curves relative to the black curve are due to the randomness in the sampling. The expected amount of deviation at any index is proportional to the standard deviation of the gray values preceding that index and inversely proportional to the square root of the sample size (taken as 50 in this example). – whuber Apr 24 '14 at 15:07