Earth Mover's distanceのR計算

• 2次元の分布としてのメタボロームデータのことをやっている
• EMdist,Wass...についてはこちらを
• そのほかとしてはTransportation problemであって、その解法はハンガリアンメソッドだという
• ハンガリアンメソッドは(重み付き)割り当てマッチング(こちら)でもあって、こちらにもつながる
• EMDのCコードはこちら
• Rにて
  • 行列型データの場合

library(emdist)
n <- 10
m <- 20
A <- matrix(rnorm(n*m),n,m)
A <- A/sum(A)
B <- matrix(rnorm(n*m),n,m)
B <- B/sum(B)
emd2d(A,B)

• • 行列サイズが大きくなると回らない
  • 大きいときは一部を切り取ってみよう。ある値より大きい地点数がそれほど多くないなら…

# 大きな行列を作る
n <- 1000
m <- 2000
A <- matrix(rnorm(n*m),n,m)
B <- matrix(rnorm(n*m),n,m)
# カットオフを決めるために２行列の値のベクトルを作って
tandem.AB <- c(A,B)
# そのtop 100をカットオフとする
s.t.AB <- sort(tandem.AB,decreasing = TRUE)
cutoff <- s.t.AB[100]
# カットオフ以上の値の番地をとり出す
A. <- which(A >= cutoff,arr.ind=TRUE)
B. <- which(B >= cutoff,arr.ind=TRUE)
# 選ばれた地点の値を取り出す
A.w <- A[A.]
B.w <- B[B.]
# 第１列を値、第２，３列を番地にした行列を作る
A.2 <- cbind(A.w,A.)
B.2 <- cbind(B.w,B.)
# 併せて１にする
A.2[,1] <- A.2[,1]/sum(A.2[,1])
B.2[,1] <- B.2[,1]/sum(B.2[,1])
# emdを計算する
emd(A.2,B.2)

• • 2軸の扱いに重みをつけるなら

k1 <- 1
k2 <- 10
A.3 <- A.2
B.3 <- B.2
A.3[,2] <- A.3[,2]*k1
A.3[,3] <- A.3[,3]*k2
B.3[,2] <- B.3[,2]*k1
B.3[,2] <- B.3[,3]*k3
emd(A.3,B.3)
# とか、下駄をはかせて
A.4 <- A.2
B.4 <- B.2
A.4[,2] <- A.4[,2]+k1
A.4[,3] <- A.4[,3]+k2
B.4[,2] <- B.4[,2]+k1
B.4[,2] <- B.4[,3]+k3
emd(A.4,B.4)