相依性度量在fMRI数据分析, 基因数据分析, 信号处理, 因果推断, 图像识别等领域都起着重要作用.相依性分析可以说是统计分析的第一步, 因此引起了许多学者的注意.最早提出的相依性度量是Pearson相关系数(PCC), PCC简单易算且高效.随着应用的发展和研究的深入, 人们渐渐地发现PCC的不足, 比如它只能识别线性相关函数, 而且对异常点非常敏感.为了弥补这些缺点, 人们相继提出了别的相依性度量:最大相关系数^[5](Maximal Correlation Coefficient, MCC)、互信息^{[2, 6-7]}(Mutual Information, MI)、核互信息^[8](Kernel Mutual Information, KMI)、广义方差(the Generalized Variance)、标准相关系数(Canonical Correlation Analysis, CCA)、核标准相关系数^[9](Kernel Canonical Correlation Analysis, KCCA)、Hilber-Schmidt独立性指数^[10](Hilber-Schmidt Independence Criterion, HSIC)、Schweizer-Wolff度量^[11](Schweizer-Wolff measure)、距离相关系数^[12](Distance based correlation, dcor)、最大平均散度^[13](Maximum Mean Discrepancy, MMD)、耦合最大平均散度^[14](Copula Maximum Mean Discrepancy, CMMD)、HHG^[3]、基于信息的标准相关系数^[15](Canonical Correlation Analysis based on Information, ICCA)、曲线相关系数^[16](correlation coefficient for variable distributed along a curve, curve correlation)、随机化相关系数^[17](Randomized Dependence Coefficient, RDC)、耦合相关系数^[4](Copula Dependence Coefficient, CDC)等等.据我们了解, 目前并没有文献给出这些方法的综合评价, 本文从检验功效, 稳健性等角度分析这些方法的优劣性, 是目前为止最为完整的综述性文章.

度量两个变量之间的相依性具有巨大的挑战性, 学者们并不满足于现有的方法^[7]. Linfoot^[18]首先给出了相关系数PCC和互信息MI之间的关系, 这对以后的研究产生了深远的影响.从那时起, 信息理论就成为了相依性度量中的基本工具. Rényi- $\alpha$^[5]和Tsallis- $\alpha$ ^[19]互信息就是继PCC后最受关注的相依性度量.互信息具有一些很好的性质, 比如对任意的严格单调函数 $f$和 $G$, 有 $MI(X, Y)=MI(f(X), G(Y))$.然而不同变量之间的互信息并不可比较.假设有三个变量 $X, Y, Z$是我们所关心的, 那我们可能就会问 $X$和 $Y$之间的关联比 $X$和 $Z$之间的关联强吗? PCC并不能回答这样的问题, 因为两个变量之间的PCC等于0, 并不意味着这两个变量就是相互独立的.互信息也无法回答这个问题, 因为互信息依赖于变量本身的信息，由此学者们提出了最大信息系数^[1](Maximal Information Coefficient, MIC), MIC本质上就是正则化的互信息.

Renyi^[5]指出相依性度量 $D(X, Y)$应该具有以下性质: (A) $D(X, Y)$对于任意的非常数随机变量 $X, Y$都有定义; (B) $D(X, Y)=D(Y, X)$; (C) $0\leq D(X, Y)\leq 1$; (D) $D(X, Y)=0$, 当且仅当 $X$和 $Y$相互独立; (E) $D(X, Y)=1$, 当且仅当存在一个函数 $f$, 满足 $Y=f(X)$; (F)对于Borel可测函数 $f(x)$和 $g(x)$, 我们有 $D(f(X), g(Y))=D(X, Y)$; (F')对于单调的 $f(x)$和 $g(x)$, 我们有 $D(f(X), g(Y))=D(X, Y)$; (G) $D(X, Y)=|\rho(X, Y)|$, 当 $X, Y$的联合分布是正态分布; 因此满足这些条件的相依性度量就是一个完美的相依性度量.显然, 相关性系数只符合条件(A), (B), (C), (G).条件(F)其实非常严格, 通常将函数 $f(x)$和 $g(x)$弱化到严格单调函数即性质(F').除了这些性质之外, 我们还希望它是一致的, 稳健的, 适用于高维情况.把相依性度量从二维推广到高维非常困难.比如计算高维变量之间的互信息, 需要先估计高维概率密度函数.而核密度估计方法有维数灾问题.

在很多情况下, 我们通常想知道偏相关性, 即:去除第三个变量的影响之后, 两个变量之间的相关性.对于PCC, 我们有偏相关系数.偏相关系数继承了相关系数的所有缺点.因此同样的, 我们也希望能在这个应用中引入其他的相依性度量来克服这些缺点.蒋杭进等^[20]给出了偏相关性度量的统一框架, 在这个框架下, 所有的相依性度量都可以推广到偏相依性度量.他们还给出了偏MIC在变量选择中的应用.因此在本文中, 我们将略过此方面的讨论.

本文的内容安排:第二部分详细分析了相关系数的性质, 并逐步引入与之相关的相依性度量, 对这一类相依性度量进行了比较分析.文章的第三部分主要介绍基于独立性条件的相依性度量, 并对这一类相依性度量进行了比较分析.第四部分则介绍基于学习理论的相依性度量.第五部分综合讨论所有的相依性度量和非线性度量问题.

相关系数(PCC)是最早提出的相依性度量.它因计算简单方便且容易解释而被广泛应用于各个领域的数据分析中, 比如fMRI, Omics数据分析, 基因数据分析. PCC的定义为

$ \begin{equation} \rho(X, Y)=\frac{cov(X, Y)}{\sqrt{var(X)var(Y)}}, \end{equation} $

(2.1)

其中 $cov(X, Y)$是随机变量 $X, Y$的协方差, $var(X)$和 $var(Y)$分别是随机变量 $X$和 $Y$方差.显然, $\rho(X, Y)\in [-1, 1]$.通常, 我们有

a. $|\rho(X, Y)|=1$意味着 $X, Y$是线性相关的;

b. $\rho(X, Y)=0$意味着 $X, Y$不是线性相关的;

c. $\rho(X, Y)<0$意味着 $X, Y$是线性负相关/相互抑制的;

d. $\rho(X, Y)>0$意味着 $X, Y$是线性正相关/相互促进的.

然而, 在某些情况下, (c)和(d)并不正确.我们将给出一个例子来说明.

假设有两种不同的食品: A和B, 它们的价格由一个趋势成分和另一个随机波动决定. A和B的随机波动是一样的(因为相互促进), 但是A和B的趋势成分是相反的.也就是说A的价格增加多少, B的价格就会相应的降低多少(见图 1).从图中可以看出A的价格和B的价格是负相关的.但是A的波动价格和B的波动价格则是正相关的.换句话说相关系数只反应了一个整体趋势的相关性, 并不能反应局部的相互作用.也就是说相互作用关系并不等同于相关性.相互作用的结果是相关量的改变, 在我们的例子中就是价格的变化, 因此我们应该基于价格的变化量来度量相互作用强度.分别记 $PA_t$和 $PB_t$为A和B在时刻 $t$的价格, 他们之间的相互作用强度定义为 $B(PA_t)$和 $B(PB_t)$之间的相关性, 即

$ \begin{equation} IC(x_t, y_t)=\rho(B(x_t), B(y_t)), \end{equation} $

(2.2)

其中 $B(x_t)=x_t-x_{t-1}$.

在我们的例子中 $\rho(A, B)=-0.945$, 然而IC $(A, B)=1$.

在引言中, 我们提到相关系数有以下缺点: (1)容易受异常点影响; (2)无法检测非线性函数关系.

相关系数是不稳健的.令 $X=\{x_1, x_2, \cdots, x_{100}\}\sim U(0, 1)$, $Y=\{y_1, y_2, \cdots, y_{100}\}$ $\sim U(0, 1)$, 且 $X$和 $Y$相互独立, 我们有 $\rho(X, Y)=0.0129\pm 0.10$.如果我们将 $X$和 $Y$的第1个值替换为 $t$, 即, $X[1]=Y[1]=t$, 随着 $t$由1变到100, $\rho(X, Y)$则由0变到1, 正如图 2所示.

图 1中, 价格A和价格B之间是负相关的, 相关系数为 $\rho(A, B)=-0.945$, 但是它们之间是相互促进的, 相互作用强度IC $(A, B)=1$.

图 2中, 我们取 $X$为 $[0, 1]$上均匀分布, 样本长度为100, $Y$同样是长度为100的均匀分布随机样本, 即: $X=\{x_1, x_2, \cdots, x_{100}\} \sim U(0, 1)$, $Y=\{y_1, y_2, \cdots, y_{100}\} \sim U(0, 1)$, 且 $X$和 $Y$相互独立, 我们有 $\rho(X, Y)=0.0129\pm 0.10$.如果我们将 $X$和 $Y$的第1个替换为 $t$, 随着 $t$由1变到100, $\rho(X, Y)$则由0变到1.

相关系数无法识别一些非线性函数.取 $X=1:100$, $Y=(X-50)^2$, 有 $\rho(X, Y)=0.0387$ ( $p$值=0.7022).然而二次函数关系是常见的, 这时相关系数给我们一个错误的结论.因此, 我们就会问什么时候用相关系数呢?以下定理回答了这个问题.

定理2.1  (Lancaster) 假设 $X$和 $Y$的联合分布是正态的, 且它们之间的相关系数为 $\rho$, 变换 $X'=X'(X)$, $Y'=Y'(Y)$将 $X, Y$变为新的随机变量, 且 $var(X')$, 和 $var(Y')$都是有限的, 那么新的随机变量之间的相关系数小于等于 $\rho$, 也就是说 $\rho(X', Y')^2 \leq \rho^2$.

假设有 $n$个观察值 $\{(x_{i}, y_{i}):i=1, 2, \cdots, n\}$.我们不失一般性地假设 $x_{1}<x_{2}<\cdots <x_{n}$, { $x_{1}, x_{2}, \cdots, x_{n}$}和{ $y_{1}, y_{2}, \cdots, y_{n}$}的秩分别为 $R_{x}=${ $1, 2, 3, \cdots, n$}和 $R_{y}=${ $R_{1}, R_{2}, \cdots, R_{n}$}. Spearman秩相关系数(Scor)被定义为 $R_{x}$和 $R_{y}$的相关系数, 即

$ \begin{eqnarray*} S &=& \frac{\sum(i-(n+1)/2)(R_{i}-(n+1)/2)}{\sqrt{\sum(i-(n+1)/2)^{2}\sum(R_{i}-(n+1)/2)^{2}}} \nonumber\\ &=& 1-\frac{6}{n(n^{2}-1)}\sum(i-R_{i})^{2} . \end{eqnarray*} $

类似于Scor, Kendall's $\tau$ (Kcor)统计量定义为

$ \begin{equation} K=1-\frac{4}{n(n-1)}\sum\limits_{i<j}\phi(R_{i}-R_{j}), \end{equation} $

(2.3)

其中, $\phi(x)$=1当 $x>0$, $\phi(x)=-1$当 $x<0$, $\phi(x)=0$当 $x=0$.

可以从定义看出这两种相依性度量都无法检测出具有周期性的函数关系, 如 $Y=\sin(X)$.

曲线相关系数^[16]用于度量两个沿着曲线分布的随机变量间的相关性.具体而言, 我们假设 $(X, Y)$是二维随机变量, $\Sigma$是其协方差矩阵.由矩阵的谱分解得到

$ \begin{equation} \Sigma= \Bigg( \begin{array}{cc} \sigma_{X}^2 ~& \sigma_{XY} \\ \sigma_{YX} ~& \sigma_{Y}^2 \end{array} \Bigg) =\Bigg( \begin{array}{cc} \cos\alpha~&-\sin\alpha \\ \sin\alpha~&\cos\alpha\\ \end{array} \Bigg)^T \Bigg( \begin{array}{cc} \lambda_1~&0 \\ 0 ~& \lambda_2\\ \end{array} \Bigg) \Bigg( \begin{array}{cc} \cos\alpha ~& -\sin\alpha \\ \sin\alpha ~& \cos\alpha\\ \end{array} \Bigg), \end{equation} $

(2.4)

其中 $\alpha$是第一主成分和横轴之间的夹角.另外我们知道第一主成分和第二主成分之间是垂直, 从而

$ \begin{equation}\label{vars1} V(X)=\sigma_X^2=\lambda_1\cos^2\alpha+\lambda_2\sin^2\alpha, \end{equation} $

(2.5)

$ \begin{equation}\label{vars2} V(Y)=\sigma_Y^2=\lambda_1\sin^2\alpha+\lambda_2\cos^2\alpha, \end{equation} $

(2.6)

$ \begin{equation}\label{vars3} cov(X, Y)=\sigma_{XY}=(\lambda_1-\lambda_2)\cos\alpha\sin\alpha, \end{equation} $

(2.7)

$ \begin{equation} \rho(X, Y)=\frac{\sigma_{XY}}{\sigma_X\sigma_Y}=\frac{(\lambda_1-\lambda_2)\cos\alpha\sin\alpha}{(\lambda_1\cos^2\alpha+\lambda_2\sin^2\alpha)^{\frac{1}{2}}(\lambda_1\sin^2\alpha+\lambda_2\cos^2\alpha)^{\frac{1}{2}}}. \end{equation} $

(2.8)

平面中的随机变量 $(X, Y)$可以表示为曲线 $c(I): I \in {\Bbb R } \mapsto {\Bbb R}^2$的点加上垂直的随机噪声, 即对任意的 $s \in I$, $(X, Y)=\chi_c(S, T)=c(s)+tv(s)$, 其中 $c'(s)^Tv(s)=0$, $c'(s)^Tc(s)=1$, $c(s)$是曲线的弧长参数形式.类似于相关系数, 令 $\alpha(s)$为 $c'(s)$和坐标轴之间的夹角, 类比于(2.5)式, (2.6)式和(2.7)式, 我们就得到了局部方差和局部协方差.

$ \begin{equation} LV_{X}(s)=V(S)\cos^2\alpha(s)+V(T|S=s)\sin^2\alpha(s), \end{equation} $

(2.9)

$ \begin{equation} LV_{Y}(s)=V(S)\sin^2\alpha(s)+V(T|S=s)\cos^2\alpha(s), \end{equation} $

(2.10)

$ \begin{equation} LCov_{(X, Y)}(s)=(V(S)-V(T|S=s))\cos\alpha(s)\sin\alpha(s), \end{equation} $

(2.11)

现在我们就可以定义在 $c(s)$处的局部相关系数

$ \begin{equation} LCor_{(X, Y)}(s)=\frac{LCov_{(X, Y)}(s)}{\sqrt{LV_{Y}(s)LV_{X}(s)}} . \end{equation} $

(2.12)

因此曲线相关系数定义为

$ \begin{equation} CorGC(X, Y)=\sqrt{E_{S}[LCor_{(X, Y)}(S)^2]}. \end{equation} $

(2.13)

曲线相关系数的计算是基于主曲线^[16], 但是有的情形下无法计算得到主曲线, 从而无法得到曲线相关系数.

标准相关系数是两组变量之间的Pearson相关系数.令 $X=(X_1, X_2, \cdots, X_n)$, $Y=(Y_1, Y_2, \cdots, Y_m)$, 标准相关系数定义为

$ \begin{equation} CCA(X, Y)= \sup\limits_{\alpha, \beta} \ \ \rho(\alpha'X, \beta'Y)=\sup\limits_{\alpha, \beta} \frac{\alpha'\Sigma_{XY}\beta}{\sqrt{\alpha'\Sigma_{XX}\alpha}\sqrt{\beta'\Sigma_{YY}\beta}}, \end{equation} $

(2.14)

其中, $X_i, i=1, 2, \cdots, n$和 $Y_j, j=1, 2, \cdots, m$是随机变量, $\Sigma_{XY}=cov(X, Y)$, $\Sigma_{XX}=cov(X, X)$, $\Sigma_{YY}=cov(Y, Y)$.

更进一步, 我们有

$ \begin{equation} CC(X, Y)= \lambda(\Sigma_{XX}^{-\frac{1}{2}} \Sigma_{XY} \Sigma_{YY}^{-1}\Sigma_{YX}\Sigma_{XX} ^{-\frac{1}{2}} ), \end{equation} $

(2.15)

其中 $\lambda(A)$是矩阵 $A$的最大特征值.本质上, $CC(X, Y)$是以下问题的解.

$ \begin{equation}\label{gep} \left( \begin{array}{cc} 0&\Sigma_{XY} \\ \Sigma_{YX}&0\\ \end{array} \right)\beta = \rho \left( \begin{array}{cc} \Sigma_{XX}&0 \\ 0&\Sigma_{YY} \\ \end{array} \right)\beta. \end{equation} $

(2.16)

ICCA定义为

$ \begin{equation} ICCA(X, Y)=\rho(\alpha'X, \beta'Y), \end{equation} $

(2.17)

其中

$ \begin{equation} (\alpha, \beta)= \arg\max \sup\limits_{a, b} E \log \frac{p(a'X, b'Y)}{p(a'X)p(b'Y)}. \end{equation} $

(2.18)

关于ICCA的详细参见文献[15].

假设 $x_1$和 $x_2$是 $\chi={\Bbb R}^p$空间中的随机变量, 令 $K_1$和 $K_2$是Mercer核, $\Phi_1$和 $\Phi_2$分别是 $K_1$和 $K_2$的特征映射. $\Phi_1(x_1)$与 $\Phi_2(x_2)$核相关系数 $\rho_{F}$定义为

$ \begin{equation} \rho_{F}=\max\limits_{(f_1, f_2)\in F_1 \times F_2}\rho(\langle \Phi_1(x_1), f_1\rangle, \langle \Phi_2(x_2), f_2\rangle), \end{equation} $

(2.19)

其中, $\langle \Phi_1(x_1), f_1\rangle = f_1(x_1)$和 $\langle \Phi_2(x_2), f_2\rangle = f_2(x_2)$.

KC还可以推广到核标准相关系数(KCCA), 参见文献[9, 21].

MCC^[22]定义为

$ \begin{equation} MCC(X, Y)=\sup\limits_{f, g \in L^2}\{\rho(f(X), g(Y))\}, \end{equation} $

(2.20)

其中, $f$和 $g$是borel可测函数.从KC的定义可以看出, 其实KC是一种特殊的MCC: KC的定义中, $L^2$空间取为再生核希尔伯特空间.因此我们有 $KC(X, Y) \leq MCC(X, Y)$.然而MCC并不具有性质(E), 因为MCC也很容易受到异常点的影响.为了克服这一困难, 蒋杭进等提出了耦合相依系数(Copula dependence coefficient, CDC)^[4].

另外计算MCC是非常困难的, 学者们给出了一种计算方法: Alternating Conditional Expectations (ACE)^[23], 更多关于它的讨论参见文献[22, 24-25].

正如前面所介绍, 耦合相依系数(Copula dependence coefficient, CDC)的提出是为了弥补MCC不稳健的缺点^[4].下面先简单介绍一下几个基本的定理.

定理2.2^[26]  设随机变量 $X$的累积分布函数为 $F$, 那么随机变量 $U=F(X)$的分布函数是 $[0, 1]$的均匀分布.

定理2.3^[26]  设随机变量向量 $X=(X_1, X_2, \cdots, X_p)$且有边缘分布函数 $F_i$, $1\leq i \leq p$.那么 $X$的联合累积分布由以下式子给出

$ \begin{equation} F(X_1, X_2, \cdots, X_p)=C_X(F_1(X_1), F_2(X_2), \cdots, F_k(X_p)), \end{equation} $

(2.21)

其中 $C_X$是 $X$的耦合函数.

经验分布是单变量的累积分布函数 $F_1, \cdots, F_p$的无偏估计

$ \begin{equation} F_{n, j}(x)=\frac{1}{n} \sum\limits_{i=1}^{n}I(X_i^j \leq x)\ \ (j=1, 2, \cdots, p). \end{equation} $

(2.22)

令

$ \begin{equation} [U_{1}^{i}, U_{2}^{i}, \cdots, U_{p}^{i}]=[F_{n, 1}(X_1^{i}), F_{n, 2}(X_2^{i}), \cdots, F_{n, p}(X_p^{i})], \end{equation} $

(2.23)

那么多变量的耦合分布为

$ \begin{equation} C_X^{n}(u_1, u_2, \cdots, u_p)=\frac{1}{n}\sum\limits_{i=1}^{n} I(U_{1}^{i}\leq u_1, U_{2}^{i}\leq u_2, \cdots, U_{p}^{i}\leq u_p). \end{equation} $

(2.24)

给定随机向量 $X=(X_1, X_2, \cdots, X_p)$, $Y=(Y_1, Y_2, \cdots, Y_q)$, 根据定理, 我们有

$ \begin{equation} F(X, Y)=C_{XY}(F_X(X), F_Y(Y)), \end{equation} $

(2.25)

其中 $F_X$和 $F_Y$分别是 $X$和 $Y$的边缘分布.

定理2.4^[26]  假设 $X$和 $Y$是连续的随机变量.那么 $X$和 $Y$相互独立当且仅当

$ \begin{equation} C_{XY}(F_X(X), F_Y(Y))=F_X(X)F_Y(Y). \end{equation} $

(2.26)

这样就引出了CDC^[4]的定义

$ \begin{equation}\label{cdc} CDC(X, Y)=MCC(F_X(X), F_Y(Y)) . \end{equation} $

(2.27)

本质上而言, CDC就是稳健版的MCC, 因此CDC继承了MCC的优良性质, 而且弥补了它的不足, 具体内容参见文献[4].

在前面的章节中, 我们给出了Spearman秩相关系数, Kendall's $\tau$, 曲线相关系数, CCA, ICCA, KC, KCCA和MCC的定义, 并简要地分析了它们的性质.本节中, 我们先比较这些方法, 到文章最后我们将给出更详细的讨论.

为了对这些方法的检验功效有个大致的了解, 我们用6种不同的函数类型做了模拟分析:线性(Linear), 二次型(Quadratic), 三次函数(Cubic), Sine, Cosine和Power 0.25, 具体的函数表达式见附录.

我们从均匀分布 $U(0, 1)$中抽取500个点作 $X$, 并令 $Y=f(X)+L\times N(0, 1)$, 这里 $L$是噪声水平.给定 $L$, 我们就得到了一组观察数据 $(X, Y)$, 计算其相依性的值 $DM_{obs}(X, Y)$.为了计算它的p值, 我们需要按照以下步骤通过模拟得到, 在噪声水平 $L$下, 各个相依性度量在原假设下的分布: (1)从均匀分布 $U(0, 1)$中抽取500个点作 $X$; (2) $Y=f(X)+L\times N(0, 1)$; (3)重新从均匀分布 $U(0, 1)$中抽取500个点作 $X'$; (4)计算观察数据 $(X', Y)$的相依性的值.重复步骤(1)-(4), 并得到了噪声水平为 $L$时, 该相依性度量在原假设的分布 $F_0$.详细的步骤在下面给出.

模拟算法A:计算相依性度量 $DM$在噪声水平 $L$下, 对函数类型 $f$的检验功效

1.给定噪声水平 $L$, 计算相依性度量 $DM$在原假设下的经验分布.

a1.从均匀分布 $U(0, 1)$中抽取500个点作 $X$.

a2. $Y=f(X)+L\times N(0, 1)$.

a3.重新从均匀分布 $U(0, 1)$中抽取500个点作 $X'$.

a4.计算观察数据 $(X', Y)$的相依性的值.

a5.重复步骤(a1)-(a4) 500次, 并得到了噪声水平为 $L$时, 该相依性度量在原假设的分布 $F_0$.

2.计算噪声水平 $L$下, 相依性度量 $DM$的检验功效

b1.从均匀分布 $U(0, 1)$中抽取500个点作 $X$.

b2. $Y=f(X)+L \times N(0, 1)$.

b3.计算其相依性的值 $DM_{obs}(X, Y)$.

b4.计算 $p$值, $1-F_0(DM_{obs}(X, Y))$.

b5.重复步骤(b1)-(b4) 500次.在置信水平0.05下, 计算拒绝原假设的次数 $N$.检验功效即为 $N/500$.

模拟结果见图 3.对于线性函数, 从图 3可以看出, pcor(Pearson相关系数)的功效是最高的, 这和定理2.1是一致的.对于二次函数, scor(Spearman秩相关系数), kcor(Kentall's $\tau$), pcor的检验功效都很低.正如前文讨论的那样, Pearson相关系数无法识别一些非线性函数关系, 而scor, kcor只能很好地识别单调函数类, 对于非单调的函数, 他们的检验功效也是很低的.对于这种函数类型, 曲线相关系数(curvecor)的功效比pcor, scor, kcor要高一些, 但是要低于最大相关系数(MCC/ACE).对于三次函数, Sine, Cosine, 我们有类似的结论.对于Power 0.25, 模拟结果看上去很奇怪.首先这是一个单调函数, 而且是一个可以用线性函数很好逼近的单调函数, 因此pcor, scor, kcor的检验功效非常类似, 然而MCC/ACE的表现却是差强人意.

总而言之, 如果存在一个线性函数 $g(x)=ax+b$可以很好地逼近函数 $f(x)$, 那么pcor, scor, kcor的检验功效是最高的.相比之下曲线相关系数的表现却显得很鸡肋, 因为从模拟结果上看, 相比于其他的方法它并没有明显的优势.

我们称随机变量 $X$和 $Y$是相互独立的, 当且仅当以下条件中有一条成立.

C1.对任意的集合 $A$和 $B$有 $P(\{X \in A\} \cap\{ Y \in B\})= $ $P(X \in A) P(Y \in B)$.

C2.对任意的Borel函数 $f$和 $g$有, $E[f(X)g(Y)]=E[f(X)]E[g(Y)]$.

C3. $E_{P(X, Y)}[f]=E_{P(X)P(Y)}[f]$, 其中 $E_{P(X, Y)}[f]$是 $f$关于 $P(X, Y)$的期望, $E_{P(X, Y)}[f]=\int f dP(x, y)$, 其中 $P(X, Y)$是 $(X, Y)$的联合分布函数, $P(X), P(Y)$分别是 $X, Y$的分布函数.

C4. $f(x, y)=f(x)f(y)$, 其中 $f(x), f(y)$分别是 $X$和 $Y$的密度函数, $f(x, y)$是 $X$和 $Y$的联合密度函数.

接下来, 我们逐一介绍基于这四个条件的相依性度量.

设 $d(\cdot, \cdot)$是一个距离度量, 比如: $d(x_{i}, x_{j})$是 $x_{i}$和 $x_{j}$之间的距离.我们可以把空间中的点按照参照点 $(x_{0}, y_{0})$分成四个部分, 如表 1所示.

特别地, 对于第 $i$和 $j$个观察值, 考虑以下两个二值随机变量: $I[d(x_{i}, X)\leq d(x_{i}, x_{j})]$和 $I[d(y_{i}, Y)\leq d(y_{i}, y_{j})]$, 其中 $I[.]$是示性函数, 我们得到一个新的 $2\times2$列联表, 并定义新的统计量.

令

$ \begin{equation} S(i, j)=\frac{(N-2)(A_{12}(i, j)A_{21}(i, j)-A_{11}(i, j)A_{22}(i, j))^{2}}{A_{., 1}(i, j)A_{., 2}(i, j)A_{1, .}(i, j)A_{2, .}(i, j)}, \end{equation} $

(3.1)

$S(i, j)$是Pearson卡方独立性检验统计量.

$ \begin{equation} T=\sum\limits_{i=1}^{N}\sum\limits_{j=1, j\neq i}^{N} S(i, j) . \end{equation} $

(3.2)

因为很难得到 $T$的渐渐分布, 因此要数值方法估计统计量 $T$的渐进分布, 而从计算 $p$值.因此HHG的计算量很大.

两个分布 $P$和 $Q$之间的最大平均散度(MMD)定义为

$ \begin{equation} MMD(P, Q)=\sup\limits_{f \in F}|E_{P}f-E_{Q}f|, \end{equation} $

(3.3)

其中 $F$为函数空间, 在计算时通常取再生核希尔伯特空间.

随机变量 $X$和 $Y$之间的协方差定义为

$ \begin{equation} cov(X, Y)=E_{X, Y}[XY]-E_{X}[X]E_{Y}[Y]. \end{equation} $

(3.4)

类似地, 函数空间中随机变量的协方差为

$ \begin{equation} \langle f, C_{XY}g\rangle _{F}=E_{X, Y}[f(X)g(Y)]-E_{X}[f(X)]E_{Y}[g(Y)], \end{equation} $

(3.5)

其中 $f \in F$, $g \in G$, $C_{XY}: G \rightarrow F $被称为协方差算子.

限制协方差定义为(COCO^[21])

$ \begin{equation} COCO(P_{X, Y};F, G)= \sup \{\langle f, C_{XY}g\rangle _{F}: f \in F, g \in G\} . \end{equation} $

(3.6)

因此我们需要给定函数空间 $F$和 $G$, 通常我们取再生核希尔伯特空间, 从而我们称COCO为核协方差.

Hilbert-Schmidt独立指数(HSIC^[10])定义为

$ \begin{equation} HSIC(P_{X, Y};F, G)=\parallel C_{XY} \parallel _{HS}^2, \end{equation} $

(3.7)

其中 $\parallel \cdot \parallel _{HS}^2$是Hilbert-Schmidt(HS)范数.

如果 $HSIC(P_{X, Y};F, G)=0$, $X$和 $Y$是统计独立的, 具体的细节请参见文献[10].

互信息^[27]

$ I(X, Y)=H(X)+H(Y)-H(X, Y)=H(X)-H(X|Y)=H(Y)-H(Y|X), $

其中 $H(X)$是香农熵.从定义我们可以看出 $I(X, Y)$取值于[0, + $\infty$), 并且 $I(X, Y)=0$当且仅当随机变量 $X$和 $Y$相互独立.为了满足性质(G), Linfoot^[18]建议取变换 $f(x)=\sqrt{1-{\rm e}^{-2x}}$, $I'(X, Y)=f(I(X, Y))$.

另外当 $X$和 $Y$之间有完全相关时, $I(X, Y)=H(Y)=H(X)$.这就意味着对于有不同联合分布却完全相关的 $(X, Y)$可以有不同的互信息.这就使得互信息具有不可比性.

对于高维正态分布情况, 我们有

$ \begin{equation} I(X;Y)=-\frac{1}{2}\log{\frac{|\Sigma|}{|\Sigma_{XX}||\Sigma_{YY}|}}=-\frac{1}{2}\log{\frac{|\Sigma|}{|B|}}, \end{equation} $

(3.8)

其中

$ \begin{equation} \Sigma=\left( \begin{array}{cc} \Sigma_{XX}~~&\Sigma_{XY} \\ \Sigma_{YX} ~~& \Sigma_{YY} \\ \end{array}\right) \end{equation} $

(3.9)

是 $X$和 $Y$的协方差矩阵,

$ \begin{equation} B=\left( \begin{array}{cc} \Sigma_{XX} ~&0 \\ 0~&\Sigma_{YY} \\ \end{array}\right). \end{equation} $

(3.10)

最小二乘互信息(Least-squares mutual information^[28])

$ \begin{equation} LSMI(X, Y)=\bigg(\int\!\!\!\int (k(x, y)-1)^{2}dF(x)dF(y)\bigg)^{1/2}, \end{equation} $

(3.11)

其中, $k(x, y)=\frac{f(x, y)}{f(x)f(y)}$是密度比.

最大信息系数^[1]可以用来检测各种各样的函数关系, 比如 $y=\sin(x)$, $y=\cos(x)$, $y=x^{2}$.当两个变量 $X$和 $Y$之间有严格的函数关系时, $MIC(X, Y)=1. $当两个变量 $X$和 $Y$是相互独立的, $MIC(X, Y)=0. $

$ \begin{equation} MIC = \frac{I(X, Y)}{H(X)}=1-\frac{H(X|Y)}{H(X)}, \end{equation} $

(3.12)

或

$ \begin{equation} MIC = \frac{I(X, Y)}{H(Y)}=1-\frac{H(Y|X)}{H(X)} . \end{equation} $

(3.13)

从定义可以看出, $0\leq MIC \leq 1$.另外, $H(Y|X)$是除去 $X$的信息后, $Y$中含有的信息量, 换句话说, $H(Y|X)$是 $Y$中含有的相对于 $X$的噪声信息.因此MIC可以理解成 $Y$中可以被 $X$解释的信息.相关学者指出用 $MIC-\rho^2$度量 $X$和 $Y$之间的非线性程度^[1].

对于任意的定义在 ${\Bbb R}^{p}\times {\Bbb R}^{q}$函数 $\gamma$, $L^2$中的加权范数 $\|.\|_{w}$定义为

$ \begin{equation} \|\gamma(t, s)\|_{w}^{2}=\int_{{\Bbb R}^{p+q}}|\gamma(t, s)|^{2}w(t, s){\rm d}t{\rm d}s, \end{equation} $

(3.14)

其中 $w(t, s)$是任意的权重函数.我们可以选择适当的权重函数 $w(t, s)$来定义新的相依性度量,

$ \begin{eqnarray*} \nu^{2}(X, Y, w) &=& \|f_{X, Y}(t, s)-f_{X}(t)f_{Y}(s)\|_{w}^{2} \nonumber\\ &=& \int_{{\Bbb R}^{p+q}}|f_{X, Y}(t, s)-f_{X}(t)f_{Y}(s)|^{2}w(t, s){\rm d}t{\rm d}s . \end{eqnarray*} $

我们有 $\nu^{2}(X, Y, w)=0$, 当且仅当 $X$和 $Y$是相互独立的.同时, 我们定义

$ \begin{eqnarray*} \nu^{2}(X, w) &=& \nu^{2}(X, X, w) \nonumber\\ &=& \|f_{X, X}(t, s)-f_{X}(t)f_{X}(s)\|_{w}^{2} \nonumber\\ &=& \int_{{\Bbb R}^{p+q}}|f_{X, X}(t, s)-f_{X}(t)f_{X}(s)|^{2}w(t, s){\rm d}t{\rm d}s, \end{eqnarray*} $

并类似地定义 $\nu^{2}(Y, w)$, 从而引出距离相关性(Dcor)的定义

$ \begin{equation} R^{2}(X, Y)=\frac{\nu^{2}(X, Y, w)}{\nu^{2}(X, w)\nu^{2}(Y, w)} \end{equation} $

(3.15)

当 $\nu^{2}(X, w)\nu^{2}(Y, w)>0$, 否则, $R^{2}(X, Y)=0. $更多的细节请参考文献[12].

在前面的章节中, 我们简单地给出了HHG, MMD, MIC, MI, LSMI, KMI, dCo, COCO, HISC的定义, 并分析了它们的性质.在这一小节中, 我们将首先讨论HSIC, COCO, KMI, MMD, HISCcop的核选择问题.然后再比较HHG, MIC, MI, LSMI, dcor的检验功效.

我们先讨论HSIC, COCO, MMD的核函数选择问题. HSIC是LSMI的逼近^[28], LSMI的功效要比HSIC高很多, 但是LSMI的计算相当耗时, 并不适用于大数据情形.其中HSIC, MMD, dcor是相互等价的^[29]. COCO只是一种特殊的MCC, 因此我们不讨论COCO的核选择问题, 只讨论HSIC的核函数选择问题.常见的核函数有: (1)径向基核函数(Radial Basis kernel, 也称为"Gaussian"核函数); (2)多项式核函数(Polynomial kernel); (3)线性核函数(Linear kernel); (4)超抛物线核函数(Hyperbolic tangent kernel); (5)拉普拉斯核函数(Laplacian kernel); (6) Bessel核函数; (7) ANOVA RBF核函数; (8)样条核函数(Spline kernel).

我们将用模拟算法A来比较基于不同的核函数的HSIC在8种不同的函数关系下的功效.这8种不同的函数关系为:二次函数(Quadratic), 三次函数(Cubic), Sine, Cosine, 圆(Circle), 阶梯函数(Step), Bell和双波叠加(Two Waves).具体的函数表达式在附录中给出, 模拟结果见图 4.

图 4中, 函数类型为: Quadratic, Cubic, Sine, Cosine, Circle, Step, Bell and Two Waves, 它们是按从左到右, 从上到下的顺序排列的.从图 4看出, 基于不同的核的HSIC的表现是很不一样的, 基于Laplacian核的HSIC的功效最高, 虽然它对于阶梯函数的功效并不是最优, 但是这也从另外一个侧面说明我们根本无法确定在什么情况下应该用什么样的核函数.这就大大限制了HSIC的使用.对于一致核函数(universal kernels^[30]), HSIC $(X, Y)=0$意味着 $X$和 $Y$是独立的.形如 $K(x, y)=g(\parallel x- y \parallel ^2)$的核函数是一致的^[30], 这就意味着径向基核函数, 拉普拉斯核函数是一致的, 其他的却不是.

我们比较HHG, MIC, MI, dcor在8种不同的函数关系下的功效.这8种不同的函数关系为:二次函数(Quadratic), 三次函数(Cubic), Sine, Cosine, 圆(Circle), 阶梯函数(Step), Bell和双波叠加(Two Waves).具体的函数表达式在附录中给出, 模拟结果见图 5, 本模拟中我们还是使用模拟算法A.为了方便比较, 我们同时也给出了相关系数的结果.

从图 5中, 可以看出: (1) MIC和dcor的表现并不令人满意.在上述8种函数关系中, 它们的表现没有优于其他办法的; (2) HHG在函数关系二次函数(Quadratic), 三次函数(Cubic), Sine, Cosine, 圆(Circle), Bell下的表现都比其他的方法要好; (3) MI在双波叠加的情况下有最高的功效; (4)对于阶梯函数, dcor和相关系数的表现非常相似.

给定数据集 $D=\{(x_i, y_i), i=1, 2, \cdots, n\}$, 学习理论的目标是找到一个函数 $f$, 满足 $Y=f(X)+\epsilon$.通常我们通过最小化残差平方和 $\Sigma_i(y_i-f(x_i))^2$来找到这样的 $f$.从而度量 $X$和 $Y$之间的相依性问题就转变成了是否存在 $f$使得 $Y=f(X)+\epsilon$.若不存在这样的函数, 则 $X$和 $Y$是相互独立的, 否则就是相互关联的.因此我们可以把相依性度量定义为

$ \begin{equation} DM(x, y)=\rho(y, \hat{f}(x)), \end{equation} $

(4.1)

其中 $\hat{f}$是 $f$的估计.

因此我们的关键是估计 $f$.通常有两种方法可以得到 $f$的估计: (1)用光滑化方法; 比如滑动平均, 加权平均, 小波去噪等. (2)在某个函数空间内估计 $f$.比如在线性函数空间下, 此时 $DM(x, y)$就是相关系数了.

估计 $f$之前, 必须先对 $f$做一些假设, 比如 $f$属于线性空间.因此相依性度量应该定义为:

$ \begin{equation} DM(x, y)=\sup \{\rho(y, f(x)), f \in H(x) \}, \end{equation} $

(4.2)

其中 $H(x)$是由 $x$张成的函数空间.

如果我们取 $H(x)$为线性空间, 即: $H(x)=\{ax+b: a, b \in {\Bbb R} \}$, 此时 $DM(x, y)$即为Pearson相关系数.

更一般地, 我们可以取 $H(x)$为可再生核希尔伯特空间, 多项式空间.取不同的空间, 可能会得到不同的相依性度量值.如果空间取得太小不能很好地逼近真实的 $f$, 就会导致错误.但是如果空间取得太大, 却会给计算带来困难.

随机相关系数是最近才提出来的用于度量两组变量间相关性的方法.下面的定理给出了一种 $f$的随机化估计方法.

定理4.1  (Rachimi-Brecht) 设 $p$是 $\Omega$上的分布函数, $|\phi(x;w)| \leq 1$.令

$ \begin{equation} F=\{f(x)=\int_{\Omega} \alpha(w) \phi(x;w) dw \mid |\alpha(w)| \leq Cp(w) \} \end{equation} $

(4.3)

从 $p$中抽取独立的样本 $w_1, w_2, \cdots, w_k$.给定 $0<\delta<0.5$, $c$是L-Lipschitz的损失函数, 数据集 $\{x_i, y_i\}_{i=1}^n$是来自 $P(X, Y)$的独立样本.设最小化损失 $c(f_k(x), y)$的参数为 $\alpha_{1}, \alpha_{2}, \cdots, \alpha_{k}$, 其中 $f_k(x)=\Sigma_{i=1}^k \alpha_i \phi(x;w_i)$, 则我们有

$ \begin{equation} E_p[c(f_k(x), y)]-\min\limits_{f \in F} E_p[c(f(x), y)] \leq \Theta \bigg\{\bigg(\frac{1}{\sqrt{n}} + \frac{1}{\sqrt{k}}\bigg) LC \sqrt{\log \frac{1}{\delta}}\bigg\} \end{equation} $

(4.4)

至少以概率 $1-2\delta$成立.

根据定理4.1, $f_k(x)=\Sigma_{i=1}^k \alpha_i \phi(x;w_i)$是 $f$的一个最优估计.

给定数据集 $X=(X_1, \cdots, X_d)'$= $(x_1, x_2, \cdots, x_n), $令

$ \begin{equation} \Phi(X;k, s)=\left( \begin{array}{ccc} \phi(w_1'x_1+b_1) &~ \cdots ~& \phi(w_k'x_1+b_k) \\ \vdots&\vdots&\vdots \\ \phi(w_1'x_n+b_1)&\cdots&\phi(w_k'x_n+b_k) \\ \end{array} \right)', \end{equation} $

(4.5)

其中, $w_i \sim N(0, sI), b_i \sim U[-\pi, \pi]$, $\phi(x;w)=\sin(w'x)$或者 $\cos(w'x)$.

随机相关系数^[17]

$ \begin{equation} RDC(X, Y;k, s)=\sup\limits_{\alpha, \beta} \rho(\alpha^T \Phi(X;k, s), \beta^T \Phi(Y;k, s)). \end{equation} $

(4.6)

在前面的章节, 我们回顾了一下现有的相依性度量方法, 并从三个角度概括了这些方法: (1)相关系数的推广; (2)基于独立性条件的相依性度量; (3)学习理论下的相依性度量.

表 3中给出了本文所总结的相依性度量的性质.我们可以看出并没有全美的相依性度量, CDC看上去是最好的相依性度量, 但面对高维数据却困难重重.

基于前面的讨论, 我们对现有的相依性度量方法有了很好的了解.从各章节后面的讨论中, 我们知道CDC, RDC, MI, HHG, 和基于拉普拉斯核的HSIC是各小类中最好的.下面我们将用模拟算法A分析它们在8种不同的函数关系下的检验功效.这8种函数关系分别为:二次函数(Quadratic), 三次函数(Cubic), Sine, Cosine, 圆(Circle), 阶梯函数(Step), Bell和双波叠加(Two Waves).具体的函数表达式在附录中给出, 模拟结果见图 6.

从结果中看出, HSIC的表现是最差的, 因为在这8种函数关系下, 它的检验功效一直都处于比较低的水平. MI则在双波情况下有比较明显的优势, 但在其他情况下也没有明显优势. RDC, CDC和HHG在这8中函数关系下的表现则非常相似.但是考虑到它们的性质, 我们认为CDC是最好的相依性度量.对于高维情形文献[4]已经给出了详细的比较, 这里就不再重复叙述了.

实际应用中, 研究人员都非常关心一个问题:用线性函数拟合数据够了吗?这触动了学者们提出非线性度量(non-linearity measurement, NLM). NLM具有以下性质: (N1) $NLM(X, Y)$ $\geq0$; (N2)若 $X$和 $Y$之间是线性函数关系, 那么 $NLM(X, Y)=0$; (N3)若 $X$和 $Y$之间是非线性函数关系, 那么 $NLM(X, Y) >0$, 且 $NLM(X, Y)$的值越大非线性的程度就越高.

学者们提出用 $MIC-\rho^2$作为一个非线性度量^[1], 即: $NLM_{MIC}(X, Y)=MIC(X, Y)-\rho(X, Y)^2$, 下标MIC为了凸显出非线性度量依赖于相依性度量MIC.显然 $NLM_{MIC}(X, Y)$并不满足性质(N1)^[1].类似地, 我们给出更一般的非线性度量

$ \begin{equation} NLM_{DM}(X, Y)=DM(X, Y)-|\rho(X, Y)|^{\alpha} . \end{equation} $

(5.1)

注意: (1)这里的相依性度量DM必须满足性质(C), 这样可以保证非线性度量的定义是有意义的; (2) $\alpha$是根据DM决定的, 对于MIC, $\alpha=2$, 然而对于CDC, $\alpha=1$; (3)非线性的判断条件: $DM(X, Y)>DM_0$, 且 $NLM_{DM}(X, Y)>NLM_0$, 这里 $DM_0$和 $NLM_0$是给定的阈值.

下面我们用简单的模拟来比较 $NLM_{MIC}$, $NLM_{CDC}$.

模拟算法B: (B0)令 $L=a:step:b$; (B1)从均匀分布中抽取500个点作为 $X$; (B2)计算 $Y=f(X)+L\times N(0, 1)$; (B3)计算NLM;

我们考虑两种不同的 $f$, (1) $f(x)=10x$; (2) $y=0.2x+(x-0.3)^2$.对于函数(1)我们取 $L=0:0.05:30$, 对函数(2)我们取 $L=0:0.01:10$.模拟结果见图 7.

为了保持统一, 我们取 $NLM_{CDC}= CDC^2-\rho^2$.从结果中我们可以看到, $NLM_{CDC}$要优于 $NLM_{MIC}$.

• 对于函数(1)理论上非线性度量值为0.但是从结果中看到, $ NLM_{CDC}\leq NLM_{MIC}$, 而且我们也可以发现在某些情况下 $NLM_{MIC}<0$.当噪声强度比较高时( $\rho(X, Y)^2\leq 0.4$), $NLM_{MIC}$明显要大于 $NLM_{CDC}$, 且 $NLM_{MIC}>0.2$, 因此很可能会给出错误的判断.

• 对于函数(2)理论上非线性度量要大于0.然而在噪声强度较弱时( $\rho(X, Y)^2>0.4$), $NLM_{MIC}$要明显小于 $NLM_{CDC}$.当噪声强度较强时却有相反的现象.因此当信号较强时, $NLM_{CDC}$要优于 $NLM_{MIC}$.

对于函数(1)理论上非线性度量值为0.但是从结果中看到, $ NLM_{CDC} \leq NLM_{MIC}$, 而且我们也可以发现在某些情况下 $NLM_{MIC}<0$.当噪声强度比较高( $\rho(X, Y)^2\leq 0.4$)的时候, $NLM_{MIC}$明显要大于 $NLM_{CDC}$, 且 $NLM_{MIC}>0.2$, 因此很可能会给出错误的判断.对于函数(2)理论上非线性度量要大于0.然而在噪声强度较弱时( $\rho(X, Y)^2>0.4$), $NLM_{MIC}$要明显小于 $NLM_{CDC}$.当噪声强度较强时却有想反的现象.因此当信号较强时, $NLM_{CDC}$要优于 $NLM_{MIC}$.

函数的定义

1. Linear: $y = x+\varepsilon$;

2. Quadratic: $y = 4(x-.5)^2+\varepsilon$;

3. Cubic: $y = 28(x-0.5)^3-48(x-.5)^2-12x+\varepsilon$;

4. Sine: $y = \sin(3\pi x)+\varepsilon$;

5. Cosine: $y = \cos(3\pi x)+\varepsilon$;

6. Power 0.25: $y = x^{1/4}+\varepsilon$;

7. Circle: $y = w+\sqrt{1 - (2x - 1)^2}+\varepsilon$;

8. Step $y = (x > 0.5)+\varepsilon$;

9. Bell $y = e^{-1/2x^2}+\varepsilon$;

10. Two Waves $y = (2w-1) \sin(3\pi x)+\varepsilon$;

11. Logarithm $y = \log(x)+\varepsilon$.

其中 $P(w=0)=P(w=1)=0.5$.