Earth mover's distance implementation for circular distributions?

Question

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)

Answer 1

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")

plots made by

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