The kantorovich
package has two main features:
With the help of the rccd
and gmp
packages, the kantorovich
package can return the exact values of the extreme joinings and of the Kantorovich distance.
As an example, take $\mu$ and $\nu$ the uniform probability measures on a finite set having three elements.
mu <- nu <- c(1/3, 1/3, 1/3)
The ejoinings
function returns the extreme joinings of $\mu$ and $\nu$. In this case these are the $6!$ permutation matrices:
library(kantorovich) ejoinings(mu, nu)
Since mu
and nu
were unnamed, the vector names c(1,2,3)
has been automatically assigned to them.
The Kantorovich distance between $\mu$ and $\nu$ is relative to a given distance on the state space of $\mu$ and $\nu$, represented by their vector names. By default, the kantorovich
package takes the discrete $0\mathrm{-}1$ distance. Obviously the Kantorovich distance is $0$ on this example, because $\mu=\nu$.
kantorovich(mu, nu)
Note the message returned by both the ejoinings
and the kantorovich
functions. In order to get exact results, use rational numbers with the gmp
package:
library(gmp) mu <- nu <- as.bigq(c(1,1,1), c(3,3,3)) # shorter: as.bigq(c(1,1,1), 3) ejoinings(mu, nu)
Let us try an example with a user-specified distance.
Let's say that the state space of $\mu$ and $\nu$ is ${a, b, c}$,
and then we use c("a","b","c")
as the vector names.
mu <- as.bigq(c(1,2,4), 7) nu <- as.bigq(c(3,1,5), 9) names(mu) <- names(nu) <- c("a", "b", "c")
The distance can be specified as a matrix.
Assume the distance $\rho$ is given by $\rho(a,b)=1$, $\rho(a,c)=2$ and $\rho(b,c)=4$.
The bigq
matrices offered by the gmp
package do not handle dimension names.
But, in our example, the distance $\rho$ takes only integer values,
therefore one can use a numerical matrix:
M <- matrix( c( c(0, 1, 2), c(1, 0, 4), c(2, 4, 0) ), byrow = TRUE, nrow = 3, dimnames = list(c("a","b","c"), c("a","b","c"))) kantorovich(mu, nu, dist=M)
If the distance takes rational values, one can proceed as before with a character matrix:
M <- matrix( c( c("0", "3/13", "2/13"), c("1/13", "0", "4/13"), c("2/13", "4/13", "0") ), byrow = TRUE, nrow = 3, dimnames = list(c("a","b","c"), c("a","b","c"))) kantorovich(mu, nu, dist=M)
One can enter the distance as a function. In such an example, this does not sound convenient:
rho <- function(x,y){ if(x==y) { return(0) } else { if(x=="a" && y=="b") return(1) if(x=="a" && y=="c") return(2) if(x=="b" && y=="c") return(4) return(rho(y,x)) } } kantorovich(mu, nu, dist=rho)
Using a function could be more convenient in the case when the names are numbers:
names(mu) <- names(nu) <- 1:3
But one has to be aware that there are in character mode:
names(mu)
Thus, one can define a distance function as follows, for example with $\rho(x,y)=\frac{|x-y|}{1+|x-y|}$:
rho <- function(x,y){ x <- as.numeric(x); y <- as.numeric(y) return(as.bigq(abs(x-y), 1+abs(x-y))) } kantorovich(mu, nu, dist=rho)
The kantorovich
package also handles the case when mu
and nu
have different lengths, such as this example:
mu <- as.bigq(c(1,2,4), 7) nu <- as.bigq(c(3,1), 4) names(mu) <- c("a", "b", "c") names(nu) <- c("b", "c") ejoinings(mu, nu) kantorovich(mu, nu)
Note the caution message. The kantorovich
package has to handle the fact that mu
is given on the set ${a,b,c}$ while nu
is given on the set ${b,c}$. It detects that the second set is included in the first one, and then treats nu
on the bigger set ${a,b,c}$ by assigning $\nu(a)=0$. To avoid this caution message, the user has to enter this $0$ value:
nu <- as.bigq(c(0,3,1), 4) names(nu) <- c("a", "b", "c")
The kantorovich
package provides three other functions to compute the Kantorovich distance:
kantorovich_lp
, which uses the lp_solve solver with the help of the lpSolve
package;
kantorovich_glpk
, which uses the GLPK solver with the help of the Rglpk
package.
kantorovich_CVX
, which uses the ECOS solver with the help of the CVXR
package.
Contrary to the kantorovich
function, these two functions do not take care of the names of the two vectors mu
and nu
representing the two probability measures $\mu$ and $\nu$, and the distance to be minimized on average must be given as a matrix only, not a function.
A better precision is achieved by kantorovich_glpk
. For instance, take the previous example for which we found $13/63$ as the exact Kantorovich distance:
mu <- c(1,2,4)/7 nu <- c(3,1,5)/9 M <- matrix( c( c(0, 1, 2), c(1, 0, 4), c(2, 4, 0) ), byrow = TRUE, nrow = 3) kanto_lp <- kantorovich_lp(mu, nu, dist=M) kanto_glpk <- kantorovich_glpk(mu, nu, dist=M) kanto_CVX <- kantorovich_CVX(mu, nu, dist=M)
Then kantorovich_lp
and kantorovich_CVX
do not return the better decimal approximation of $13/63$:
print(kanto_lp, digits=22) print(kanto_glpk, digits=22) print(kanto_CVX, digits=22) print(13/63, digits=22)
But kantorovich_CVX
is the fastest one, and it handles the case when the
marginal probability measures mu
and nu
have a large support.
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.