coupling: optimal transport; returns p-Wasserstein distance

In kisungyou/DAS: Distance-bAsed Statistical methods

Description

c(x,y) = dxy^p; ground cost/metric

Usage

coupling(
  dxy,
  p = 1,
  wx,
  wy,
  method = c("networkflow", "shortsimplex", "revsimplex", "primaldual")
)

Examples

## create two small datasets from bivariate normal
X = matrix(rnorm(5*2),ncol=2) # 5 obs. for X
Y = matrix(rnorm(5*2),ncol=2) # 5 obs. for Y

## compute cross-distance between X and Y
dXY = array(0,c(5,5))
for (i in 1:5){
  vx = as.vector(X[i,])
  for (j in 1:5){
    vy  = as.vector(Y[j,])
    dXY[i,j] = sqrt(sum((vx-vy)^2))
  }
}

## compute the distance and report
output = coupling(dXY, p=2) # 2-Wasserstein distance
image(output$coupling, main=paste("distance=",round(output$distance,4),sep=""))

## Not run: 
## create two datasets from bivariate normal
## let's try to see the evolution of 2-Wasserstein distance
nmax = 1000
X = matrix(rnorm(nmax*2),ncol=2) # obs. for X
Y = matrix(rnorm(nmax*2),ncol=2) # obs. for Y

## compute cross-distance between X and Y
dXY = array(0,c(nmax,nmax))
for (i in 1:nmax){
  vx = as.vector(X[i,])
  for (j in 1:nmax){
    vy  = as.vector(Y[j,])
    dXY[i,j] = sqrt(sum((vx-vy)^2))
  }
}

## compute
xgrid = 2:nmax
ygrid = rep(0,nmax-1)
for (i in 1:(nmax-1)){
  pXY = dXY[1:(i+1),1:(i+1)]
  ygrid[i] = coupling(pXY, p=2)$distance
  print(paste("Iteration ",i+1,"/",nmax," Complete..",sep=""))
}

## visualize
plot(xgrid, ygrid, "b", lwd=1, main="Evolution of 2-Wasserstein Distances",
     xlab="number of samples", ylab="distance", pch=18)

## End(Not run)

kisungyou/DAS documentation built on Jan. 6, 2020, 7:09 a.m.