CN109119166A - 一种基于矩阵传递熵的时间序列因果性分析方法及其计算机装置 - Google Patents
一种基于矩阵传递熵的时间序列因果性分析方法及其计算机装置 Download PDFInfo
- Publication number
- CN109119166A CN109119166A CN201810608100.0A CN201810608100A CN109119166A CN 109119166 A CN109119166 A CN 109119166A CN 201810608100 A CN201810608100 A CN 201810608100A CN 109119166 A CN109119166 A CN 109119166A
- Authority
- CN
- China
- Prior art keywords
- entropy
- causality
- matrix
- time series
- analysis method
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
Classifications
-
- G—PHYSICS
- G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
- G16H—HEALTHCARE INFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR THE HANDLING OR PROCESSING OF MEDICAL OR HEALTHCARE DATA
- G16H50/00—ICT specially adapted for medical diagnosis, medical simulation or medical data mining; ICT specially adapted for detecting, monitoring or modelling epidemics or pandemics
- G16H50/70—ICT specially adapted for medical diagnosis, medical simulation or medical data mining; ICT specially adapted for detecting, monitoring or modelling epidemics or pandemics for mining of medical data, e.g. analysing previous cases of other patients
Landscapes
- Engineering & Computer Science (AREA)
- Data Mining & Analysis (AREA)
- Health & Medical Sciences (AREA)
- Medical Informatics (AREA)
- Public Health (AREA)
- Biomedical Technology (AREA)
- Databases & Information Systems (AREA)
- Pathology (AREA)
- Epidemiology (AREA)
- General Health & Medical Sciences (AREA)
- Primary Health Care (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明公开了一种基于矩阵传递熵的时间序列因果性分析方法及其计算机装置和应用，通过Parzen窗变量的概率密度估计方法直接计算二阶Renyi熵，并将其推广到变量X的α阶Renyi熵，结合Gram矩阵定义矩阵传递熵，然后在传统的基于传递熵的因果性分析方法的基础上，通过矩阵传递熵计算因果性分析指标，最后按照因果性判断标准分析两个时间序列的因果性。与格兰杰因果性分析方法相比，该方法克服了“回归模型”的缺陷，能够用于具有非线性因果关系的时间序列中，与传统的基于传递熵的因果性分析方法相比，该方法鲁棒性好，易于计算而且计算复杂度低，适用于计算复杂度低的场合；将该分析方法应用于临床应用中用于探索患有呼吸暂停综合症病人的呼吸和心跳之间的因果性，具有良好的应用效果。
Description
技术领域
本发明属于时间序列分析领域；涉及一种基于矩阵传递熵的时间序列因果性分析方法；本发明还涉及能够实施上述分析方法的计算机装置。
背景技术
相关性分析描述两个变量之间的相互影响程度。而因果性分析揭示一个变量是如何对另一个变量进行作用的，即阐述变量之间的信息传递，具有方向性，在探索变量之间的关系、机器学习的特征提取和大脑效应网络的构建等领域中具有重要的应用价值。
目前，研究者们常采用格兰杰因果性分析方法探索变量之间的因果性。格兰杰因果性分析方法为2003年诺贝尔经济学奖得主Clive W.J.Granger所开创，用于分析经济变量之间的格兰杰因果关系。该方法基于线性回归模型，使用过去一些时刻点上所有信息的最佳最小二乘预测的方差来分析变量之间的因果性，易于理解并且计算复杂度低，能准确地分析具有线性因果关系且夹杂高斯噪声的时间序列之间的因果性。然而在实际应用中采集到的信号往往混有大量的非高斯噪声，而且变量之间的影响往往是非线性的，例如功能性核磁共振信号(fMRI)、脑电(EMG)和肌电(EEG)等信号。这使得格兰杰因果性分析方法的性能急剧下降。
2000年，Thomas Schreiber基于信息论中香农熵提出了新方法——传递熵用于分析两个时间序列之间的因果性，并应用于临床应用中探索患有呼吸暂停综合症病人的呼吸和心跳之间的因果性。变量Y到X的传递熵反映了Y的信息的加入对X的不确定性大小的改变，即Y传递给X信息量的大小，因此传递熵可以作为衡量因果性的指标。传递熵易于理解，仅考虑变量间的信息量传递，而不需要假定变量间具有特定形式的关系，并对非高斯噪声不敏感，因此具有比格兰杰因果性更好的适用性，尤其是对于夹杂非高斯噪声的非线性变量。然而在计算熵时，需要通过有限的数据估计变量的概率密度函数。通常变量的概率密度函数的估计是近似的，尤其是变量的联合概率密度函数和条件概率密度函数；同时概率密度函数的估计也存在计算复杂度高的问题。
发明内容
本发明提供了一种基于矩阵传递熵的时间序列因果性分析方法及其计算机装置；与格兰杰因果性分析方法相比，该方法克服了“回归模型”的缺陷，能够用于具有非线性因果关系的时间序列中，与传统的基于传递熵的因果性分析方法相比，该方法鲁棒性好，易于计算而且计算复杂度低，适用于计算复杂度低的场合。
本发明的技术方案是：一种基于矩阵传递熵的时间序列因果性分析方法，包括步骤S1，确定两个时间序列X和Y，并且生成关于X的向量自回归模型，并且采用贝叶斯信息准则确定X的嵌入维度；步骤S2，结合步骤S1中X的嵌入维度基于Silverman准则确定高斯核宽度；步骤S3，基于矩阵传递熵计算因果性分析指标，计算Y到X的因果性；其具体过程是：基于Parzen窗变量的概率密度估计方法计算变量X的二阶Renyi熵，并将其推广到变量X的α阶Renyi熵，结合Gram矩阵定义X的矩阵传递熵，得到Y到X的矩阵传递熵为：其中k为嵌入维度；计算时取α＝2；步骤S4，重复步骤S1-S3，计算得到X到Y的矩阵传递熵MTEX→Y；步骤S5，按照因果性判断标准分析并得到两个时间序列的因果性。
更进一步的，本发明的特点还在于：
其中步骤S3中X的二阶Renyi熵为：其中N为时间序列长度，Gσ(·)为高斯核函数，σ为高斯核宽度。
其中步骤S3中X的矩阵熵为Sα(X)＝-log[tr(Aα)]，其中A为变量X的Gram矩阵；X和Y的联合矩阵熵为其中B为Y的Gram矩阵，且变量X和Y的条件矩阵熵为Sα(X|Y)＝Sα(X，Y)-Sα(Y)；计算时取α＝2。
其中步骤S1中贝叶斯信息准则为：其中r为自回归模型误差向量，N为时间序列长度；所述嵌入维度k通过贝叶斯信息准则确定，其中P为预设的最大的嵌入维度。
其中步骤S2中的Silverman准则为：其中d为时间序列的个数，s1 2为第一个时间序列的方差，s2 2为第二个时间序列的方差。
其中步骤S5中因果性判断标准为：若MTEY→X＞＞0且MTEY→X/MTEX→Y＞1，则Y是X的因；若MTEX→Y＞＞0且MTEY→X/MTEX→Y＜1，则X是Y的因；若MTEY→X＞＞0，MTEX→Y＞＞0且MTEY→X/MTEX→Y≈1，则认为X和Y互为因果。
本发明还提供了一种计算机装置，该计算机装置能够实施上述的基于矩阵传递熵的时间序列因果性分析方法。
与现有技术相比，本发明的有益效果是：由于格兰杰因果性分析方法不能有效地分析夹杂着非高斯噪声且具有非线性关系的两个时间序列之间的因果性，而传递熵方法的计算复杂度高。本发明提出了基于矩阵传递熵的因果性分析方法，该方法具有良好的普适性，适合应用于非高斯噪声系统的因果性分析场合，能解决格兰杰因果性分析方法对叠加非高斯噪声且具有非线性因果关系的两个时间序列分析不准确和传递熵的计算复杂度高的问题，具有重要的研究意义和广泛的应用价值。
更进一步的，该方法通过Parzen窗变量的概率密度估计方法直接计算二阶Renyi熵，并将其推广到变量X的α阶Renyi熵，结合Gram矩阵定义矩阵传递熵，然后在传统的基于传递熵的因果性分析方法的基础上，通过矩阵传递熵计算因果性分析指标，最后按照因果性判断标准分析两个时间序列的因果性。与格兰杰因果性分析方法相比，该方法克服了“回归模型”的缺陷，可用于具有非线性因果关系的时间序列中；与传统的基于传递熵的因果性分析方法相比，该方法鲁棒性好，易于计算而且计算复杂度低，适用于要求方法的计算复杂度低的场合
附图说明
图1为本发明的流程示意图；
图2为本发明矩阵传递熵和现有传递熵的计算复杂度对比图；
图3为采用本发明矩阵传递熵在不同高斯核宽度下计算具有因果关系的时间序列的因果性分析指标图；
图4为采用现有传递熵在不同高斯核宽度下计算具有因果关系的时间序列的因果性分析指标图；
图5为采用本发明矩阵传递熵在不同高斯核宽度下计算具有非线性因果关系的时间序列的因果性分析指标图；
图6为采用现有递熵在不同高斯核宽度下计算具有非线性因果关系的时间序列的因果性分析指标图。
具体实施方式
下面结合附图和具体实施例对本发明的技术方案进一步说明。
本发明提供了一种基于矩阵传递熵的时间序列因果性分析方法，通过Parzen窗变量的概率密度估计方法直接计算二阶Renyi熵，并将其推广到变量X的α阶Renyi熵，结合Gram矩阵定义矩阵传递熵，然后在传统的基于传递熵的因果性分析方法的基础上，通过矩阵传递熵计算因果性分析指标，最后按照因果性判断标准分析两个时间序列的因果性。
如图1所示，该分析方法的具体过程包括以下步骤：
步骤S1，确定两个时间序列X和Y，并且生成关于X的向量自回归模型，同时采用贝叶斯信息准则确定X的嵌入维度；其中贝叶斯信息准则公式为：其中r为自回归模型误差向量；X的嵌入维度表示为：其中P为预设的最大的嵌入维度。
步骤S2，结合步骤S1中X的嵌入维度k，并且基于Silverman准则确定X的高斯核宽度；具体的，Silverman准则公式为：其中d为时间序列的个数，s1 2为第一个时间序列的方差，s2 2为第二个时间序列的方差，σ为核宽度。
步骤S3，基于矩阵传递熵计算因果性分析指标，计算Y到X的因果性；具体过程是，首先采用二阶Renyi熵计算变量X的熵，得到其中α＞0，α≠1，一般取α＝2；其中p(x)为概率密度函数；同时，Parzen窗变量的概率密度估计方法为其中Gσ(·)为高斯函数，σ为核宽度，N为样本量；从而得到基于Parzen窗变量的概率密度估计方法和二阶Renyi熵的X熵的计算公式为其中N为时间序列的长度。
然后结合Gram矩阵定义X的矩阵传递熵，设定Gram矩阵为A，其中从而得到X的二阶矩阵熵S2(X)＝-log[tr(A2)]，其中tr(·)为矩阵的迹函数；一般地，X的矩阵熵Sα(X)＝-log[tr(Aα)]；同理得到X和Y的联合矩阵熵为得到X和Y的条件矩阵熵为Sα(X|Y)＝Sα(X，Y)-Sα(Y)，其中B为变量Y的Gram矩阵。
基于传递熵的计算方法，得到Y到X的矩阵传递熵为计算时取α＝2。
步骤S4，重复上述步骤S1-S3，计算得到X到Y的矩阵传递熵MTEX→Y。
步骤S5，按照因果性判断标准分析并得到两个时间序列的因果性；具体的因果性判断标准为：若MTEY→X＞＞0且MTEY→X/MTEX→Y＞1，则Y是X的因；若MTEX→Y＞＞0且MTEY→X/MTEX→Y＜1，则X是Y的因；若MTEY→X＞＞0，MTEX→Y＞＞0且MTEY→X/MTEX→Y≈1，则认为X和Y互为因果。
如图2所示，在仿真实验下矩阵传递熵和传递熵计算复杂度的对比，该仿真实验中采用线性因果模型：其中N1服从均值为0方差为1的高斯分布，N2服从参数为[1.6，0，0.001，0]的Levy alpha-stable分布；从图2中可以看出矩阵传递熵的计算复杂度远远低于传递熵，因此将本发明的方法将矩阵传递熵引入到因果性分析领域。
考虑具有因果关系的两个时间序列和X是Y的因，其由上述线性因果模型生成，分别采用矩阵传递熵和传递熵计算时间序列的因果性分析指标分别如图3和图4所示，为在不同高斯核宽度下计算具有因果关系的时间序列的因果性分析指标图，其中图3和图4中的竖线是基于Silverman准则确定的高斯和宽度。从图3和图4中能够看出，基于矩阵传递熵和传递熵均能分析出时间序列的因果性，X是因，Y是果。但是矩阵传递熵的计算结果更平滑，传递熵存在抖动，说明传递熵相对于矩阵传递熵的不稳定性。
考虑具有非线性因果关系的两个时间序列和X是Y的因，其由模型生成，其中服从均值为0方差为1的高斯分布，服从参数为[1.6，0，0.001，0]的Levy alpha-stable分布，N＝1000为时间序列的长度，分别采用矩阵传递熵和传递熵计算时间序列的因果性分析指标，分别得到如图5和图6所示。
图5和图6分别采用矩阵传递熵和传递熵在不同高斯核宽度下计算具有非线性因果关系的时间序列的因果性分析指标图，竖线是基于Silverman准则确定的高斯核宽度。从图5和图6中能够看出，基于矩阵传递熵和传递熵均能分析出具有非线性因果关系的时间序列的因果性，X是因，Y是果。同样的，矩阵传递熵的计算结果更平滑，而传递熵存在抖动，说明了传递熵的不稳定性。
然而使用格兰杰因果性分析方法不能检测出具有非线性因果关系的时间序列。因此本发明基于矩阵传递熵的因果性分析方法的计算复杂度低且鲁棒性高，可用于具有非线性因果关系的时间序列分析中，具有极高的应用价值。
本发明还提供了一种计算机装置，该计算机装置能够实施上述的基于矩阵传递熵的时间序列因果性分析方法。
Claims (7)
1.一种基于矩阵传递熵的时间序列因果性分析方法，其特征在于，包括以下步骤：
步骤S1，确定两个时间序列X和Y，并且生成关于X的向量自回归模型，并且采用贝叶斯信息准则确定X的嵌入维度；
步骤S2，结合步骤S1中X的嵌入维度基于Silverman准则确定高斯核宽度；
步骤S3，基于矩阵传递熵计算因果性分析指标，计算Y到X的因果性；其具体过程是：基于Parzen窗变量的概率密度估计方法计算变量X的二阶Renyi熵，并将其推广到变量X的α阶Renyi熵，结合Gram矩阵定义X的矩阵传递熵，得到Y到X的矩阵传递熵为：
其中k为嵌入维度；计算时取α＝2；
步骤S4，重复步骤S1-S3，计算得到X到Y的矩阵传递熵MTEX→Y；
步骤S5，按照因果性判断标准分析并得到两个时间序列的因果性。
2.根据权利要求1所述的基于矩阵传递熵的时间序列因果性分析方法，其特征在于，所述步骤S3中X的二阶Renyi熵为：其中N为时间序列长度，Gσ(·)为高斯核函数，σ为高斯核宽度。
3.根据权利要求1和2任意一项所述的基于矩阵传递熵的时间序列因果性分析方法，其特征在于，所述步骤S3中X的矩阵熵为Sα(X)＝-log[tr(Aα)]，其中A为变量X的Gram矩阵；X和Y的联合矩阵熵为其中B为Y的Gram矩阵，且变量X和Y的条件矩阵熵为Sα(X|Y)＝Sα(X,Y)-Sα(Y)；计算时取α＝2。
4.根据权利要求2所述的基于矩阵传递熵的时间序列因果性分析方法，其特征在于，所述步骤S1中贝叶斯信息准则为：其中r为自回归模型误差向量，N为时间序列长度；所述嵌入维度k通过贝叶斯信息准则确定，其中P为预设的最大的嵌入维度。
5.根据权利要求2所述的基于矩阵传递熵的时间序列因果性分析方法，其特征在于，所述步骤S2中的Silverman准则为：其中d为时间序列的个数，s1 2为第一个时间序列的方差，s2 2为第二个时间序列的方差。
6.根据权利要求1所述的基于矩阵传递熵的时间序列因果性分析方法，其特征在于，所述步骤S5中因果性判断标准为：若MTEY→X＞＞0且MTEY→X/MTEX→Y＞1，则Y是X的因；若MTEX→Y＞＞0且MTEY→X/MTEX→Y＜1，则X是Y的因；若MTEY→X＞＞0，MTEX→Y＞＞0且MTEY→X/MTEX→Y≈1，则认为X和Y互为因果。
7.一种计算机装置，其特征在于，该计算机装置能够实施权利要求1所述的基于矩阵传递熵的时间序列因果性分析方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810608100.0A CN109119166A (zh) | 2018-06-13 | 2018-06-13 | 一种基于矩阵传递熵的时间序列因果性分析方法及其计算机装置 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810608100.0A CN109119166A (zh) | 2018-06-13 | 2018-06-13 | 一种基于矩阵传递熵的时间序列因果性分析方法及其计算机装置 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN109119166A true CN109119166A (zh) | 2019-01-01 |
Family
ID=64822192
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810608100.0A Pending CN109119166A (zh) | 2018-06-13 | 2018-06-13 | 一种基于矩阵传递熵的时间序列因果性分析方法及其计算机装置 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109119166A (zh) |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110766314A (zh) * | 2019-10-21 | 2020-02-07 | 中国民航信息网络股份有限公司 | 一种因果关系分析方法及装置 |
| CN111008363A (zh) * | 2019-11-21 | 2020-04-14 | 西安交通大学 | 多变量因果驱动的复杂机电系统服役安全态势评估方法 |
| CN113269336A (zh) * | 2021-07-19 | 2021-08-17 | 中国民用航空总局第二研究所 | 航班事件因果检测方法、装置、电子设备及可读存储介质 |
-
2018
- 2018-06-13 CN CN201810608100.0A patent/CN109119166A/zh active Pending
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110766314A (zh) * | 2019-10-21 | 2020-02-07 | 中国民航信息网络股份有限公司 | 一种因果关系分析方法及装置 |
| CN111008363A (zh) * | 2019-11-21 | 2020-04-14 | 西安交通大学 | 多变量因果驱动的复杂机电系统服役安全态势评估方法 |
| CN111008363B (zh) * | 2019-11-21 | 2021-11-19 | 西安交通大学 | 多变量因果驱动的复杂机电系统服役安全态势评估方法 |
| CN113269336A (zh) * | 2021-07-19 | 2021-08-17 | 中国民用航空总局第二研究所 | 航班事件因果检测方法、装置、电子设备及可读存储介质 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107481264B (zh) | 一种自适应尺度的视频目标跟踪方法 | |
| Kong et al. | L2RM: Low-rank linear regression models for high-dimensional matrix responses | |
| CN109002845A (zh) | 基于深度卷积神经网络的细粒度图像分类方法 | |
| Clayton | Some odds ratio statistics for the analysis of ordered categorical data | |
| CN109119166A (zh) | 一种基于矩阵传递熵的时间序列因果性分析方法及其计算机装置 | |
| CN109190537A (zh) | 一种基于掩码感知深度强化学习的多人物姿态估计方法 | |
| Qin et al. | Simulating and Predicting of Hydrological Time Series Based on TensorFlow Deep Learning. | |
| Yu et al. | Physformer++: Facial video-based physiological measurement with slowfast temporal difference transformer | |
| Mandic et al. | On the characterization of the deterministic/stochastic and linear/nonlinear nature of time series | |
| CN103226595B (zh) | 基于贝叶斯混合公共因子分析器的高维数据的聚类方法 | |
| CN111009321A (zh) | 一种机器学习分类模型在青少年孤独症辅助诊断中的应用方法 | |
| Delicado et al. | A small sample comparison of maximum likelihood, moments and L-moments methods for the asymmetric exponential power distribution | |
| CN109271876A (zh) | 基于时间演化建模和多示例学习的视频动作检测方法 | |
| CN109685334B (zh) | 一种新的基于多尺度理论的水文模型模拟评估方法 | |
| CN108959188A (zh) | 基于量化最小误差熵准则的格兰杰因果判辨方法 | |
| CN103034837A (zh) | 特征参数与脉象要素的关联 | |
| US9613123B2 (en) | Data stream processing | |
| Zhou et al. | Causality detection with matrix-based transfer entropy | |
| Yin et al. | Multiscale multifractal detrended cross-correlation analysis of traffic flow | |
| Paynabar et al. | Informative sensor and feature selection via hierarchical nonnegative garrote | |
| Bonett et al. | Confidence intervals for a ratio of binomial proportions based on paired data | |
| Luo et al. | Learning differential operators for interpretable time series modeling | |
| Ng et al. | Analysis of vector autoregressions in the presence of shifts in mean | |
| Hassani et al. | Does noise reduction matter for curve fitting in growth curve models? | |
| CN116975600A (zh) | 一种多尺度、多模态对比学习的生物医学时序信号应用算法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| RJ01 | Rejection of invention patent application after publication | ||
| RJ01 | Rejection of invention patent application after publication |
Application publication date: 20190101 |