Is there an R or python package to calculate wasserstein metric between negative binomial distributions? As the title says I am looking for an R or python package which can calculate wasserstein distance (aka earthmovers distance between) two lists (vectors) of sampled values from a negative binomial distribution?
Relevant publication: Barbour, A. D., Gan, H. L., & Xia, A. (2015). Stein factors for negative binomial approximation in Wasserstein distance. Bernoulli, 21(2), 1002-1013.
Explanation of the Wasserstien distance: https://rpubs.com/FJRubio/NWD
 A: In R, you can use emd() in the emdist package. In Python, there is scipy.stats.wasserstein_distance. Related: Computing Wasserstein Distance.

I am not entirely sure about emdist, though. Looking at Wikipedia under Examples/One-dimensional distributions, we see a simple formula we should be able to code immediately:
$$ W_1 = \int_{\mathbb{R}} |F_1(x)-F_2(x)|\,dx, $$
where the $F_i$ are the cumulative distribution functions of your two distributions. Since you already have samples, this can just be evaluated using the ecdf() function. However, the emdist::emd() function and this direct evaluation give me two quite different results, see below. So it would be good to check the correctness of either approach with a simpler example, because either Wikipedia, or the emdist package, or I have an error.
set.seed(1)
samples <- list(rnbinom(100,5,0.5), rnbinom(120,4,0.4))
par(mfrow=c(2,1),mai=c(.5,.5,.5,.1))
for ( ii in 1:2 ) {
    plot(as.numeric(names(table(samples[[ii]]))),table(samples[[ii]]),
        main=paste("Sample",ii),las=1,type="h",xlab="",ylab="",lwd=5,col="grey",
        xlim=range(unlist(samples)),
        ylim=range(as.numeric(unlist(sapply(samples,table),names))))
}

wasserstein <- 0
for ( xx in 0:max(unlist(samples)) ) {
    wasserstein <- wasserstein+abs(ecdf(samples[[1]])(xx)-ecdf(samples[[2]])(xx))
}
(wasserstein <- wasserstein/(max(unlist(samples))+1))
# 0.04361111    

library(emdist)
emd(cbind(1,samples[[1]]),cbind(1,samples[[2]]))
# 0.18


