在数据集上计算连续随机变量的信息熵和互信息--k-近邻估计方法

写在前面

• 信息熵(entropy)的原始定义是离散(discrete)的，后来发展了在连续域上的微分熵(differential entropy)。然而，通常在给定的数据集上，无法知道连续变量的概率分布，其概率密度函数也就无法获得，不能够用微分熵的计算公式。那要如何计算呢？
• 一种常见的方式是直方图法，它将连续变量的取值离散化，通过变量范围内划分bins，将不同的变量取值放入一个个bins中，然后统计其频率，继而使用离散信息熵的计算公式进行计算。然而，每个bin应该取多大的范围是很难确定的，通常需要反复计算获得最优的解。
• 一种无参熵估计法(non-parametric entropy estimation)可以避免划分bins来计算熵值，包括了核密度估计(kernel density estimator, KDE)和k-近邻估计(k-NN estimator)。
• 相比之下，直方图法不够精确，而核密度估计法运算量太大，k-近邻估计成为了普遍使用的一种计算连续随机变量的熵值方式。

一、计算离散随机变量的熵

计算离散随机变量的信息熵可以直接使用如下公式进行计算。
$$H(x)=-\sum _{i=1}^{N}p(x_i)log(p(x_i))$$
其中$N$是变量$x$在数据集中的样本数量，$p(x_i)$是$x_i$值在整个数据集中出现的频率。

二、k-近邻计算连续随机变量的熵

1.该方法最早由Kozachenko & Leonenko[1]在1987年提出的，基于该方法的一系列改进公式称为经典k-近邻(classical k-NN)。可以参考网站：Kozachenko-Leonenko estimator。在$D$维空间中的经典k-近邻估计计算公式[2]为：
$$H(x)\approx \psi (N)-\psi (k)+log(c_D)+\frac{D}{N}\sum _{i=1}^{N}log({\epsilon }_i)$$
其中，$\psi$是Digamma函数，${\epsilon }_i$是$x_i$到它第$k$个邻居的欧式距离(Euclidean distance)，$x_i$是变量$x$的一个采样点，$c_D=\frac{\pi \frac{D}{2}}{\Gamma (1+\frac{D}{2})}$，而$\Gamma$是Gamma函数。其他关于该公式的详细描述可以参考原文：https://homepages.tuni.fi/jaakko.peltonen/online-papers/lu2020aaai.pdf。

2.Gamma函数是阶乘函数$x!$在实数范围上的拓展，而Digamma函数是Gamma函数的对数的导数。其具体的定义式可以自行搜索。Gamma函数和Digamma函数虽然定义比较复杂，但是它们的特殊值可以很方便地使用递推公式来计算。

• Gamma函数递推计算如下：
  (1)$\Gamma (\frac{1}{2})=\sqrt{\pi }$,
  (2)$\Gamma (x+1)=x\Gamma (x)$.
• Digamma函数递推计算如下：
  (1)$\psi (1)=-\gamma$,
  (2)$\psi (x+1)=\psi (x)+\frac{1}{x}$.
  其中，$\gamma$是欧拉常数，近似值为0.57721 56649 01532 86060 65120。

三、计算两个离散随机变量的互信息

互信息是基于熵来计算的，用来表示两个随机变量的熵的重叠部分。当两个变量独立的时候，它们的互信息为0。两个离散随机变量的互信息计算公式为：
$$I(x;y)=\sum _{i,j}p(i,j)log\frac{p(i,j)}{p_x(i)p_y(j)}$$
其中，$p(i,j)$是$x=x_i$和$y=y_j$同时在整个数据集中出现的频率，$p_x(i)$是$x=x_i$在整个数据集中出现的频率，$p_y(j)$是$y=y_j$在整个数据集中出现的频率。

四、k-近邻计算两个连续随机变量的互信息

两个连续随机变量的互信息计算也是基于k-近邻估计，由Kraskov等人[3]在2004年提出。其中提出了两种基于k-近邻的计算的方法，此处仅给出第一种方法的公式：
$$I(x;y)\approx \psi (k)-<\psi (n_x+1)+\psi (n_y+1)>+\psi (N)$$
其中，$<...>$意为求$N$个点的均值。该公式的更详细描述可以参见[3]或者博客：机器学习特征选择：传统互信息、k-nearest neighbor互信息。

• 所提出的两种计算方法在大多数情况下所得到的结果是类似的，但一般来说，第一种方法有较小的统计误差和较大的系统误差。
• 而对于高维互信息量的计算来说，第二种方法会更加适合。
• 该计算方法在python的sklearn.feature_selection库的源码中使用，详细的实现可以参考其源码，调用方式为：

from sklearn.feature_selection import mutual_info_regression

• 两种计算互信息的方法均可以拓展到高维空间。这里给出第一种方法的高维拓展公式。假设有$m$个变量，则它们的互信息量为：
  $$I(x_1,x_2,...,x_m)\approx \psi (k)+(m-1)\psi (N)-<\psi (n_{x_1})+\psi (n_{x_2})+...+\psi (n_{x_m})>$$
  该公式的计算和二维的公式计算方式类似，关于计算的其他详细过程描述可以参考原文：https://journals.aps.org/pre/pdf/10.1103/PhysRevE.69.066138。

五、k-近邻计算一个连续随机变量和一个离散随机变量的互信息

该方法是基于k-近邻的两个连续随机变量互信息估计的拓展，由Ross[4]在2014年提出，该计算方法也在python的sklearn.feature_selection库的源码中使用，详细的实现可以参考其源码。其计算公式为：
$$I(x;y)\approx \psi (N)-<\psi (N_x)>+\psi (k)-<\psi (m)>$$
其中，$N_{x_i}$是和$x_i$有相同离散值的点数，$m_i$是在$N_{x_i}$个点中确定了第k个邻居后，在该邻居之内的点数。

一个计算例子如上图：

1. $x$是离散变量，$y$是连续变量，计算箭头指向的点的k近邻。
2. 在$y$变量中找它对应的x变量是红色取值的所有点（包括箭头指向点本身），点的个数即为$N_x$（此处$N_x=6$）。
3. 然后取离箭头指向点最近的第$k$个红色点（此处$k=3$）。
4. 以第k个红色点为界，取在以它为半径的范围中的所有点，统计点的个数（包括第k个红色点但不包括箭头指向点本身），即为$m$值（此处$m=6$）。

六、k的取值

一般来说，k越小，统计误差越大，系统误差越小，k越大则相反。通常k取3，这个也是sklearn.feature_selection库中默认的k值。

七、基于熵和互信息的其他信息度量计算

1. 联合熵：$$H(x,y)=H(x)+H(y)-I(x;y)$$
2. 条件熵：$$H(x|y)=H(x)-I(x;y)$$
3. 对称不确定性：$$SU(x,y)=\frac{2I(x;y)}{H(x)+H(y)}$$

但值得注意的是，如果x和y一个是在离散域上定义，一个是在连续域上定义，那么它们虽然可以用上面给的公式计算，但是由于它们的计算原理不同，可能效果会没有那么理想。但也不失为一种在信息领域的度量方式吧。

参考文献

[1]Kozachenko L F, Leonenko N N. Sample estimate of the entropy of a random vector[J]. Problemy Peredachi Informatsii, 1987, 23(2): 9-16.
[2]Lu C, Peltonen J. Enhancing Nearest Neighbor Based Entropy Estimator for High Dimensional Distributions via Bootstrapping Local Ellipsoid[C]. Proceedings of the AAAI Conference on Artificial Intelligence. 2020, 34(04): 5013-5020.
[3]Kraskov A, Stögbauer H, Grassberger P. Estimating mutual information[J]. Physical review E, 2004, 69(6): 066138.
[4]Ross B C. Mutual information between discrete and continuous data sets[J]. PloS one, 2014, 9(2): e87357.