多元排列熵 Multivariate Permutation Entropy(含c++代码)

熵(Entropy)

信息论中熵的概念首次被香农提出，目的是寻找一种高效/无损地编码信息的方法：以编码后数据的平均长度来衡量高效性，平均长度越小越高效；同时还需满足“无损”的条件，即编码后不能有原始信息的丢失。这样，香农提出了熵的定义：无损编码事件信息的最小平均编码长度。

香农信息熵 Shannon entropy

香农信息熵是由香农提出的一个概念，它描述了信息源各可能事件发生的不确定性。这个概念在信息论中扮演着重要的角色，解决了对信息的量化度量问题。
一条信息的信息量大小和它的不确定性有直接的关系。比如说，我们要搞清楚一件非常非常不确定的事，或是我们一无所知的事情，就需要了解大量的信息。相反，如果我们对某件事已经有了较多的了解，我们不需要太多的信息就能把它搞清楚。所以，从这个角度，我们可以认为，信息量的度量就等于不确定性的多少。
信息量是对信息的度量，我们考虑一个离散的随机变量 x ，当我们观察到这个变量的具体值的时候，我们接收到了多少信息呢？
多少信息用信息量来衡量，我们接受到的信息量跟具体发生的事件有关。
信息的大小跟随机事件的概率有关。越小概率的事情发生了，产生的信息量越大，如湖南地震；越大概率的事情发生了产生的信息量越小，如太阳从东边升起来了（肯定发生， 没什么信息量）。
因此一个具体事件的信息量应该是随着其发生概率而递减的。

香农借鉴了热力学的概念，把信息中排除了冗余后的平均信息量称为“信息熵”，并给出了计算信息熵的数学表达式。

$H(x)=-\sum p(x_i)lo{g}_2(p(x_i)),i=1,2,..,n$。

其中，x表示信息，$x_i(i=1,2,..,n)$表示x的各种可能取值，$p(x_i)$表示x取值为$x_i$的概率，H的单位是比特。这个公式可以用来计算信息的不确定性，即信息熵。信息熵的提出解决了对信息的量化度量问题。
香农熵在生物信息领域基因表达分析中也有广泛的应用，如一些或一个基因在不同组织材料中表达情况己知，但如何确定这些基因是组织特异性表达，还是广泛表达的，那我们就来计算这些基因在N个样本中的香农熵，结果越趋近于log2(N)，则表明它是一个越广泛表达的基因，结果越趋近于0则表示它是一个特异表达的基因。

排列熵（Permutation Entropy，PE）

是用于衡量时间序列复杂程度的指标。排列熵是一种常用的非线性时序分析方法，用于对多元时间序列进行复杂性分析。
时间序列分析之排列熵（Permutation Entropy）
对于某个长度为n的排列x，其元素分别为$x_1,x_2,...,x_n$

1. 规定一个嵌入维度m（即m-neighborhood）和时间延迟t，进行相空间重构

2. 得到k个子序列，$k=n-(m-1)t$，每个子序列分别为：

(1) $x_1,x_{1+t},...,x_{1+(m-1)t}$

(2) $x_2,x_{2+t},...,x_{2+(m-1)t}$

(3) …

(4) $x_k,x_{k+t},...,x_{k+(m-1)t}$

3. 并把其转换为大小关系的排列（k个，共有m!种可能性）

4. 计算每种大小关系排列的概率P，P(排列)=该排列出现次数/k，

5. 计算这些概率的香农信息熵

按照步骤举个例子，便于理解：

x={2,4,5,6,3,7,1}，其长度n=7

1. 设嵌入维度m=3（3-neightborhood），时间延迟t=1（没有skip）

2. 得到k=n-(m-1)t=5个子序列，即：

(1) 2,4,5

(2) 4,5,6

(3) 5,6,3

(4) 6,3,7

(5) 3,7,1

3. 转换为大小关系的排列，分别为：
   有人说，“该数排第几”和“排在这个位置的是第几个数”这两种理解方法并不影响计算结果！
   原算法里的描述是“该数排第几”，但是这个例子里用的方法是“排在这个位置的是第几个数”。

(1) 1,2,3

(2) 1,2,3

(3) 2,3,1

(4) 2,1,3

(5) 2,3,1

4. 以上排列共有3种，分别为2次(1,2,3)，2次(2,3,1）和1次(2,1,3)，这些排列的概率分别为：

(1) P(1,2,3) = 2/5

(2) P(2,3,1) = 2/5

(3) P(2,1,3) = 1/5

5. 计算信息熵，得到$Hpe(3)=0.4\times lo{g}_{2}2.5+0.4\times lo{g}_{2}2.5+0.2lo{g}_{2}5=1.5219$

排列熵作为衡量时间序列复杂程度的指标，越规则的时间序列，它对应的排列熵越小；越复杂的时间序列，它对应的排列熵越大。但是这样的结果是建立在合适的 m的选择的基础上的，如果 m 的选取很小，如1或者2的话，那么它的排列空间就会很小（1!、2!）。经过研究表明，这个 m 的选取还是要根据实际情况来决定，一般而言，Bandt and Pompe 建议的取值是m = 3 , . . . , 7

多尺度排列熵多尺度熵（MultiScale Permutation Entropy）

关于多尺度排列熵的一些思考

多元排列熵（Multivariate Permutation Entropy,MPE或MvPE)

多元排列熵（Multivariate Permutation Entropy,MPE或MvPE)是排列熵的扩展，由于 EEG 每个通道的数据并非独立，这样的扩展非常必要。
考虑EEG通道的时间窗口大小为T秒，其采样频率为$f_s=\frac{1}{T}$；因此，每个窗口将包括$（f_{s}T）$个样本，即数据点。
对于每个通道$i\in [1,m]$，每个$h\in [1,n=d!]$（即对于每个“基序”），计数所有时间$s\in [1,f_{s}T-d]$，其中通道时间对$(i,s)$提供基序j。
将计数除以mT后获得的频率$p_{i,j}$是矩阵的项$P_t(m,n)＝p_{i,j}$，反映了基序在长度为T的时间片中的分布。
它保持${\sum }_{i=1}^m{\sum }_{j=1}^{d！}{p}_{i,j}=1$。
根据该程序，原始多元时间序列被转换为一个时间相关矩阵，相关统计数据和可以容易地提取熵。
特别地，计算边际相对值很容易描述基序分布的频率，如：
$p_j={\sum }_{i=1}^m{p}_{i,j},j=1,...,d!$，d表示多变量排列熵的跨通道复杂性可以是计算为$p_j$的排列熵：$H_{MPE}(s)=-{\sum }_{j=1}^{d!}{p}_j{\log }_2{p}_j$

通过相同的矩阵，也可以计算单通道多元排列熵，如下所示：

$H_E(i,s)=-{\sum }_{j=1}^{d!}m{p}_{i,j}{\log }_2(m{p}_{i,j}),i=1,2,...,m$

可计算出的一个有趣的量是多元排列熵和通过平均m个单通道排列熵得到的曲线之间的均方差。这个量被称为偶然性。当且仅当单通道分布重合时，它消失。如果它们是高度“相似”的，那么时间序列的整体复杂度有时是两项的总和：信道的平均复杂度和依赖于信道之间的不均匀性的休息。然而，对突发性对多元排列熵的影响的彻底分析超出了目前工作的范围。排列熵是一种对噪声（特别是高频噪声）稳健的测量方法：作为多元排列熵的平均操作实质上的轻微变化，它可以帮助吸收数据采集的一些不确定性。

多元排列熵在脑电图信号处理中很有用，因为如果它是在“遥远的”通道上计算的，即在不同的半球和/或不同的区域，它可以通过突出长期空间（非线性）相关性来提取跨通道的规律。

多元排列熵的源代码

百度

1.多元时序的permutation entropy python代码
2.基于有序模式的度量对多变量时间序列进行非线性分析（Matlab代码实现）
debug
1.python读取.mat文件，可使用scipy.io包内的loadmat()函数读取，由于读出的data是字典类型，需要对其处理，利用numpy进行矩阵转置存入X中，转为ndarray类型。

data = scio.loadmat("signal01.mat")
X=data['data'].T

用python读写.mat文件——使用scipy库的scipy.io和h5py库
2.ndarry

import numpy as np
array_a = np.array([[1, 2, 3], [4, 5, 6]])
print(array_a)
print('ndarray的维度: ', array_a.ndim)
print('ndarray的形状: ', array_a.shape)
print('ndarray的元素数量: ', array_a.size)
print('ndarray中的数据类型: ', array_a.dtype)
print(type(array_a))

3.ValueError: setting an array element with a sequence错误是因为你在尝试将一个序列赋值给一个数组元素，而数组的元素应该是单个的值而不是序列
4.matplotlib scatter散点图函数中参数c的使用
5.IndexError: too many indices for array
出错原因如：第一个数组的形状为（3,4），第二个数组的形状为（3，）。第二个数组缺少第二个维度。np.array无法使用此输入来构造矩阵（或长度类似的数组）。它只能创建一个列表数组。

github上的代码

MpePy库

在github上关于多元排列熵有两个python代码实现，对其分别进行验证。
一个是mpePy库，marisa mohr写的，链接，而且是老师给推荐的，里面函数有：

Calculation of different multivariate permutation entropy, e.g.,
计算不同的多元排列熵，
- pooled permutation entropy,
-合并排列熵，
- multivariate weighted permutation entropy,
-多元加权排列熵，
- multivariate multiscale permutation entropy,
-多元多尺度排列熵，
- multivariate ordinal pattern permutation entropy, and
-多元有序模式排列熵，以及
- multivariate permutation entropy based on principal component analysis
-基于主成分分析的多元排列熵

这就有点不会了，这里有如此多multivariate permutation entropy的实现，但都有定语，这可怎磨办，先下库看看。

pip install mpepy

下好了以后：

Requirement already satisfied: 
mpepy in c:\users\sylvia\appdata\local\packages\pythonsoftwarefoundation.python.3.10_qbz5n2kfra8p0\localcache\local-packages\python310\site-packages (0.0.2)
[notice] A new release of pip is available: 23.2.1 -> 23.3.1
[notice] To update, run: C:\Users\Sylvia\AppData\Local\Microsoft\WindowsApps\PythonSoftwareFoundation.Python.3.10_qbz5n2kfra8p0\python.exe -m pip install --upgrade pip

顺便更新一波pip

Requirement already satisfied: pip in c:\users\sylvia\appdata\local\packages\pythonsoftwarefoundation.python.3.10_qbz5n2kfra8p0\localcache\local-packages\python310\site-packages (23.2.1)
Collecting pip
  Obtaining dependency information for pip from https://files.pythonhosted.org/packages/47/6a/453160888fab7c6a432a6e25f8afe6256d0d9f2cbd25971021da6491d899/pip-23.3.1-py3-none-any.whl.metadata
  Downloading pip-23.3.1-py3-none-any.whl.metadata (3.5 kB)
Downloading pip-23.3.1-py3-none-any.whl (2.1 MB)
   ---------------------------------------- 2.1/2.1 MB 2.7 MB/s eta 0:00:00
Installing collected packages: pip
  Attempting uninstall: pip
    Found existing installation: pip 23.2.1
    Uninstalling pip-23.2.1:
      Successfully uninstalled pip-23.2.1
  WARNING: The scripts pip.exe, pip3.10.exe and pip3.exe are installed in 
  'C:\Users\Sylvia\AppData\Local\Packages\PythonSoftwareFoundation.Python.3.10_qbz5n2kfra8p0\LocalCache\local-packages\Python310\Scripts' which is not on PATH.
  Consider adding this directory to PATH or, if you prefer to suppress this warning, use --no-warn-script-location.
Successfully installed pip-23.3.1

发现下的pip的路径不在环境变量path里，加进去具体步骤。
1.还是报错Import "mpepy" could not be resolved
针对问题：python下好了第三方库 但vscode报错
解决方法如下，但我还是没解决，我猜应该是我之前移动了c盘的一些文件夹位置导致的。。Visual Studio Code(vscode)安装路径无法更改，如何处理
2.vscode报Import "numpy" could not be resolved，而我是有numpy库的，网上说引起问题的原因是安装了多个python版本，vscode无从选择，选好python interpreter就好了。
3.于是我选择不使用mpepy库，而是直接把实现这个库的函数整合到我的代码里，却发现她给的样例测试数据是pandas.DataFrame类型的，用ndarray不匹配，字典也不匹配。建了一个DataFrame类型的对象：Python Pandas中dataframe常用操作（创建、读取写入、切片等）

mfbm=pd.DataFrame({'data':[2,4,5,6,3,7,1]})

意料之中地继续报错：

\code1\multivariate_permutation_entropy.py:102: RuntimeWarning: divide by zero encountered in log
  ppe = ppe + (freq * np.log(freq))
\code1\multivariate_permutation_entropy.py:102: RuntimeWarning: invalid value encountered in double_scalars
  ppe = ppe + (freq * np.log(freq))
Traceback (most recent call last):
  File "multivariate_permutation_entropy.py", line 329, in <module>
    multivariate_weighted_permutation_entropy(mfbm, order = 3 , delay = 1)
  File "multivariate_permutation_entropy.py", line 201, in multivariate_weighted_permutation_entropy
    freq_var = sum_variance_per_perm / total_variance
ZeroDivisionError: division by zero

warning是存在数字太小的原因，取对数后结果趋于负无穷大。
error是不能除0。
问题这是库，我也不造咋改。
随手看了一下里面的代码，发现很多轮子。numpy笔记整理 multivariate_normal（多元正态分布采样）
还是看看第二个家伙的代码吧

Multivariate-multiscale-permutation-entropy

alberto-ara写的，一位来自mcgill麦吉尔大学的研究员，据说是加拿大的哈佛，看上去很靠谱。github链接
起手先pip install pyentrp，很快就下好了可以使用，这就很奇怪，不知道为什么之前的mpepy库一直配不好。
一共两个函数，第一个是def mpe(mts, m, d)

This function computes the multivariate permutation entropy of a multivariate time series.
该函数计算多变量时间序列的多变量排列熵。
The algorithm follows the formulations of Morabito et al. (2012) and it is based for the most part on
该算法遵循Morabito等人的公式。（2012），在很大程度上基于
Nikolay Donets' pyEntropy code (2012).
Nikolay Donets的pyEntropy代码（2012）。
INPUTS are:
输入是：
- mts: multivariate time series of shape "N series x N samples"
-mts：“N系列×N个样本”形状的多变量时间序列
- m: order of possible permutation motifs
-m：可能置换基序
- d: time-lag
-d：时间滞后
OUTPUTS are:
输出为：
- pe_channel: single series entropy
-pe_channel：单序列熵
- pe_cross: cross-series entropy
-pe_cross：交叉级数熵

很好这个看起来就是我们要的，不过我的测试数据是1000*1的向量，其实并没有多个样本。问题不大， 先拿小一点数据的去尝试。
1.vscode默认编码是utf-8，但cmd的默认编码方式是gbk，因此，如果不调vscode的编码，会导致终端输出乱码。解决方式就是在写代码之前调好编码方式，以免后续更多工作。
2. bool next_permutation(iterator start,iterator end)
当当前序列不存在下一个排列时，函数返回false，否则返回true。对数组num中的前n个元素进行全排列，并改变num数组的值。
3.memset初始化count数组时出现问题，手动实现初始化。

c++代码已开源至github：sylvia的github