Online algorithm for mean absolute deviation and large data set I have a little problem that is making me freaking out.
I have to write procedure for an online acquisition process of a multivariate time series.
At every time interval (for example 1 second), I get a new sample, which is basically a floating point vector of size N.
The operation I need to do is a little bit tricky:


*

*For each new sample, I compute the percentuals for that sample (by normalizing the vector so that the elements will sum to 1).

*I calculate the average percentuals vector in the same way, but using the past values.

*For each past value, I compute the absolute deviation of the percentuals vector related to that sample with the global average percentuals vector computed at step 2. This way, the absolute deviation is always a number between 0 (when the vector is equal to the average vector) and 2 (when it is totaly different).

*Using the average of the deviations for all the previous samples, I compute the mean absolute deviation, which is again a number between 0 and 2.

*I use the mean absolute deviation to detect if a new sample is compatible with the other samples (by comparing its absolute deviation with the mean absolute deviation of the whole set computed at step 4).
Since every time a new sample is collected the global average changes (and so the mean absolute deviation changes as well), is there a way to compute this value without scanning the whole data set multiple times? (one time for computing the global average percentuals, and one time for collecting the absolute deviations).
Ok, I know it's absolutely easy to calculate the global averages without scanning the whole set, since I just have to use a temporary vector for storing the sum of each dimension, but what about the mean absolute deviation? Its calculation includes the abs() operator, so I need to access to all the past data!
Thanks for your help.
 A: If you can accept some inaccuracy, this problem can be solved easily by binning counts. That is, pick some largeish number $M$ (say $M = 1000$), then initialize some integer bins $B_{i,j}$ for $i = 1\ldots M$ and $j = 1\ldots N$, where $N$ is the vector size, as zero. Then when you see the $k$th observation of a percentual vector, increment $B_{i,j}$ if the $j$th element of this vector is between $(i-1)/M$ and $i/M$, looping over the $N$ elements of the vector. (I am assuming your input vectors are non-negative, so that when you compute your 'percentuals', the vectors are in the range $[0,1]$. )
At any point in time, you can estimate the mean vector from the bins, and the mean absolute deviation. After observing $K$ such vectors, the $j$th element of the mean is estimated by 
$$\bar{X}_j = \frac{1}{K} \sum_i \frac{i - 1/2}{M} B_{i,j},$$ and the $j$th element of the mean absolute deviation is estimated by $$\frac{1}{K} \sum_i | \bar{X_j} - \frac{i - 1/2}{M} | B_{i,j}$$
edit: this is a specific case of a more general approach where you are building an empirical density estimate. This could be done with polynomials, splines, etc, but the binning approach is the easiest to describe and implement. 
A: I've used the following approach in the past to calculate absolution deviation moderately efficiently (note, this a programmers approach, not a statisticians, so indubitably there may be clever tricks like shabbychef's that might be more efficient). 
WARNING: This is not an online algorithm.  It requires O(n) memory.  Furthermore, it has a worst case performance of O(n), for datasets like [1, -2, 4, -8, 16, -32, ...] (i.e. the same as the full recalculation). [1] 
However, because it still performs well in many use cases it might be worth posting here.  For example, in order to calculate the absolute deviance of 10000 random numbers between -100 and 100 as each item arrives, my algorithm takes less than one second, while the full recalculation takes over 17 seconds (on my machine, will vary per machine and according to input data).  You need to maintain the entire vector in memory however, which may be a constraint for some uses.  The outline of the algorithm is as follows:


*

*Instead of having a single vector to store past measurements, use three sorted priority queues (something like a min/max heap).  These three lists partition the input into three: items greater than the mean, items less than the mean and items equal to the mean.

*(Almost) every time you add an item the mean changes, so we need to repartition.  The crucial thing is the sorted nature of the partitions which means that instead of scanning every item in the list to repartion, we only need to read those items we are moving.  While in the worst case this will still require O(n) move operations, for many use-cases this is not so.

*Using some clever bookkeeping, we can make sure that the deviance is correctly calculated at all times, when repartitioning and when adding new items.


Some sample code, in python, is below.  Note that it only allows items to be added to the list, not removed.  This could easily be added, but at the time I wrote this I had no need for it.  Rather than implement the priority queues myself, I have used the sortedlist from Daniel Stutzbach's excellent blist package, which use B+Trees internally.
Consider this code licensed under the MIT license.  It has not been significantly optimised or polished, but has worked for me in the past.  New versions will be available here.  Let me know if you have any questions, or find any bugs.
from blist import sortedlist
import operator

class deviance_list:
    def __init__(self):
        self.mean =  0.0
        self._old_mean = 0.0
        self._sum =  0L
        self._n =  0  #n items
        # items greater than the mean
        self._toplist =  sortedlist()
        # items less than the mean
        self._bottomlist = sortedlist(key = operator.neg)
        # Since all items in the "eq list" have the same value (self.mean) we don't need
        # to maintain an eq list, only a count
        self._eqlistlen = 0

        self._top_deviance =  0
        self._bottom_deviance =  0

    @property
    def absolute_deviance(self):
        return self._top_deviance + self._bottom_deviance

    def append(self,  n):
        # Update summary stats
        self._sum += n
        self._n +=  1
        self._old_mean =  self.mean
        self.mean =  self._sum /  float(self._n)

        # Move existing things around
        going_up = self.mean > self._old_mean
        self._rebalance(going_up)

        # Add new item to appropriate list
        if n >  self.mean:
            self._toplist.add(n)
            self._top_deviance +=  n -  self.mean
        elif n == self.mean: 
            self._eqlistlen += 1
        else:
            self._bottomlist.add(n)
            self._bottom_deviance += self.mean -  n


    def _move_eqs(self,  going_up):
        if going_up:
            self._bottomlist.update([self._old_mean] *  self._eqlistlen)
            self._bottom_deviance += (self.mean - self._old_mean) * self._eqlistlen
            self._eqlistlen = 0
        else:
            self._toplist.update([self._old_mean] *  self._eqlistlen)
            self._top_deviance += (self._old_mean - self.mean) * self._eqlistlen
            self._eqlistlen = 0


    def _rebalance(self, going_up):
        move_count,  eq_move_count = 0, 0
        if going_up:
            # increase the bottom deviance of the items already in the bottomlist
            if self.mean !=  self._old_mean:
                self._bottom_deviance += len(self._bottomlist) *  (self.mean -  self._old_mean)
                self._move_eqs(going_up)


            # transfer items from top to bottom (or eq) list, and change the deviances
            for n in iter(self._toplist):
                if n < self.mean:
                    self._top_deviance -= n -  self._old_mean
                    self._bottom_deviance += (self.mean -  n)
                    # we increment movecount and move them after the list
                    # has finished iterating so we don't modify the list during iteration
                    move_count +=  1
                elif n == self.mean:
                    self._top_deviance -= n -  self._old_mean
                    self._eqlistlen += 1
                    eq_move_count +=  1
                else:
                    break
            for _ in xrange(0,  move_count):
                self._bottomlist.add(self._toplist.pop(0))
            for _ in xrange(0,  eq_move_count):
                self._toplist.pop(0)

            # decrease the top deviance of the items remain in the toplist
            self._top_deviance -= len(self._toplist) *  (self.mean -  self._old_mean)
        else:
            if self.mean !=  self._old_mean:
                self._top_deviance += len(self._toplist) *  (self._old_mean -  self.mean)
                self._move_eqs(going_up)
            for n in iter(self._bottomlist): 
                if n > self.mean:
                    self._bottom_deviance -= self._old_mean -  n
                    self._top_deviance += n -  self.mean
                    move_count += 1
                elif n == self.mean:
                    self._bottom_deviance -= self._old_mean -  n
                    self._eqlistlen += 1
                    eq_move_count +=  1
                else:
                    break
            for _ in xrange(0,  move_count):
                    self._toplist.add(self._bottomlist.pop(0))
            for _ in xrange(0,  eq_move_count):
                self._bottomlist.pop(0)

            # decrease the bottom deviance of the items remain in the bottomlist
            self._bottom_deviance -= len(self._bottomlist) *  (self._old_mean -  self.mean)


if __name__ ==  "__main__":
    import random
    dv =  deviance_list()
    # Test against some random data,  and calculate result manually (nb. slowly) to ensure correctness
    rands = [random.randint(-100,  100) for _ in range(0,  1000)]
    ns = []
    for n in rands: 
        dv.append(n)
        ns.append(n)
        print("added:%4d,  mean:%3.2f,  oldmean:%3.2f,  mean ad:%3.2f" %
              (n, dv.mean,  dv._old_mean,  dv.absolute_deviance / dv.mean))
        assert sum(ns) == dv._sum,  "Sums not equal!"
        assert len(ns) == dv._n,  "Counts not equal!"
        m = sum(ns) / float(len(ns))
        assert m == dv.mean,  "Means not equal!"
        real_abs_dev = sum([abs(m - x) for x in ns])
        # Due to floating point imprecision, we check if the difference between the
        # two ways of calculating the asb. dev. is small rather than checking equality
        assert abs(real_abs_dev - dv.absolute_deviance) < 0.01, (
            "Absolute deviances not equal. Real:%.2f,  calc:%.2f" %  (real_abs_dev,  dv.absolute_deviance))

[1] If symptoms persist, see your doctor.
A: There is also a parametric approach. Ignoring the vector nature of your data, and looking only at the marginals, it suffices to solve the problem: find an online algorithm to compute the mean absolute deviation of scalar $X$. If (and this is the big 'if' here) you thought that $X$ followed some probability distribution with unknown parameters, you could estimate the parameters using some online algorithm, then compute the mean absolute deviation based on that parametrized distribution. For example, if you thought that $X$ was (approximately) normally distributed, you could estimate its standard deviation, as $s$, and the mean absolute deviation would be estimated by $s \sqrt{2 / \pi}$ (see Half Normal Distribution).
A: MAD(x) is just two concurrent median computation, each of which can be made online through the binmedian algorithm.
You can find the associated paper as well as C and FORTRAN code online here.
(this is just the use of a clever trick on top of Shabbychef's clever trick, to save on memory).
Addendum:
There exist a host of older multi-pass methods for computing quantiles. A popular approach is to maintain/update a deterministically sized reservoir of observations randomly selected from the stream and recursively compute quantiles (see this review) on this reservoir. This (and related) approach are superseded by the one proposed above.
A: The following provides an inaccurate approximation, although the inaccuracy will depend on the distribution of the input data.  It is an online algorithm, but only approximates the absolute deviance.  It is based on a well known algorithm for calculating variance online, described by Welford in the 1960s.  His algorithm, translated into R, looks like:
M2 <- 0
mean <- 0
n <- 0

var.online <- function(x){
    n <<- n + 1
    diff <- x - mean
    mean <<- mean + diff / n
    M2 <<- M2 + diff * (x - mean)
    variance <- M2 / (n - 1)
    return(variance)
}

It performs very similarly to R's builtin variance function:
set.seed(2099)
n.testitems <- 1000
n.tests <- 100
differences <- rep(NA, n.tests)
for (i in 1:n.tests){
        # Reset counters
        M2 <- 0
        mean <- 0
        n <- 0

        xs <- rnorm(n.testitems)
        for (j in 1:n.testitems){
                v <- var.online(xs[j])
        }

        differences[i] <- abs(v - var(xs))

}
summary(differences)
     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
0.000e+00 2.220e-16 4.996e-16 6.595e-16 9.992e-16 1.887e-15 

Modifying the algorithm to calculate absolute deviation simply involves an additional sqrt call.  However, the sqrt introduces inaccuracies that are reflected in the result:
absolute.deviance.online <- function(x){
    n <<- n + 1
    diff <- x - mean
    mean <<- mean + diff / n
    a.dev <<- a.dev + sqrt(diff * (x - mean))
    return(a.dev)
}

The errors, calculated as above, are much greater than for the variance calculation:
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
0.005126 0.364600 0.808000 0.958800 1.360000 3.312000 

However, depending on your use case, this magnitude of error might be acceptable.

