Estimating rate of occurrence of an event with exponential smoothing and irregular events

Question

Imagine that I have a set of measurements of x that are taken by many processes x₀ ... x_N at times t₀ ... t_N. Let's assume that at time t I want to make an estimate of the current value of x, based on the assumption that there is no long term trend I know about and that x can be predicted from an algorithm such as exponential smoothing. As we have many processes, and N can get very large, I can't store more than a few values (e.g. the previous state).

One approach here would be to adapt the normal exponential smoothing algorithm. If samples are taken regularly, I would maintain an estimator y_n such that:

y_n = α.y_n-1 + ( 1 - α ). x_n

This approach is not great where the sampling is irregular as many samples together would have a disproportionate influence. Therefore this formula could be adapted to:

y_n = α_n.y_n-1 + ( 1 - α_n ). x_n

where

α_n = e^{-k.(t_n - t_n-1)}

IE dynamically adjusting the smoothing constant depending on the interval between the previous two samples. I'm happy with this approach and it seems to work. It's the first answer given here, and a good summary of these sorts of techniques are given by Eckner in this 2012 paper (PDF).

Now, my question is as follows. I want to adapt the above to estimate the rate of an occurrence. Occasionally an event will occur. Using a similar exponential technique, I want to get an estimate of the rate that event occurs.

Two obvious strategies would be:

• To use the first or the second technique using the delay between the last two events as the data series x_n.
• To use the first or the second technique using the reciprocal of the delay between the last two events (i.e. the rate) as the data series x_n.

Neither of these turn out to be good strategies as far as I can tell. Firstly, take an event that occurs every 500ms (on the one hand) and an event that occurs with a 200ms delay and an 800ms delay on the other. Clearly these both occur twice a second, so the rate estimate given should be the same. Ignoring the time from the last sample seems foolhardy, so I'll concentrate on the second strategy. Using the delay (rather than the reciprocal) does not turn out to be a good predictor because the simulating the 200ms/800ms sample stream produces an estimate of about 1.5 (on the basis the average of reciprocals is not the reciprocal of the average).

But, far more significantly, neither strategy copes with what is actually happening in practice, which is that suddenly all the events stop for a long while. The 'latest' value of y is thus the value at the last event, and the estimate of the rate thus never gets calculated. Hence the rate appears constant. Of course if I was analysing the data retrospectively, this wouldn't be an issue, but I'm analysing it in real time.

I realise another way to do this would be to run some thread periodically (e.g. every 10 seconds) and count the number of occurrences in this 10 second interval. This is very resource heavy my end as the statistics are not needed often, and I am loathe to run a thread that polls everything due to mutex issues. Therefore I'd like to (somehow) use an algorithm which adjusts the state read by (e.g.) the time since the last sample was taken. This seems a reasonable approach as if the performance is measured at times chosen independently of the samples, the measurement time will on average be half way through the period between samples, so a very crude unsmoothed estimate of the rate would be half the reciprocal of the time since the last sample. To complicate things further, my measurement time will not be independent of the samples.

I have a feeling this has a simple answer but it's eluding me. I have a feeling that the correct route is to assume the the events are Poisson distributed, and derive an estimate for λ based on the interval since the last sample and some form of moving average, but my statistics is too rusty to make this work.

There is a near dupe of this question here but the answer doesn't seem to be very satisfactory (I hope I explained why). I'd add that a Kalman filter seems way to heavyweight given I have one variable to estimate and know nothing about it. There are a number of other near dupes, most of which either suggest keeping large bins of values (not realistic here from a memory point of view) or do not address the two issues above.

Gunboat

Can you keep a list of the say last 10 or so time stamps the event occurred? – Jonjona

How do the statistics of your events look? The reason I'm saying this is that spooky-long periods over which "all the events stop" can often indicate pathological statistics, like power-law statistics, for which these estimators never work well. – Necessarily

@DanielBrückner I could keep a list of a (fixed) small number of events (e.g. 10), and say "how long did those take to occur". However, that wouldn't help with the issue where no more events come (at all) and I need to estimate the rate. – Gunboat

@RamanShah the events should be come in at a roughly constant (probably Poisson distributed in the short term) rate from each of a number of sources, leading to a nice steady average. What I'm trying to measure is (a) big drop offs in these indicating severe connectivity problems, (b) a reasonably accurate 'rate' in steady state conditions so there is something to measure performance against. You could consider it a bit like interface counters on networking equipment (with many many interfaces). – Gunboat

So to clarify, you want your estimator to decay to zero when you haven't seen an event for a long time? I think I'm starting to understand your data and what you want from it, but I'm not certain. – Necessarily

Here is a paper dealing with your situation of estimating the parameter of an inhomogeneous Poisson process. They also mention the problem of the rate going to zero. – Jonjona

@DanielBrückner thanks. Looks very useful. I shall read that. In the meantime you gave me an idea which I shall post below. – Gunboat

Once you start looking for "estimate inhomogeneous poisson process` countless papers occur. This seems a well studied problem. There is even at least one library. – Jonjona

@DanielBrückner unless I'm misreading that paper, it appears to require that you store the timestamps of all the data. Are you suggesting that you apply it merely to the timestamps of the past N occurrences? Doesn't that mean the 'window' of the moving average would contract when the data becomes more frequent? My issue is the data could be much less than once a second, or tens of times a second, and I'd quite like the averages to have similar-ish time periods. – Gunboat

I made my first comment before I discovered the paper and I did not read it entirely. My idea was to keep the last n events and when you want to know the current rate you calculate it from the last n occurrences and the time since the last occurrence. So you will not have a variable containing the current rate but you will calculate it when needed and maybe cache the result for some time in case it is frequently used. – Jonjona

You wrote: Ignoring the time from the last sample seems foolhardy. This is actually a valid method: est:=(est*exp(-dt/T)+dt)/(exp(-dt/T)+1), where est is the estimate of the period (inverse of frequency), T the window width, and dt the time since the previous event. You can correct the bias for this one, too. – Scarecrow

Answer 1

First, if you assume that the occurrence rate of the events itself is constant (or that you're only interested in its long-term average), then you can simply estimate it as:

        λ* = N / (t − t₀)

where t is the current time, t₀ is the start of observations, N is the number of events observed since t₀ and λ* is the estimate of the true frequency λ.

At this point, it's useful to note that the estimation formula given above may be reformulated as the integral:

        λ* = integral( δ_event(τ) dτ ) / integral( 1 dτ )

where the variable of integration τ ranges from t₀ to t, and δ_event(τ) = sum( δ(τ − t_i), i = 1 .. N ) is a sum N of Dirac delta functions, with a single delta-peak at the occurrence time t_i of each event i.

Of course, this would be a completely useless way to calculate λ*, but it turns out to be a conceptually useful formulation. Basically, the way to view this formula is that the function δ_event(τ) measures the instantaneous rate at which the number of events increases at time τ, while the second integrand, which is just the constant 1, measures the rate at which time increases over time (which, of course, is simply one second per second).

---

OK, but what if the frequency λ itself may change over time, and you want to estimate its current value, or at least its average over a recent period?

Using the ratio-of-integrals formulation given above, we can obtain such an estimate simply by weighing both integrands by some weighing function w(τ) which is biased towards recent times:

        λ*_recent = integral( δ_event(τ) w(τ) dτ ) / integral( w(τ) dτ )

Now, all that remains is to pick a reasonable w(τ) such that these integrals simplify to something easy to calculate. As it turns out, if we choose an exponentially decaying weighing function of the form w(τ) = exp(k(τ − t)) for some decay rate k, the integrals simplify to:

        λ*_recent = sum( exp(k(t_i − t)), i = 0 .. N ) k / ( 1 − exp(k(t₀ − t)) )

In the limit as t₀ → −∞ (i.e., in practice, when the total observation time (t − t₀) is much larger than the weight decay timescale 1/k), this further simplifies to just:

        λ*_recent = k sum( exp(k(t_i − t)), i = 0 .. N )

Alas, naïvely applying this formula would still require us to remember all the event times t_i. However, we can use the same trick as for calculating usual exponentially weighted averages — given the weighted average event rate λ*_recent(t') at some earlier time t', and assuming that no new events have occurred between t' and t, we can calculate the current weighted average event rate λ*_recent(t) simply as:

        λ*_recent(t) = exp( k(t' − t) ) λ*_recent(t')

Further, if we now observe a new event occurring at exactly time t, the weighted average event rate just after the event becomes:

        λ*_recent(t) = k + exp( k(t' − t) ) λ*_recent(t')

---

Thus, we get a very simple rule: all we need to store is the time t_last of the previous observed event, and the estimated recent event rate λ*_last just after said event. (We may initialize these e.g. to t_last = t₀ and λ*_last = 0; in fact, with λ*_last = 0, the value of t_last makes no difference, although for non-zero λ*_last it does.)

Whenever a new event occurs (at time t_new), we update these values as:

        λ*_last ← k + exp( k(t_last − t_new) ) λ*_last
        t_last ← t_new

and whenever we wish to know the recent event rate average at the current time t, we simply calculate it as:

        λ*(t) = exp( k(t_last − t) ) λ*_last

---

Ps. To correct for the initial bias towards the (arbitrary) initial value of t_last, we can add back the 1 / ( 1 − exp(k(t₀ − t)) ) correction term that we simplified out earlier when we assumed that t ≫ t₀. To do that, simply start from t_last = 0 at t = t₀, update t_last as above, but calculate the estimated recent event rate average at time t as:

        λ*_corr(t) = exp( k(t_last − t) ) λ*_last / ( 1 − exp(k(t₀ − t)) )

(Here, t₀ denotes the time at which you start measuring events, not the occurrence of the first event.)

This will eliminate the initial bias towards zero, at the cost of increasing the early variance. Here's an example plot showing the effects of the correction, for k = 0.1 and a true mean event rate of 2:

The red line shows λ*(t) without the initial bias correction (starting from λ*(t₀) = 0), while the green line shows the bias-corrected estimate λ*_corr(t).

---

Pps. As the plot above shows, λ*, as calculated above, will not a be continuous function of time: it jumps up by k whenever an event occurs, and decays exponentially towards zero when events do not occur.

If you'd prefer a smoother estimate, you can calculate an exponentially decaying average of λ* itself:

        λ**(t) = integral( λ*(τ) exp(k₂(τ − t)) dτ ) / integral( exp(k₂(τ − t)) dτ )

where λ* is the exponentially decaying average event rate as calculated above, k₂ is the decay rate for the second average, and the integrals are over −∞ < τ ≤ t.

This integral can also be calculated by a step-wise update rule as above:

        λ**_last ← W(Δt) λ*_last + exp( −k₂ Δt ) λ**_last
        λ*_last ← k₁ + exp( −k₁ Δt ) λ*_last
        t_last ← t_new

where k₁ and k₂ are the decay rates for the first and second averages, Δt = t_new − t_last is the elapsed time between the events, and:

        W(Δt) = k₂ ( exp( −k₂ Δt ) − exp( −k₁ Δt ) ) / (k₁ − k₂)

if k₁ ≠ k₂, or

        W(Δt) = k Δt exp( −k Δt )

if k₁ = k₂ = k (the latter expression arising from the former as the limit when (k₁ − k₂) → 0).

To calculate the second average for an arbitrary point in time t, use the same formula:

        λ**(t) = W(Δt) λ*_last + exp( −k₂ Δt ) λ**_last

except with Δt = t − t_last.

---

As above, this estimate can also be bias-corrected by applying a suitable time-dependent scaling factor:

        λ**_corr(t) = λ**(t) / (1 - S(t − t₀))

where:

        S(Δt) = ( k₁ exp( −k₂ Δt ) − k₂ exp( −k₁ Δt ) ) / (k₁ − k₂)

if k₁ ≠ k₂, or

        S(Δt) = (1 + k Δt) exp( −k Δt )

if k₁ = k₂ = k.

The plot below shows the effects of this smoothing. The red and green lines show λ*(t) and λ*_corr(t) as above, while the yellow and blue lines show λ**(t) and λ**_corr(t), as calculated with k₁ = 0.1 (as above) and k₂ = 0.2:

Culottes

Fantastic; thanks. I'd pretty much got to the exponential 'scale down' of the previous estimate on each re-estimate, but had missed the way to scale that for later time <i>t</i>. – Gunboat

Yup, lovely! I've used this ratio-of-integrals formulation before for physics data processing, but the efficient updating and extraction of the average on demand is really nice. – Necessarily

Just to add, I put this in the real code (rather than just a simulator) and it works really well. My only adaption was to start with a slightly larger k and reduce it over the first 50 or so samples to the desired value. This gets the estimate closer (but no doubt with a larger SD) earlier in the run, as I have no good initial estimate of lamda to use initially. – Gunboat

@abligh: You can actually fix the initial bias simply by scaling the estimate by 1 / (1 − exp(k(t0 − t))); see the edit I just made above. – Culottes

@IlmariKaronen thanks. I need to think about how to use that because some of my sample streams aren't necessarily running at all when the program starts, so for them t0 should be when the program starts, not the first sample, as clearly if the program runs for 10 minutes before the first sample, then gets 2 samples in one second, the best estimate of the rate is not 2! – Gunboat

@abligh: Yes, t0 should be the time at which you start counting events; the first actual event will generally occur some time later, at t1. (Ps. I added some extra detail to my answer above about using second-order smoothing to get rid of the jumps in the average when an event occurs.) – Culottes

amazing answer. *bow* – Kauai

Have you please considered the interesting question of estimating the period instead of the frequency, and then passing to frequency? – Scarecrow

@AntonShepelev: Not in much detail, but at a glance, that seems backwards — if you're interested in the period (or rather the mean time between events, since they may well not be periodic), I'd recommend estimating the frequency (i.e. mean event rate) and taking its inverse to get an estimate of the period. The reason is that the expected number of events within a (weighted) time interval is a linear function of the frequency, but not of the period, and thus, conversely, a weighted integral of the observed events over time can be used as an estimator of the frequency, like in my answer above. – Culottes

Answer 2

You could try this:

Keep an estimator zn so that at each event:

z_n = (z_n-1+κ).e^{-κ.(t_n-t_n-1)}

This will converge towards the event rate in s^-1. A sligtly better estimator is then (as there is still an error/noise related if you compute the estimate just before or just after an event) :

w_n = z_n.e^-κ/(2.z_n)

In your example it will converge to 2s^-1 (the inverse of 500ms)

The constant κ is responsible for the smoothing and is in s^-1. Small values will smooth more. If your event rate is roughly of seconds, a value of 0.01s-1

for κ is a good start.

This method has a starting bias, and z₀ could be set to an estimate of the value for faster convergence. Small values of κ will keep the bias longer.

There are much more powerful ways of analyzing poisson-like distributions, but they often require large buffers. Frequency analysis such as Fourier transform is one.

Houdan

This is basically the same as the formula I gave, except that you're tracking the average as it was just before the most recent event, whereas I'm tracking the average just after the most recent event. Of course, either way works; they just differ by k. Anyway, +1 for being quicker. :-) – Culottes

At least you prove it, I don't :) (I did it just intuitively and checked on a random serie...) – Houdan

After the last edit mine is more accurate as it compensates for the bias related to the interval – Houdan

It appears that you're calculating something slightly different than me now: my λ* is a simple exponentially weighted average of the recent event rate, and will decay exponentially towards zero if no new events occur, whereas you seem to be projecting forwards from the last event to calculate an "average of the average" (and your estimate only changes when events happen). I suppose I could add some kind of smoothing term to my formula too, if the discontinuities in the average at each event are considered undesirable; I'll have to think about it a bit more... – Culottes

Can you explain the derivation of w sub n, and how you would adjust if estimating after the time of the last sample? – Gunboat