# Earth mover's distance implementation for circular distributions?

I am interested in finding an implementation (preferably in R but not necessary) to calculate earth mover's distance between two empirical distributions of circular data. My data are time points on a 24-hour clock under different treatments. I am not too familiar with circular statistics, so it would also be very helpful for me to get a brief explanation of how the calculation of EMD differs for circular distributions versus scale data.

# sample data

These two vectors represent the times of day, rounded to the nearest hour, that an individual ant from a particular species was observed under two different temperature treatments.

x1 <- c(23, 3, 0, 8, 3, 10, 2, 23, 0, 23, 10, 12, 9, 0, 3, 23, 23,
10, 9, 23, 3, 23, 3, 2, 21, 0, 0, 2, 16, 23, 3, 21, 22, 21, 10,
23, 23, 23, 2, 3, 23, 3, 10, 10, 3, 23, 10, 2, 23, 10, 15, 3,
0, 10, 10, 13, 10, 18, 2, 8, 23, 2, 10, 0, 23, 23, 23, 2, 21,
16, 0, 23, 3, 10, 0, 3, 10, 23, 23, 3, 3, 10, 21, 0, 3, 15, 19,
0, 10, 21, 2, 21, 8, 23, 8, 9, 21, 3, 3, 3)

x2 <- c(20, 6, 10, 7, 7, 16, 20, 12, 6, 18, 15, 10, 10, 6, 10, 21,
10, 5, 12, 20, 12, 6, 0, 19, 0, 16, 20, 0, 5, 16, 6, 20, 12,
15, 7, 6, 16, 10, 12, 7, 5, 19, 16, 20, 5, 6, 15, 10, 16, 16,
20, 10, 20, 7, 7, 10, 20, 20, 10, 19, 0, 10, 10, 19, 12, 6, 19,
20, 0, 10, 12, 0, 7, 6, 10, 15, 16, 6, 6, 19, 10, 21, 20, 0,
15, 6, 10, 20, 7, 6, 10, 0, 16, 19, 6, 20, 10, 7, 12, 7)


Algorithms for computing EMD is discussed at Coupling and Total variational distance and Calculate Earth Mover's Distance for two grayscale images, and also over at math SE. There are some R packages, some fast examples (with your data):

x1.circ <- circular::circular(x1, units="hours", template="clock24")
x2.circ <- circular::circular(x2, units="hours", template="clock24")


par(mfrow=c(1, 2))
circular::plot.circular(x1.circ, stack=TRUE)
circular::plot.circular(x2.circ, stack=TRUE)


Then for the EMD distance using R package emdist:

calc_weight <- function(x) { # a vector of hours
tab <- table(factor(x,  levels=as.character(0:23)),
useNA="ifany")

dimnames(tab) <- NULL
mat <- cbind( weights=tab/sum(tab), points=0:23 )
mat
}

A <- calc_weight(x1)
B <- calc_weight(x2)

hourdist <- function(A, B) sum(pmin(  (A-B)%%24, (B-A)%%24 ) )

emdist::emd(A, B, dist=hourdist)
[1] 2.63


The package transport can give more information:

costm <- outer(0:23, 0:23, FUN=function(x, y) pmin( (y-x)%%24, (x-y)%%24 ) )
transport::transport(A[, 1], B[, 1], costm)

from to mass
1     1  1 0.08
2     1  4 0.02
3     1 12 0.01
4     2  2 0.04
5     2  4 0.05
6     3  3 0.14
7     3  4 0.03
8     3  6 0.01
9     4  5 0.04
10    5  6 0.03
11    6  5 0.14
12    6  6 0.02
13    7  6 0.01
14    8  6 0.01
15    9  7 0.02
16   10  8 0.02
17   11  8 0.01
18   12  9 0.01
19   13  7 0.02
20   13 11 0.06
21   14  7 0.01
22   15  8 0.06
23   15 10 0.07
24   15 11 0.08
25   15 12 0.01