# Average and standard deviation of timestamps (time wraps around at midnight)

I have lots of sensor data with timestamps like "2014-09-09 16:10:45" and accompanying sensor readings. To get some insight into these I want to find "unusual" events by looking at the average and standard deviation of the time part of the timestamp. How can I handle the wrap-around of time on midnight?

A made up example: Imagine power readings being influenced by people turning on machines in the morning (a power sensor would notice increasing values) and turning them off in the evening (decreasing values). I want to find sensor readings that are unusual. This would be decreasing sensor readings in time periods where readings usually increase and increasing readings in time periods where readings usually decrease.

My idea was to extract the time part of the timestamps (eg. 12:55:10), convert that to seconds (the day has 86400) and then, split by the tendency of readings (eg only looking at increasing readings) calculate the average and standard deviation. If I then take the time window from "average second of the day minus standard deviation" to "average second of the day plus standard deviation" (maybe using twice the standard deviation) I would have the typical periods and every increasing reading outside this time window would be unusual.

The problem: Time wraps around at midnight! A reading at 00:15:00 would actually be really close to 23:50:00 in reality, but "far away" in calculation. This surely skews the statistics unless everything happens mid-day. Is there a standard practice to handle this? Can you give me ideas? I am totally stumped at the moment. I would love to stay in PostgresQL but as that is not a requirement I did not tag that. Anything helps!

Below is some example data, I have about 200-300 readings per sensor. You can see that in this example the increases happen in the morning.

"Timestamp as %Y-%m-%d %H:%M:%S";"Day of the year";"Second of the day";"Tendency of reading"
"2014-03-01 14:45:00";60;53100;-0.030
"2014-03-03 08:18:00";62;29880;0.150
"2014-03-03 14:17:00";62;51420;-0.120
"2014-03-03 16:37:00";62;59820;-0.030
"2014-03-04 08:11:00";63;29460;0.150
"2014-03-04 10:21:00";63;37260;-0.150
"2014-03-04 16:12:00";63;58320;-0.030
"2014-03-05 08:04:00";64;29040;0.150
"2014-03-05 14:42:00";64;52920;-0.060
"2014-03-05 17:27:00";64;62820;-0.030
"2014-03-06 08:29:00";65;30540;0.090
"2014-03-06 12:06:00";65;43560;-0.030
"2014-03-06 13:49:00";65;49740;-0.120
"2014-03-07 08:21:00";66;30060;0.150
"2014-03-07 10:27:00";66;37620;-0.030
"2014-03-07 11:27:00";66;41220;0.030
"2014-03-07 13:46:00";66;49560;-0.060
"2014-03-07 16:59:00";66;61140;-0.030
"2014-03-07 18:52:00";66;67920;-0.030
"2014-03-08 08:47:00";67;31620;0.120

• This seems like a duplicate of stats.stackexchange.com/questions/109704 (currently unanswered, alas), but perhaps the conversation in the comments there might help get you started? – whuber Sep 9 '14 at 15:03
• Thank you so very much, I was not aware of the term circular statistics so that was a great pointer. I'll see what I manage to scrap together. – inc42 Sep 9 '14 at 16:35

Let's use the simplification you suggest: only use the data from positive readings and disregard the value of the reading, so we are left with a single set of circular data. You can use the circular dispersion, as whuber suggested, possibly multiplied by some constant to determine how much of the data should be seen as an outlier. A good text that is slightly easier to understand than the Wikipedia page would be Statistical Analysis of Circular data, by N.I. Fisher (1995).

I'll give some more straightforward formulas than the Wiki page, and give some sample code.

The dispersion can be calculated as (due to Fisher (p.32-34)):

1. Denote the data by $\boldsymbol\theta = \{\theta_1, \dots, \theta_n\}.$ An estimate of the mean $\bar\theta$ can be calculated with

2. Calculate $\bar{R} = \frac{\sqrt{S^2 + C^2}}{n}$.

3. Calculate the dispersion as suggested by whuber. I'm not sure why, but Wikipedia's definition seems to differ slightly from Fisher. I will use Fisher's:

• $\hat\delta = \frac{1 - \left[ (1/n) \sum_{i=1}^{n} \cos 2 (\theta_i - \hat\mu) \right]}{2\bar{R}^2}.$
4. Then, choose some constant $c$ (1 is probably fine, but you may fine-tune). Then, the interval is given by

• $\left[\hat\mu - c \hat\delta, \hat\mu + c \hat\delta \right]$.

I know you want to avoid R, but just to show how to calculate this in code, here is some basic R code anyway, which also generates a plot:

n  <- 200
th <- runif(n, 0.5 * pi, 1.5 * pi)

plot(cos(th), sin(th), xlim=c(-1, 1), ylim = c(-1, 1))

S <- sum(sin(th))
C <- sum(cos(th))

mu_hat <- atan2(S, C)

R_bar  <- sqrt(S^2 + C^2) / n

delta_hat <- (1 - sum(cos(2 * (th-mu_hat)))/n) / (2 * R_bar^2)

constant  <- 0.8

CI <- mu_hat + c(-1, 1) * constant * delta_hat

lines(x = c(0, cos(CI[1])),
y = c(0, sin(CI[1])), col="green")

lines(x = c(0, cos(CI[2])),
y = c(0, sin(CI[2])), col="blue")


In a final note, it may still be better to use the additional information provided by the value of the reading, and not only the sign, because it may provide better estimates. However, the simplification of only using the sign makes the problem much more manageable. If anyone has a good solution that incorporates the readings, I would love to know!

• This is a good answer (+1). When defining the mean, I believe your intention is to average, rather than simply sum, the values of $\sin(\theta)$ and $\cos(\theta)$. – Sycorax says Reinstate Monica Sep 11 '14 at 14:28
• Thanks. In fact, $\bar{S} = S/n$ and $\bar{C} = C/n$ are commonly useful, as they provide the coordinates to the centroid as you calculated it. However, here, it is irrelevant whether we sum or average, because $\text{atan2}(x, y) = \text{atan2}(cx, cy)$, for each constant $c \neq 0$.} – Kees Mulder Sep 11 '14 at 14:44
• This is fantastic and that you wrote it in easier steps that Wikipedia helps me a lot. Thank you very much! – inc42 Sep 11 '14 at 18:07
• I struggle with the magical constant. Unfortunately pages 33 and 34 are excluded on Google Books (weird!) so I cannot read up on it myself. When testing with the actual data a constant of 1 was often way too narrow. See i.imgur.com/RY0Ousa.png for example for the data pastebin.com/Jn7f5Tq1 . In this case I would need a constant of about 5-6 to get a useful "those other points are at the edges" selection. What exactly does the constant mean statistically? – inc42 Oct 1 '14 at 20:00
• The dispersion should be somewhat comparable to a standard deviation. With the standard deviation, we can state, for example, that in a normal distribution only 2.5% of the data lies more than 1.96 standard deviations above the mean. Here, we have not assumed any model (except symmetry), but if we would assume any distribution, the calculated interval here would correspond to some percentage of data we would expect outside this interval. Usually this percentage would be low, so that the outliers are really extreme. So, set $c$ such that most non-outlier data falls within the interval. – Kees Mulder Oct 2 '14 at 11:58

This might be kind of dumb, but I just plotted the points as if they lie in a disc: the angle is taken to be time, and I assumed unit radius. The angle of the centroid of this point cloud is the mean angle of your data, i.e. the mean time of day which respects the "wrap around at midnight" property. This value also has the nice property of maximizing the MLE for the location parameter of a von Mises distribution.

Since all of our points lie roughly between 5 a.m. and 5 p.m, it shouldn't be surprising that their mean is near noon.

But, at least for the moment, this is all of the circular statistics that I understand! I wish I could give you some more help, but I'm still puzzling through the wikipedia articles.

sec.radians <- 2*pi*sec/(60*60*24)

theta   <- seq(0, 2*pi, by=0.01)
lines(cos(theta), sin(theta), col="red", lty="dashed")

landmarks   <- c(2*pi, 3*pi/2, pi, pi/2)
text(0.5*cos(landmarks), 0.5*sin(landmarks), c("midnight", "6 a.m.", "noon", "6 p.m."))

theta.mean  <- atan2(centroid$y, centroid$x)

• @inc42 I've revised it to be a little more coherent. I'm playing around with the R package circular, and I'll let you know if I manage to do anything helpful. – Sycorax says Reinstate Monica Sep 9 '14 at 17:23