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, 2014 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. Sep 9, 2014 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
Sep 11, 2014 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$.} Sep 11, 2014 at 14:44
• This is fantastic and that you wrote it in easier steps that Wikipedia helps me a lot. Thank you very much! Sep 11, 2014 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? Oct 1, 2014 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. Oct 2, 2014 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."))

points(centroid, col="purple", lwd=5)
theta.mean  <- atan2(centroid$y, centroid$x)

• If you were to ignore the intensity and plot the points on the circle, then computed the centroid, you would reconstruct the basic location summary used in circular statistics. Scaling them by intensity seems of questionable value for the purpose of finding outliers. See the rest of that linked article for information on circular dispersion, which could be used to detect outliers.
– whuber
Sep 9, 2014 at 16:23
• @whuber I am not a smart man.
– Sycorax
Sep 9, 2014 at 16:36
• Thanks, this is very helpful! The readings would be irrevelant, they are only used to divide the timestamps into "this timestamp is of an increasing reading" and "this one is of a decrease". All points will have the same radius (only interested in the sec) and only the angle is used for statistics. I think I will be able to work from your example and whuber's link seems perfect (although I have lots of trouble with mathematical terms and writing, we'll see). Sep 9, 2014 at 17:13
• @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
Sep 9, 2014 at 17:23
• Re "smart", I would beg to differ with that negative assessment: you have independently recreated a legitimate and effective statistical method to characterize circular data. (+1)
– whuber
Sep 9, 2014 at 17:33