首页 > 代码库 > 概率分布之间的距离度量以及python实现
概率分布之间的距离度量以及python实现
1. 欧氏距离(Euclidean Distance)
欧氏距离是最易于理解的一种距离计算方法,源自欧氏空间中两点间的距离公式。
(1)二维平面上两点a(x1,y1)与b(x2,y2)间的欧氏距离:
(2)三维空间两点a(x1,y1,z1)与b(x2,y2,z2)间的欧氏距离:
(3)两个n维向量a(x11,x12,…,x1n)与 b(x21,x22,…,x2n)间的欧氏距离:
(4)也可以用表示成向量运算的形式:
python中的实现:
方法一:
import numpy as np x=np.random.random(10) y=np.random.random(10) #方法一:根据公式求解 d1=np.sqrt(np.sum(np.square(x-y))) #方法二:根据scipy库求解 from scipy.spatial.distance import pdist X=np.vstack([x,y]) d2=pdist(X)
2. 曼哈顿距离(Manhattan Distance)
从名字就可以猜出这种距离的计算方法了。想象你在曼哈顿要从一个十字路口开车到另外一个十字路口,驾驶距离是两点间的直线距离吗?显然不是,除非你能穿越大楼。实际驾驶距离就是这个“曼哈顿距离”。而这也是曼哈顿距离名称的来源, 曼哈顿距离也称为城市街区距离(City Block distance)。
(1)二维平面两点a(x1,y1)与b(x2,y2)间的曼哈顿距离
(2)两个n维向量a(x11,x12,…,x1n)与 b(x21,x22,…,x2n)间的曼哈顿距离
python中的实现 :
import numpy as np x=np.random.random(10) y=np.random.random(10) #方法一:根据公式求解 d1=np.sum(np.abs(x-y)) #方法二:根据scipy库求解 from scipy.spatial.distance import pdist X=np.vstack([x,y]) d2=pdist(X,‘cityblock‘)
3. 切比雪夫距离 ( Chebyshev Distance )
国际象棋玩过么?国王走一步能够移动到相邻的8个方格中的任意一个。那么国王从格子(x1,y1)走到格子(x2,y2)最少需要多少步?自己走走试试。你会发现最少步数总是max( | x2-x1 | , | y2-y1 | ) 步 。有一种类似的一种距离度量方法叫切比雪夫距离。
(1)二维平面两点a(x1,y1)与b(x2,y2)间的切比雪夫距离
(2)两个n维向量a(x11,x12,…,x1n)与 b(x21,x22,…,x2n)间的切比雪夫距离
这个公式的另一种等价形式是
看不出两个公式是等价的?提示一下:试试用放缩法和夹逼法则来证明。
在python中的实现:
import numpy as np x=np.random.random(10) y=np.random.random(10) #方法一:根据公式求解 d1=np.max(np.abs(x-y)) #方法二:根据scipy库求解 from scipy.spatial.distance import pdist X=np.vstack([x,y]) d2=pdist(X,‘chebyshev‘)
4. 闵可夫斯基距离(Minkowski Distance)
闵氏距离不是一种距离,而是一组距离的定义。
(1) 闵氏距离的定义
两个n维变量a(x11,x12,…,x1n)与 b(x21,x22,…,x2n)间的闵可夫斯基距离定义为:
也可写成
其中p是一个变参数。
当p=1时,就是曼哈顿距离
当p=2时,就是欧氏距离
当p→∞时,就是切比雪夫距离
根据变参数的不同,闵氏距离可以表示一类的距离。
(2)闵氏距离的缺点
闵氏距离,包括曼哈顿距离、欧氏距离和切比雪夫距离都存在明显的缺点。
举个例子:二维样本(身高,体重),其中身高范围是150~190,体重范围是50~60,有三个样本:a(180,50),b(190,50),c(180,60)。那么a与b之间的闵氏距离(无论是曼哈顿距离、欧氏距离或切比雪夫距离)等于a与c之间的闵氏距离,但是身高的10cm真的等价于体重的10kg么?因此用闵氏距离来衡量这些样本间的相似度很有问题。
简单说来,闵氏距离的缺点主要有两个:(1)将各个分量的量纲(scale),也就是“单位”当作相同的看待了。(2)没有考虑各个分量的分布(期望,方差等)可能是不同的。
python中的实现:
import numpy as np x=np.random.random(10) y=np.random.random(10) #方法一:根据公式求解,p=2 d1=np.sqrt(np.sum(np.square(x-y))) #方法二:根据scipy库求解 from scipy.spatial.distance import pdist X=np.vstack([x,y]) d2=pdist(X,‘minkowski‘,p=2)
5. 标准化欧氏距离 (Standardized Euclidean distance )
(1)标准欧氏距离的定义
标准化欧氏距离是针对简单欧氏距离的缺点而作的一种改进方案。标准欧氏距离的思路:既然数据各维分量的分布不一样,好吧!那我先将各个分量都“标准化”到均值、方差相等吧。均值和方差标准化到多少呢?这里先复习点统计学知识吧,假设样本集X的均值(mean)为m,标准差(standard deviation)为s,那么X的“标准化变量”表示为:
标准化后的值 = ( 标准化前的值 - 分量的均值 ) /分量的标准差
经过简单的推导就可以得到两个n维向量a(x11,x12,…,x1n)与 b(x21,x22,…,x2n)间的标准化欧氏距离的公式:
如果将方差的倒数看成是一个权重,这个公式可以看成是一种加权欧氏距离(Weighted Euclidean distance)。
python中的实现:
import numpy as np x=np.random.random(10) y=np.random.random(10) X=np.vstack([x,y]) #方法一:根据公式求解 sk=np.var(X,axis=0,ddof=1) d1=np.sqrt(((x - y) ** 2 /sk).sum()) #方法二:根据scipy库求解 from scipy.spatial.distance import pdist d2=pdist(X,‘seuclidean‘)
6. 马氏距离(Mahalanobis Distance)
(1)马氏距离定义
有M个样本向量X1~Xm,协方差矩阵记为S,均值记为向量μ,则其中样本向量X到u的马氏距离表示为:
而其中向量Xi与Xj之间的马氏距离定义为:
若协方差矩阵是单位矩阵(各个样本向量之间独立同分布),则公式就成了:
也就是欧氏距离了。
若协方差矩阵是对角矩阵,公式变成了标准化欧氏距离。
python 中的实现:
import numpy as np x=np.random.random(10) y=np.random.random(10) #马氏距离要求样本数要大于维数,否则无法求协方差矩阵 #此处进行转置,表示10个样本,每个样本2维 X=np.vstack([x,y]) XT=X.T #方法一:根据公式求解 S=np.cov(X) #两个维度之间协方差矩阵 SI = np.linalg.inv(S) #协方差矩阵的逆矩阵 #马氏距离计算两个样本之间的距离,此处共有10个样本,两两组合,共有45个距离。 n=XT.shape[0] d1=[] for i in range(0,n): for j in range(i+1,n): delta=XT[i]-XT[j] d=np.sqrt(np.dot(np.dot(delta,SI),delta.T)) d1.append(d) #方法二:根据scipy库求解 from scipy.spatial.distance import pdist d2=pdist(XT,‘mahalanobis‘)
马氏优缺点:
1)马氏距离的计算是建立在总体样本的基础上的,这一点可以从上述协方差矩阵的解释中可以得出,也就是说,如果拿同样的两个样本,放入两个不同的总体中,最后计算得出的两个样本间的马氏距离通常是不相同的,除非这两个总体的协方差矩阵碰巧相同;
2)在计算马氏距离过程中,要求总体样本数大于样本的维数,否则得到的总体样本协方差矩阵逆矩阵不存在,这种情况下,用欧式距离计算即可。
3)还有一种情况,满足了条件总体样本数大于样本的维数,但是协方差矩阵的逆矩阵仍然不存在,比如三个样本点(3,4),(5,6)和(7,8),这种情况是因为这三个样本在其所处的二维空间平面内共线。这种情况下,也采用欧式距离计算。
4)在实际应用中“总体样本数大于样本的维数”这个条件是很容易满足的,而所有样本点出现3)中所描述的情况是很少出现的,所以在绝大多数情况下,马氏距离是可以顺利计算的,但是马氏距离的计算是不稳定的,不稳定的来源是协方差矩阵,这也是马氏距离与欧式距离的最大差异之处。
优点:它不受量纲的影响,两点之间的马氏距离与原始数据的测量单位无关;由标准化数据和中心化数据(即原始数据与均值之差)计算出的二点之间的马氏距离相同。马氏距离还可以排除变量之间的相关性的干扰。缺点:它的缺点是夸大了变化微小的变量的作用。
7. 夹角余弦(Cosine)
也可以叫余弦相似度。 几何中夹角余弦可用来衡量两个向量方向的差异,
import numpy as np x=np.random.random(10) y=np.random.random(10) #方法一:根据公式求解 d1=np.dot(x,y)/(np.linalg.norm(x)*np.linalg.norm(y)) #方法二:根据scipy库求解 from scipy.spatial.distance import pdist X=np.vstack([x,y]) d2=1-pdist(X,‘cosine‘)
两个向量完全相等时,余弦值为1,如下的代码计算出来的d=1。
d=1-pdist([x,x],‘cosine‘)
8. 皮尔逊相关系数(Pearson correlation)
(1) 皮尔逊相关系数的定义
前面提到的余弦相似度只与向量方向有关,但它会受到向量的平移影响,在夹角余弦公式中如果将 x 平移到 x+1, 余弦值就会改变。怎样才能实现平移不变性?这就要用到皮尔逊相关系数(Pearson correlation),有时候也直接叫相关系数。
如果将夹角余弦公式写成:
表示向量x和向量y之间的夹角余弦,则皮尔逊相关系数则可表示为:
皮尔逊相关系数具有平移不变性和尺度不变性,计算出了两个向量(维度)的相关性。
在python中的实现:
import numpy as np x=np.random.random(10) y=np.random.random(10) #方法一:根据公式求解 x_=x-np.mean(x) y_=y-np.mean(y) d1=np.dot(x_,y_)/(np.linalg.norm(x_)*np.linalg.norm(y_)) #方法二:根据numpy库求解 X=np.vstack([x,y]) d2=np.corrcoef(X)[0][1]
相关系数是衡量随机变量X与Y相关程度的一种方法,相关系数的取值范围是[-1,1]。相关系数的绝对值越大,则表明X与Y相关度越高。当X与Y线性相关时,相关系数取值为1(正线性相关)或-1(负线性相关)。
9. 汉明距离(Hamming distance)
(1)汉明距离的定义
两个等长字符串s1与s2之间的汉明距离定义为将其中一个变为另外一个所需要作的最小替换次数。例如字符串“1111”与“1001”之间的汉明距离为2。
应用:信息编码(为了增强容错性,应使得编码间的最小汉明距离尽可能大)。
在python中的实现:
import numpy as np from scipy.spatial.distance import pdist x=np.random.random(10)>0.5 y=np.random.random(10)>0.5 x=np.asarray(x,np.int32) y=np.asarray(y,np.int32) #方法一:根据公式求解 d1=np.mean(x!=y) #方法二:根据scipy库求解 X=np.vstack([x,y]) d2=pdist(X,‘hamming‘)
10. 杰卡德相似系数(Jaccard similarity coefficient)
(1) 杰卡德相似系数
两个集合A和B的交集元素在A,B的并集中所占的比例,称为两个集合的杰卡德相似系数,用符号J(A,B)表示。
杰卡德相似系数是衡量两个集合的相似度一种指标。
(2) 杰卡德距离
与杰卡德相似系数相反的概念是杰卡德距离(Jaccard distance)。杰卡德距离可用如下公式表示:
杰卡德距离用两个集合中不同元素占所有元素的比例来衡量两个集合的区分度。
(3) 杰卡德相似系数与杰卡德距离的应用
可将杰卡德相似系数用在衡量样本的相似度上。
样本A与样本B是两个n维向量,而且所有维度的取值都是0或1。例如:A(0111)和B(1011)。我们将样本看成是一个集合,1表示集合包含该元素,0表示集合不包含该元素。
在python中的实现:
import numpy as np from scipy.spatial.distance import pdist x=np.random.random(10)>0.5 y=np.random.random(10)>0.5 x=np.asarray(x,np.int32) y=np.asarray(y,np.int32) #方法一:根据公式求解 up=np.double(np.bitwise_and((x != y),np.bitwise_or(x != 0, y != 0)).sum()) down=np.double(np.bitwise_or(x != 0, y != 0).sum()) d1=(up/down) #方法二:根据scipy库求解 X=np.vstack([x,y]) d2=pdist(X,‘jaccard‘)
11. 布雷柯蒂斯距离(Bray Curtis Distance)
Bray Curtis距离主要用于生态学和环境科学,计算坐标之间的距离。该距离取值在[0,1]之间。它也可以用来计算样本之间的差异。
样本数据:
计算:
在python中的实现:
import numpy as np from scipy.spatial.distance import pdist x=np.array([11,0,7,8,0]) y=np.array([24,37,5,18,1]) #方法一:根据公式求解 up=np.sum(np.abs(y-x)) down=np.sum(x)+np.sum(y) d1=(up/down) #方法二:根据scipy库求解 X=np.vstack([x,y]) d2=pdist(X,‘braycurtis‘)
1、卡方检验
统计学上的χ2统计量,由于它最初是由英国统计学家Karl Pearson在1900年首次提出的,因此也称之为Pearson χ2,其计算公式为
其中,Ai为i水平的观察频数,Ei为i水平的期望频数,n为总频数,pi为i水平的期望频率。i水平的期望频数Ei等于总频数n×i水平的期望概率pi。当n比较大时,χ2统计量近似服从k-1(计算Ei时用到的参数个数)个自由度的卡方分布。
卡方检验经常用来检验某一种观测分布是不是符合某一类典型的理论分布(如二项分布,正态分布等)。观察频数与期望频数越接近,两者之间的差异越小,χ2值越小;如果两个分布完全一致,χ2值为0;反之,观察频数与期望频数差别越大,两者之间的差异越大,χ2值越大。换言之,大的χ2值表明观察频数远离期望频数,即表明远离假设。小的χ2值表明观察频数接近期望频数,接近假设。因此,χ2是观察频数与期望频数之间距离的一种度量指标,也是假设成立与否的度量指标。如果χ2值“小”,研究者就倾向于不拒绝H0;如果χ2值大,就倾向于拒绝H0。至于χ2在每个具体研究中究竟要大到什么程度才能拒绝H0,则要借助于卡方分布求出所对应的P值来确定(通常取p=0.05)。
在python中的实现:
# -*- coding: utf-8 -*- ‘‘‘ 卡方公式(o-e)^2 / e 期望值和收集到数据不能低于5,o(observed)观察到的数据,e(expected)表示期望的数据 (o-e)平方,最后除以期望的数据e ‘‘‘ import numpy as np from scipy.stats import chisquare list_observe=np.array([30,14,34,45,57,20]) list_expect=np.array([20,20,30,40,60,30]) #方法一:根据公式求解(最后根据c1的值去查表判断) c1=np.sum(np.square(list_observe-list_expect)/list_expect) #方法二:使用scipy库来求解 c2,p=chisquare(f_obs=list_observe, f_exp=list_expect) ‘‘‘ 返回NAN,无穷小 ‘‘‘ if p>0.05 or p=="nan": print("H0 win,there is no difference") else: print("H1 win,there is difference")
2、交叉熵
通常,一个信源发送出什么符号是不确定的,衡量它的不确定性可以根据其出现的概率来度量。概率大,出现机会多,不确定性小;反之就大。
不确定性函数f必须满足两个条件:
1)是概率P的单调递降函数;
2)两个独立符号所产生的不确定性应等于各自不确定性之和,即f(P1,P2)=f(P1)+f(P2),这称为可加性。
同时满足这两个条件的函数f是对数函数,即
在信源中,考虑的不是某一单个符号发生的不确定性,而是要考虑这个信源所有可能发生情况的平均不确定性。若信源符号有n种取值:U1…Ui…Un,对应概率为:P1…Pi…Pn,且各种符号的出现彼此独立。这时,信源的平均不确定性应当为单个符号不确定性-logPi的统计平均值(E),可称为信息熵,即
我们称H(p)为信息熵,称H(p,q)为交叉熵。
交叉熵在CNN分类中经常用到,用来作为预测值和真实标签值的距离度量。经过卷积操作后,最后一层出来的特征经过softmax函数后会变成一个概率向量,我们可以看作为是概率分布q, 而真实标签我们可以看作是概率分布p, 因此真实分布p和预测分布q的交叉熵就是我们要求的loss损失值,即
在python中的实现:
import numpy as np import tensorflow as tf fea=np.asarray([6.5,4.2,7.4,3.5],np.float32) label=np.array([1,0,0,0]) #方法一:根据公式求解 def softmax(x): return np.exp(x)/np.sum(np.exp(x),axis=0) loss1=-np.sum(label*np.log(softmax(fea))) #方法二:调用tensorflow深度学习框架求解 sess=tf.Session() logits=tf.Variable(fea) labels=tf.Variable(label) sess.run(tf.global_variables_initializer()) loss2=sess.run(tf.losses.softmax_cross_entropy(labels,logits)) sess.close()
3、相对熵(relative entropy)
又称为KL散度(Kullback–Leibler divergence,简称KLD),信息散度(information divergence),信息增益(information gain)。
相对熵是交叉熵与信息熵的差值。即
相对熵=交叉熵-信息熵
KL(p||q)=H(p,q)-H(p)
表示用分布q模拟真实分布p相比用p模拟p,所需的额外信息。
相对熵(KL散度)有两个主要的性质。如下
(1)尽管 KL 散度从直观上是个度量或距离函数,但它并不是一个真正的度量或者距离,因为它不具有对称性,即
(2)相对熵具有非负性
总结一下:
信息熵公式:
交叉熵公式:
相对熵公式:
三者的关系:
在python中的实现:
import numpy as np import scipy.stats p=np.asarray([0.65,0.25,0.07,0.03]) q=np.array([0.6,0.25,0.1,0.05]) #方法一:根据公式求解 kl1=np.sum(p*np.log(p/q)) #方法二:调用scipy包求解 kl2=scipy.stats.entropy(p, q)
4、js散度(Jensen-Shannon)
因为kl散度不具对称性,因此js散度在kl散度的基础上进行了改进:
现有两个分布p1和p2,其JS散度公式为:
在python中的实现:
import numpy as np import scipy.stats p=np.asarray([0.65,0.25,0.07,0.03]) q=np.array([0.6,0.25,0.1,0.05]) M=(p+q)/2 #方法一:根据公式求解 js1=0.5*np.sum(p*np.log(p/M))+0.5*np.sum(q*np.log(q/M)) #方法二:调用scipy包求解 js2=0.5*scipy.stats.entropy(p, M)+0.5*scipy.stats.entropy(q, M)
1、f 散度(f-divergence)
KL-divergence 的坏处在于它是无界的。事实上KL-divergence 属于更广泛的 f-divergence 中的一种。
如果P和Q被定义成空间中的两个概率分布,则f散度被定义为:
一些通用的散度,如KL-divergence, Hellinger distance, 和total variation distance,都是f散度的一种特例。只是f函数的取值不同而也。
在python中的实现 :
import numpy as np import scipy.stats p=np.asarray([0.65,0.25,0.07,0.03]) q=np.array([0.6,0.25,0.1,0.05]) def f(t): return t*np.log(t) #方法一:根据公式求解 f1=np.sum(q*f(p/q)) #方法二:调用scipy包求解 f2=scipy.stats.entropy(p, q)
2、Hellinger distance
1 定义
1.1 度量理论
为了从度量理论的角度定义Hellinger距离,我们假设P和Q是两个概率测度,并且它们对于第三个概率测度λ来说是绝对连续的,则P和Q的Hellinger距离的平方被定义如下:
这里的dP / dλ 和 dQ / dλ分别是P和Q的Radon–Nikodym微分。这里的定义是与λ无关的,因此当我们用另外一个概率测度替换λ时,只要P和Q关于它绝对连续,那么上式就不变。为了简单起见,我们通常把上式改写为:
1.2 基于Lebesgue度量的概率理论
为了在经典的概率论框架下定义Hellinger距离,我们通常将λ定义为Lebesgue度量,此时dP / dλ 和 dQ / dλ就变为了我们通常所说的概率密度函数。如果我们把上述概率密度函数分别表示为 f 和 g ,那么可以用以下的积分形式表示Hellinger距离:
上述等式可以通过展开平方项得到,注意到任何概率密度函数在其定义域上的积分为1。
根据柯西-施瓦茨不等式(Cauchy-Schwarz inequality),Hellinger距离满足如下性质:
1.3 离散概率分布
对于两个离散概率分布 P=(p1,p2,...,pn)和 Q=(q1,q2,...,qn),它们的Hellinger距离可以定义如下:
上式可以被看作两个离散概率分布平方根向量的欧式距离,如下所示:
也可以写成:
在python中的实现:
import numpy as np p=np.asarray([0.65,0.25,0.07,0.03]) q=np.array([0.6,0.25,0.1,0.05]) #方法一: h1=1/np.sqrt(2)*np.linalg.norm(np.sqrt(p)-np.sqrt(q)) #方法二: h2=np.sqrt(1-np.sum(np.sqrt(p*q)))
3、巴氏距离(Bhattacharyya Distance)
在统计中,Bhattacharyya距离测量两个离散或连续概率分布的相似性。它与衡量两个统计样品或种群之间的重叠量的Bhattacharyya系数密切相关。Bhattacharyya距离和Bhattacharyya系数以20世纪30年代曾在印度统计研究所工作的一个统计学家A. Bhattacharya命名。同时,Bhattacharyya系数可以被用来确定两个样本被认为相对接近的,它是用来测量中的类分类的可分离性。
对于离散概率分布 p和q在同一域 X,巴氏距离被定义为:
其中BC(p,q)是Bhattacharyya系数:
对于连续概率分布,Bhattacharyya系数被定义为:
从公式可以看出,Bhattacharyya系数BC(P,Q)可以和前面的Hellinger距离联系起来,此时Hellinger距离可以被定义为:
因此,求得巴氏系数之后,就可以求得巴氏距离和Hellinger距离。
在python中的实现:
import numpy as np p=np.asarray([0.65,0.25,0.07,0.03]) q=np.array([0.6,0.25,0.1,0.05]) BC=np.sum(np.sqrt(p*q)) #Hellinger距离: h=np.sqrt(1-BC) #巴氏距离: b=-np.log(BC)
4、MMD距离(Maximum mean discrepancy)
最大均值差异(Maximum mean discrepancy),度量在再生希尔伯特空间中两个分布的距离,是一种核学习方法。两个随机变量的距离为:
其中k(.)是映射,用于把原变量映射到高维空间中。X,Y表示两种分布的样本,F表示映射函数集。
基于两个分布的样本,通过寻找在样本空间上的映射函数K,求不同分布的样本在K上的函数值的均值,通过把两个均值作差可以得到两个分布对应于K的mean discrepancy。寻找一个K使得这个mean discrepancy有最大值,就得到了MMD。最后取MMD作为检验统计量(test statistic),从而判断两个分布是否相同。如果这个值足够小,就认为两个分布相同,否则就认为它们不相同。更加简单的理解就是:求两堆数据在高维空间中的均值的距离。
近年来,MMD越来越多地应用在迁移学习中。在迁移学习环境下训练集和测试集分别取样自分布p和q,两类样本集不同但相关。我们可以利用深度神经网络的特征变换能力,来做特征空间的变换,直到变换后的特征分布相匹配,这个过程可以是source domain一直变换直到匹配target domain。匹配的度量方式就是MMD。
在python中的实现,根据核函数不同,公式可能不一样,根据公式编程即可。
5、Wasserstein distance
Wasserstein 距离,也叫Earth Mover‘s Distance,推土机距离,简称EMD,用来表示两个分布的相似程度。
Wasserstein distance 衡量了把数据从分布“移动成”分布时所需要移动的平均距离的最小值(类似于把一堆土从一个形状移动到另一个形状所需要做的功的最小值)
EMD是2000年IJCV期刊文章《The Earth Mover‘s Distance as a Metric for Image Retrieval》提出的一种直方图相似度量(作者在之前的会议论文中也已经提到,不过鉴于IJCV的权威性和完整性,建议参考这篇文章)。
假设有两个工地P和Q,P工地上有m堆土,Q工地上有n个坑,现在要将P工地上的m堆土全部移动到Q工地上的n个坑中,所做的最小的功。
每堆土我们用一个二元组来表示(p,w),p表示土堆的中心,w表示土的数量。则这两个工地可表示为:
每个土堆中心pi到每个土坑中心qj都会有一个距离dij,则构成了一个m*n的距离矩阵。
那么问题就是我们希望找到一个流(flow),当然也是个矩阵[fij],每一项fij代表从pi到qj的流动数量,从而最小化整体的代价函数:
问题描述清楚了:就是把P中的m个坑的土,用最小的代价搬到Q中的n个坑中,pi到qj的两个坑的距离由dij来表示。fij是从pi搬到qj的土的量;dij是pi位置到qj位置的代价(距离)。要最小化WORK工作量。EMD是把这个工作量归一化以后的表达,即除以对fij的求和。
EMD公式:
更多关于EMD的理解请参考:
http://blog.csdn.net/zhangping1987/article/details/25368183
在python中的实现:调用opencv
import numpy as np import cv #p、q是两个矩阵,第一列表示权值,后面三列表示直方图或数量 p=np.asarray([[0.4,100,40,22], [0.3,211,20,2], [0.2,32,190,150], [0.1,2,100,100]],np.float32) q=np.array([[0.5,0,0,0], [0.3,50,100,80], [0.2,255,255,255]],np.float32) pp=cv.fromarray(p) qq=cv.fromarray(q) emd=cv.CalcEMD2(pp,qq,cv.CV_DIST_L2)
最后计算出来的emd:
emd = 160.542770
概率分布之间的距离度量以及python实现