Exponentially weighted averaging for rates

In a previous post, I talked about how to produce exponentially weighted averages of time-embedded values.  This is fine for cases where you really want averages, but it doesn't work for rates.  Happily, the technique I presented for simple averages can be easily extended to work for rates as well.

To recap a bit, an on-line algorithm for computing a simple average works by accumulating the number of observations and the sum of the observations. There is no mention of time.  For each observation $(x, t)$ we update our state $(s, w)$ according to
\begin{eqnarray*}
s_{n} &=& s_{n-1} + x \\
w_{n} &=& w_{n-1} + 1
\end{eqnarray*}
The exponentially weighted average approach that I mentioned before uses the time of each observation to discount the previous state based on the time of the previous observation.  Thus, the state $(s,w,t)$ is updated according to
\begin{eqnarray*}
\pi &=& e^{-(t-t_{n-1})/\alpha} \\
s_{n} &=& \pi s_{n-1} + x \\
w_{n} &=& \pi w_{n-1} + 1 \\
t_{n} &=& t
\end{eqnarray*}
If we were to compute simple rates without discounting, then instead of interpreting $w$ as a count, we would interpret it as a time interval. Thus we would update state $(s,w,t)$ according to with an observation $(x, t)$ according to
\begin{eqnarray*}
s_{n} &=& s_{n-1} + x \\
w_{n} &=& w_{n-1} + (t - t_n) \\
t_n &=& t
\end{eqnarray*}
By analogy with the discounted version of the simple average, we can produce an exponentially weighted rate average by updating the state according to this
\begin{eqnarray*}
\pi &=& e^{-(t-t_{n-1})/\alpha} \\
s_{n} &=& \pi s_{n-1} + x \\
w_{n} &=& \pi w_{n-1} + (t-t_n) \\
t_{n} &=& t
\end{eqnarray*}
This approach has a defect in that each data pont should really be considered to be a report of events spread out in an uncertain way over the time since the last report. As such, the interval and the events should probably both be discounted as if they had been reported uniformly over the period in question. In practice, this correction does not matter much since the averaging time constant $\alpha$ is typically large compared to the average reporting interval. Another consideration comes up when multiple sources are reporting concurrently. If we want to do proper discounting of each observed number, then we have to keep track of the time since last report for each of the sources. Since this probably won't matter and since it considerably complicates the implementation, I would rather not do it.

See https://issues.apache.org/jira/browse/MAHOUT-634 for code. This will be in Mahout before long.

Update on exponential time-embedded averages

I will be adding this code to Mahout shortly.

See https://issues.apache.org/jira/browse/MAHOUT-634 for status.

Also, if you are measuring rates, then it is common for rates to be reported from multiple sources independently.  Such an average can be computed pretty easily using this same framework if the sources report often relative to the averaging time constant.  This simple implementation just attributes each reported count as if they occurred in the interval since the most recent report from any reporter.  If the time constant is relatively long, this can work out reasonably well as long as we are careful.

If reporting intervals are longer, then the averaging is a bit trickier because we really would like to attribute the reported counts over the entire interval from the last report from the same source.  This means that we have to discount some of the counts because they are effectively kind of old.

More details shortly.

Exponential weighted averages with irregular sampling

Ted Yu on the hbase list wants to compute an exponential weighted average for the number of queries per second or other statistics that characterize the performance or state of an hbase system.

The trick is that the measurements are only available at irregular intervals. If they were sampled regularly, then the standard mixing trick would work:
\[
m_{n+1} = \mu m_n + (1-\mu) x_{n+1}
\]
where $m$ is our current estimate of the mean, $x_n$ is the $n$-th sample and $\mu$ determines how much history to use.

With unequal sample times, things become a bit more complicated. If we get lots of measurements all at once, we want to give them nearly equal weight but if we have a long gap, we want to weight the very old samples much less.

In fact, we want to weight old samples according to how old they are with exponentially decreasing weight. If we sample values $\left \lbrace x_1 \ldots x_n \right \rbrace$ at times $t_1 \ldots t_n$ then we want the weighted mean defined by
\[
m_n = {\sum_{i=1}^n x_i e^{-(t_n - t_i)/\alpha} \over \sum_{i=1}^n e^{-(t_n - t_i)/\alpha} }
\]
Here $\alpha$ plays the same role as $\mu$ did before, but on a different scale. If the evenly sampled data comes at time intervals $\Delta t$ then $\mu = e^{\Delta t / \alpha}$.

Happily, there is a very simple recurrence relationship that allows us to keep only two intermediate values while computing the value of $m_1 \ldots m_n$ in an entirely on-line fashion as the $x_i$ values arrive.

To see this, define

\begin{eqnarray*}
\pi_n &=& e^{-(t_{n+1}-t_n)/\alpha} \\
w_{n+1} &=&
\sum_{i=1}^{n+1} e^{-(t_{n+1} - t_i)/\alpha} =
1+e^{-(t_{n+1}-t_n)/\alpha} \sum_{i=1}^{n} e^{-(t_{n} - t_i)/\alpha} \\
& =& 1 + \pi w_n\\
s_{n+1} &=&
\sum_{i=1}^{n+1} x_i e^{-(t_{n+1} - t_i)/\alpha} =
x_{n+1}+e^{-(t_{n+1}-t_n)/\alpha} \sum_{i=1}^{n} x_i e^{-(t_{n} - t_i)/\alpha} \\
&=& x_{n+1} + \pi_n s_n
\end{eqnarray*}


Then note that
\[
m_{n+1} = {s_{n+1} \over w_{n+1}}
\]

This leads naturally to a procedure that has state consisting of $t, w, m$ which are updated with using new values of $t_n, x_n$ according to
\begin{eqnarray*}
\pi &=& e^{-(t_{n}-t)/\alpha} \\
w &=& 1 + \pi w \\
s &=& x_n + \pi s \\
m &=& {s \over w} \\
t &=& t_{n}
\end{eqnarray*}
Isn't that a kick!

To do this right, however, we need a test. Here are some data vectors computed for $\alpha=5$:


         t          x        pi        w          s          m
1 11.35718  1.5992071 1.0000000 1.000000  1.5992071  1.5992071
2 21.54637 -1.3577032 0.1303100 1.130310 -1.1493105 -1.0168100
3 28.91061 -0.3405638 0.2292718 1.259148 -0.6040683 -0.4797436
4 33.03586  0.7048632 0.4382129 1.551775  0.4401527  0.2836447
5 39.57767  0.3020558 0.2702621 1.419386  0.4210124  0.2966159

Writing a proper unit test is an exercise left to the reader. (but I will be happy to help)

Recent lecture on Mahout for SDForum

I gave a lecture last night on recent additions to Mahout.

The slides are here: http://www.slideshare.net/tdunning/sdforum-11042010

I can amplify this post with answers to any questions that anybody puts in the comments.

Mahout 0.4 released!

Go to the Apache Mahout site for more info.  Here is the official announcement:

We are pleased to announce release 0.4 of Mahout. Virtually every corner of the project has changed, and significantly, since 0.3. Developers are invited to use and depend on version 0.4 even as yet more change is to be expected before the next release. Highlights include:

- Model refactoring and CLI changes to improve integration and consistency

- New ClusterEvaluator and CDbwClusterEvaluator offer new ways to evaluate clustering effectiveness

- New Spectral Clustering and MinHash Clustering (still experimental)

- New VectorModelClassifier allows any set of clusters to be used for classification

- Map/Reduce job to compute the pairwise similarities of the rows of a matrix using a customizable similarity measure

- Map/Reduce job to compute the item-item-similarities for item-based collaborative filtering

- RecommenderJob has been evolved to a fully distributed item-based recommender

- Distributed Lanczos SVD implementation

- More support for distributed operations on very large matrices

- Easier access to Mahout operations via the command line

- New HMM based sequence classification from GSoC (currently as sequential version only and still experimental)

- Sequential logistic regression training framework

- New SGD classifier

- Experimental new type of NB classifier, and feature reduction options for existing one

- New vector encoding framework for high speed vectorization without a pre-built dictionary

- Additional elements of supervised model evaluation framework

- Promoted several pieces of old Colt framework to tested status (QR decomposition, in particular)

- Can now save random forests and use it to classify new data

- Many, many small fixes, improvements, refactorings and cleanup

Details on what's included can be found in the release notes.

Downloads are available from the Apache Mirrors

New Mahout release coming

The vote has started for the 0.4 Mahout release.  Lots of new stuff, but the part that I am excited about is a fairly comprehensive implementation of logistic regression suitable for large scale training and high speed classification, but there is a whole lot more.

With the 0.4 release, Mahout is moving along strongly towards the fabled 1.0 release.  At that point, we will start paying lots of attention to backwards compatibility.  That will be good, but the current wild and wooly policy is pretty handy if you have something in mind that Mahout really, really needs because we can get new things in pretty readily right now.

See http://mahout.apache.org/ for the release when it arrives and watch my twitter feed for an announcement.

Why is the sum of two uniform randoms not uniform?

Lance Norkog asks on the Mahout mailing list why adding two uniformly distributed random variables gives a pyramidal distributed value. I would normally answer on the mailing list, but here I can use lovely math notation. As I mentioned on-list, this is a very basic result that is related to the law of large numbers.

If we were to draw a picture of the joint distribution of these variables \(x\) and \(y\), we would get something that is 1 inside the \([0,1] \times [0,1]\) square and 0 outside that region.

For a given value \(\alpha\) of the sum \(x + y\), there is a diagonal line segment where \(x+y=\alpha\) and \(x\) and \(y\) are in the square. Where \(z \le 0\) or \(z \ge 2\) that intersection vanishes and for \(0 \lt z \lt 2\), that intersection varies in length. The probability of the sum having some particular value z is proportional to the length of that intersection. As you can imagine, the intersection varies in size linearly and it reaches a maximum where z = 1.

For the sum of three random variables, the situation is more complex to reason about geometrically because we need to worry about the intersection of a plane and a cube.  For more variables, the geometry is not worth the trouble.

If we tackle the problem a bit more rigorously, then the easiest way to approach the problem is to compute the cumulative distribution of various values of sums. That leads to a convolution integral over the density functions involved. Since the densities are all 1, the integration limits are the key to the value and those limits have to broken down into cases. Actually doing these integrals is a pretty rare activity since the limit is approximated so well after just a few iterations.

Just how quickly that convergence happens is can be seen by looking at the empirical distribution of the sum of three uniform deviates.  I used something very like this R code to produce the graph below:


   hist(runif(100000)+runif(100000)+runif(100000), 
        breaks=50, prob=T)
   lines(seq(-1,4,by=0.01), dnorm(seq(-1,4,by=0.01), 
        1.5, sqrt(1/4)))


In this graph, the red curve is the normal distribution with the same mean and standard deviation.  As you can see, the peak is a tiny bit too high and the tails are just a skoshi too long for the normal to be a good description of the samples of the sum.   

This is, however, just the sum of three random values.  If you sum more values, the convergence to the normal distribution is very strong and by the time you are adding six uniform random values together, the difference between the distributions is no longer visible in a graph like this and can only be detected numerically using lots of data and clever things like a Kolmogorov-Smirnov test.

The moral here is that there isn't much way to avoid this regression to the normal distribution and distorting the data to avoid it is probably pointless.

But if you are like me, being just a little more normal always made it easier to get along.