4. 现在，我们可以用类似的办法找到「距离方差」。请记住，若两个向量相同，其协方差与其方差相等。因此，距离方差可表示如下：
5. 最后，我们利用上述公式计算距离相关性。请记住，（距离）标准差与（距离）方差的平方根相等。
如果你更喜欢代码实现而非数学符号，那么请看下面的 R 语言实现：
set.seed(1234) doubleCenter <- function(x){ centered <- x for(i in 1:dim(x)[1]){ for(j in 1:dim(x)[2]){ centered[i,j] <- x[i,j] - mean(x[i,]) - mean(x[,j]) mean(x) } } return(centered) } distanceCovariance <- function(x,y){ N <- length(x) distX <- as.matrix(dist(x)) distY <- as.matrix(dist(y)) centeredX <- doubleCenter(distX) centeredY <- doubleCenter(distY) calc <- sum(centeredX * centeredY) return(sqrt(calc/(N^2))) } distanceVariance <- function(x){ return(distanceCovariance(x,x)) } distanceCorrelation <- function(x,y){ cov <- distanceCovariance(x,y) sd <- sqrt(distanceVariance(x)*distanceVariance(y)) return(cov/sd) } # Compare with Pearson's r x <- -10:10 y <- x^2 rnorm(21,0,10) cor(x,y) # --> 0.057 distanceCorrelation(x,y) # --> 0.509
任意两变量的距离相关性都在 0 和 1 之间。其中，0 代表两变量相互独立，而接近于 1 则表明变量间存在依赖关系。
如果你不想从头开始编写距离相关方法，你可以安装 R 语言的 energy 包（https://cran.r-project.org/web/packages/energy/index.html），设计此方案的研究者提供了本代码。在该程序包中，各类可用方案调用的是 C 语言编写的函数，因此有着很大的速度优势。
物理解释
关于距离相关性的表述，有着一个更令人惊讶的结果——它与布朗关联（Brownian correlation）有着确切的等价关系。
布朗关联指的是两个布朗过程之间的独立性（或依赖性）。相互依赖的布朗过程将会表现出彼此「跟随」的趋势。
让我们用一个简单的比喻来把握距离相关性的概念——请看下图中漂浮在湖面上的小纸船。
如果没有盛行风向，那么每艘船都将进行随机漂流——这与布朗运动类似。
无盛行风向时，小船随机漂流
如果存在盛行风向，那么小船漂流的方向将依赖于风的强度。风力越强，依赖性越显著。
有盛行风向时，小船倾向于同向漂流
与之类似，无关变量可以被看作无盛行风向时随机漂流的小船；相关变量可以被看作在盛行风向影响下漂流的小船。在这个比喻中，风的强弱就代表着两个变量之间相关性的强弱。
如果我们允许盛行风向在湖面的不同位置有所不同，那么我们就可以引入非线性的概念。距离相关性利用「小船」之间的距离推断盛行风的强度。
置信区间？
我们可以采取「重采样（resampling）」方法为距离相关性估计建立置信区间。一个简单的例子是 bootstrap 重采样。
这是一个巧妙的统计技巧，需要我们从原始数据集中随机抽样（替换）以「重建」数据。这个过程将重复多次（例如 1000 次），每次都计算感兴趣的统计量。
这将为我们感兴趣的统计量产生一系列不同的估计值。我们可以通过它们估计在给定置信水平下的上限和下限。
请看下面的 R 语言代码，它实现了简单的 bootstrap 函数：
set.seed(1234) bootstrap <- function(x,y,reps,alpha){ estimates <- c() original <- data.frame(x,y) N <- dim(original)[1] for(i in 1:reps){ S <- original[sample(1:N, N, replace = TRUE),] estimates <- append(estimates, distanceCorrelation(S$x, S$y)) } u <- alpha/2 ; l <- 1-u interval <- quantile(estimates, c(l, u)) return(2*(dcor(x,y)) - as.numeric(interval[1:2])) } # Use with 1000 reps and threshold alpha = 0.05 x <- -10:10 y <- x^2 rnorm(21,0,10) bootstrap(x,y,1000,0.05) # --> 0.237 to 0.546
如果你想建立统计显著性，还有另一个重采样技巧，名为「排列检验（permutation test）」。
排列检验与上述 bootstrap 方法略有不同。在排列检验中，我们保持一个向量不变，并通过重采样对另一个变量进行「洗牌」。这接近于零假设（null hypothesis）——即，在变量之间不存在依赖关系。
这个经「洗牌」打乱的变量将被用于计算它和常变量间的距离相关性。这个过程将被执行多次，然后，结果的分布将与实际距离相关性（从未被「洗牌」的数据中获得）相比较。
然后，大于或等于「实际」结果的经「洗牌」的结果的比例将被定为 P 值，并与给定的显著性阈值（如 0.05）进行比较。
以下是上述过程的代码实现：
permutationTest <- function(x,y,reps){ estimates <- c() observed <- distanceCorrelation(x,y) N <- length(x) for(i in 1:reps){ y_i <- sample(y, length(y), replace = T) estimates <- append(estimates, distanceCorrelation(x, y_i)) } p_value <- mean(estimates >= observed) return(p_value) } # Use with 1000 reps x <- -10:10 y <- x^2 rnorm(21,0,10) permutationTest(x,y,1000) # --> 0.036
最大信息系数
最大信息系数（MIC）于 2011 年提出，它是用于检测变量之间非线性相关性的最新方法。用于进行 MIC 计算的算法将信息论和概率的概念应用于连续型数据。
深入细节
由克劳德·香农于 20 世纪中叶开创的信息论是数学中一个引人注目的领域。
信息论中的一个关键概念是熵——这是一个衡量给定概率分布的不确定性的度量。概率分布描述了与特定事件相关的一系列给定结果的概率。
概率分布的熵是「每个可能结果的概率乘以其对数后的和」的负值
为了理解其工作原理，让我们比较下面两个概率分布：
X 轴标明了可能的结果；Y 轴标明了它们各自的概率
左侧是一个常规六面骰子结果的概率分布；而右边的六面骰子不那么均匀。
从直觉上来说，你认为哪个的熵更高呢？哪个骰子结果的不确定性更大？让我们来计算它们的熵，看看答案是什么。
entropy <- function(x){ pr <- prop.table(table(x)) H <- sum(pr * log(pr,2)) return(-H) } dice1 <- 1:6 dice2 <- c(1,1,1,1,2:6) entropy(dice1) # --> 2.585 entropy(dice2) # --> 2.281
不出所料，常规骰子的熵更高。这是因为每种结果的可能性都一样，所以我们不会提前知道结果偏向哪个。但是，非常规的骰子有所不同——某些结果的发生概率远大于其它结果——所以它的结果的不确定性也低一些。
这么一来，我们就能明白，当每种结果的发生概率相同时，它的熵最高。而这种概率分布也就是传说中的「均匀」分布。
交叉熵是熵的一个拓展概念，它引入了第二个变量的概率分布。
crossEntropy <- function(x,y){ prX <- prop.table(table(x)) prY <- prop.table(table(y)) H <- sum(prX * log(prY,2)) return(-H) }
两个相同概率分布之间的交叉熵等于其各自单独的熵。但是对于两个不同的概率分布，它们的交叉熵可能跟各自单独的熵有所不同。
这种差异，或者叫「散度」可以通过 KL 散度（Kullback-Leibler divergence）量化得出。
两概率分布 X 与 Y 的 KL 散度如下：
概率分布 X 与 Y 的 KL 散度等于它们的交叉熵减去 X 的熵
KL 散度的最小值为 0，仅当两个分布相同。
KL_divergence <- function(x,y){ kl <- crossEntropy(x,y) - entropy(x) return(kl) }
为了发现变量具有相关性，KL 散度的用途之一是计算两个变量的互信息（MI）。
互信息可以定义为「两个随机变量的联合分布和边缘分布之间的 KL 散度」。如果二者相同，MI 值取 0。如若不同，MI 值就为一个正数。二者之间的差异越大，MI 值就越大。
为了加深理解，我们首先简单回顾一些概率论的知识。
变量 X 和 Y 的联合概率就是二者同时发生的概率。例如，如果你抛掷两枚硬币 X 和 Y，它们的联合分布将反映抛掷结果的概率。假设你抛掷硬币 100 次，得到「正面、正面」的结果 40 次。联合分布将反映如下：
P(X=H, Y=H) = 40/100 = 0.4
jointDist <- function(x,y){ N <- length(x) u <- unique(append(x,y)) joint <- c() for(i in u){ for(j in u){ f <- x[paste0(x,y) == paste0(i,j)] joint <- append(joint, length(f)/N) } } return(joint) }
边缘分布是指不考虑其它变量而只关注某一特定变量的概率分布。假设两变量独立，二者边缘概率的乘积即为二者同时发生的概率。仍以抛硬币为例，假如抛掷结果是 50 次正面和 50 次反面，它们的边缘分布如下：
P(X=H) = 50/100 = 0.5 ; P(Y=H) = 50/100 = 0.5
P(X=H) × P(Y=H) = 0.5 × 0.5 = 0.25
marginalProduct <- function(x,y){ N <- length(x) u <- unique(append(x,y)) marginal <- c() for(i in u){ for(j in u){ fX <- length(x[x == i]) / N fY <- length(y[y == j]) / N marginal <- append(marginal, fX * fY) } } return(marginal) }
现在让我们回到抛硬币的例子。如果两枚硬币相互独立，边缘分布的乘积表示每个结果可能发生的概率，而联合分布则为实际得到的结果的概率。
如果两硬币完全独立，它们的联合概率在数值上（约）等于边缘分布的乘积。若只是部分独立，此处就存在散度。
这个例子中，P(X=H,Y=H) > P(X=H) × P(Y=H)。这表明两硬币全为正面的概率要大于它们的边缘分布之积。
联合分布和边缘分布乘积之间的散度越大，两个变量之间相关的可能性就越大。两个变量的互信息定义了散度的度量方式。
X 和 Y 的互信息等于「二者边缘分布积和的联合分布的 KL 散度」
mutualInfo <- function(x,y){ joint <- jointDist(x,y) marginal <- marginalProduct(x,y) Hjm <- - sum(joint[marginal > 0] * log(marginal[marginal > 0],2)) Hj <- - sum(joint[joint > 0] * log(joint[joint > 0],2)) return(Hjm - Hj) }
此处的一个重要假设就是概率分布是离散的。那么我们如何把这些概念应用到连续的概率分布呢？
分箱算法
其中一种方法是量化数据（使变量离散化）。这是通过分箱算法（bining）实现的，它能将连续的数据点分配对应的离散类别。
此方法的关键问题是到底要使用多少「箱子（bin）」。幸运的是，首次提出 MIC 的论文给出了建议：穷举！
也就是说，去尝试不同的「箱子」个数并观测哪个会在变量间取到最大的互信息值。不过，这提出了两个挑战：
- 要试多少个箱子呢？理论上你可以将变量量化到任意间距值，可以使箱子尺寸越来越小。
- 互信息对所用的箱子数很敏感。你如何公平比较不同箱子数目之间的 MI 值？
第一个挑战从理论上讲是不能做到的。但是，论文作者提供了一个启发式解法（也就是说，解法不完美，但是十分接近完美解法）。他们也给出了可试箱子个数的上限。
最大可用箱子个数由样本数 N 决定
至于如何公平比较取不同箱子数对 MI 值的影响，有一个简单的做法……就是归一化！这可以通过将每个 MI 值除以在特定箱子数组合上取得的理论最大值来完成。我们要采用的是产生最大归一化 MI 总值的箱子数组合。
互信息可以通过除以最小的箱子数的对数来归一化
最大的归一化互信息就是 X 和 Y 的最大信息系数（MIC）。我们来看看一些估算两个连续变量的 MIC 的代码。
MIC <- function(x,y){ N <- length(x) maxBins <- ceiling(N ** 0.6) MI <- c() for(i in 2:maxBins) { for (j in 2:maxBins){ if(i * j > maxBins){ next } Xbins <- i; Ybins <- j binnedX <-cut(x, breaks=Xbins, labels = 1:Xbins) binnedY <-cut(y, breaks=Ybins, labels = 1:Ybins) MI_estimate <- mutualInfo(binnedX,binnedY) MI_normalized <- MI_estimate / log(min(Xbins,Ybins),2) MI <- append(MI, MI_normalized) } } return(max(MI)) } x <- runif(100,-10,10) y <- x**2 rnorm(100,0,10) MIC(x,y) # --> 0.751
以上代码是对原论文中方法的简化。更接近原作的算法实现可以参考 R package minerva（https://cran.r-project.org/web/packages/minerva/index.html）。
在 Python 中的实现请参考 minepy module（https://minepy.readthedocs.io/en/latest/）。
MIC 能够表示各种线性和非线性的关系，并已得到广泛应用。它的值域在 0 和 1 之间，值越高表示相关性越强。
置信区间？
为了建立 MIC 估计值的置信区间，你可以简单地使用一个像我们之前介绍过的 bootstrap 函数。我们可以利用 R 语言的函数式编程，通过传递我们想要用作参数的函数来泛化 bootstrap 函数。
bootstrap <- function(x,y,func,reps,alpha){ estimates <- c() original <- data.frame(x,y) N <- dim(original)[1] for(i in 1:reps){ S <- original[sample(1:N, N, replace = TRUE),] estimates <- append(estimates, func(S$x, S$y)) } l <- alpha/2 ; u <- 1 - l interval <- quantile(estimates, c(u, l)) return(2*(func(x,y)) - as.numeric(interval[1:2])) } bootstrap(x,y,MIC,100,0.05) # --> 0.594 to 0.88
总结
为了总结相关性这一主题，我们来测试下各算法在人工生成数据上的处理能力。
完整代码：https://gist.github.com/anonymous/fabecccf33f9c3feb568384f626a2c07
噪声函数
set.seed(123) # Noise x0 <- rnorm(100,0,1) y0 <- rnorm(100,0,1) plot(y0~x0, pch = 18) cor(x0,y0) distanceCorrelation(x0,y0) MIC(x0,y0)
- Pearson's r = - 0.05
- 距离相关性 = 0.157
- MIC = 0.097
简单线性函数
