c(x,y) = dxy^p; ground cost/metric
1 2 3 4 5 6 7 |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 | ## 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)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.