kantorovich packageThe 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.
The ejoinings function returns the extreme joinings of
\(\mu\) and \(\nu\). In this case these are the \(6!\) permutation matrices:
## Message: You should enter mu and nu in rational with the gmp package.
## [[1]]
## 1 2 3
## 1 0.3333333 0.0000000 0.0000000
## 2 0.0000000 0.0000000 0.3333333
## 3 0.0000000 0.3333333 0.0000000
##
## [[2]]
## 1 2 3
## 1 0.3333333 0.0000000 0.0000000
## 2 0.0000000 0.3333333 0.0000000
## 3 0.0000000 0.0000000 0.3333333
##
## [[3]]
## 1 2 3
## 1 0.0000000 0.3333333 0.0000000
## 2 0.0000000 0.0000000 0.3333333
## 3 0.3333333 0.0000000 0.0000000
##
## [[4]]
## 1 2 3
## 1 0.0000000 0.3333333 0.0000000
## 2 0.3333333 0.0000000 0.0000000
## 3 0.0000000 0.0000000 0.3333333
##
## [[5]]
## 1 2 3
## 1 0.0000000 0.0000000 0.3333333
## 2 0.0000000 0.3333333 0.0000000
## 3 0.3333333 0.0000000 0.0000000
##
## [[6]]
## 1 2 3
## 1 0.0000000 0.0000000 0.3333333
## 2 0.3333333 0.0000000 0.0000000
## 3 0.0000000 0.3333333 0.0000000
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\).
## Message: You should enter mu and nu in rational with the gmp package.
## [1] 0
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)## [[1]]
## 1 2 3
## 1 "1/3" "0" "0"
## 2 "0" "0" "1/3"
## 3 "0" "1/3" "0"
##
## [[2]]
## 1 2 3
## 1 "1/3" "0" "0"
## 2 "0" "1/3" "0"
## 3 "0" "0" "1/3"
##
## [[3]]
## 1 2 3
## 1 "0" "1/3" "0"
## 2 "0" "0" "1/3"
## 3 "1/3" "0" "0"
##
## [[4]]
## 1 2 3
## 1 "0" "1/3" "0"
## 2 "1/3" "0" "0"
## 3 "0" "0" "1/3"
##
## [[5]]
## 1 2 3
## 1 "0" "0" "1/3"
## 2 "0" "1/3" "0"
## 3 "1/3" "0" "0"
##
## [[6]]
## 1 2 3
## 1 "0" "0" "1/3"
## 2 "1/3" "0" "0"
## 3 "0" "1/3" "0"
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.
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)## Big Rational ('bigq') :
## [1] 13/63
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)## Big Rational ('bigq') :
## [1] 1/63
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)## Big Rational ('bigq') :
## [1] 13/63
Using a function could be more convenient in the case when the names are numbers:
But one has to be aware that there are in character mode:
## [1] "1" "2" "3"
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)## Big Rational ('bigq') :
## [1] 37/378
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)## Caution: some names of mu and/or nu were missing or not compatible - automatic change
## [[1]]
## b c
## a "1/7" "0"
## b "1/28" "1/4"
## c "4/7" "0"
##
## [[2]]
## b c
## a "1/7" "0"
## b "2/7" "0"
## c "9/28" "1/4"
##
## [[3]]
## b c
## a "0" "1/7"
## b "5/28" "3/28"
## c "4/7" "0"
##
## [[4]]
## b c
## a "0" "1/7"
## b "2/7" "0"
## c "13/28" "3/28"
## Caution: some names of mu and/or nu were missing or not compatible - automatic change
## Big Rational ('bigq') :
## [1] 13/28
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:
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)
## [1] 0.2063492063492062544849
print(kanto_glpk, digits=22)
## [1] 0.2063492063492063377517
print(kanto_CVX, digits=22)
## [1] 0.2063492063214846794494
print(13/63, digits=22)
## [1] 0.2063492063492063377517But kantorovich_CVX is the fastest one, and it handles
the case when the marginal probability measures mu and
nu have a large support.