在因果发现中的一个常见任务是证明和验证为什么一个假设或重建的因果网络是合理的。在这里，我们展示如何利用估计的图来构建一个模型，解释数据集的滞后相关结构。具体步骤如下：

1、估计（马尔可夫等价类的）因果图
2、如果马尔可夫等价类有多个成员（存在未定向的边），选择类的一个成员，这可以自动完成
3、对从图中提取的因果父节点进行线性结构因果模型拟合
4、估计残差的噪声协方差矩阵
5、利用这种噪声结构构建线性高斯结构因果模型，并使用与数据相同的样本大小生成许多实现
6、展示原始数据的滞后相关性与生成数据的（整体均值和置信区间）滞后相关性

# Imports
import numpy as np
import matplotlib
from matplotlib import pyplot as plt
%matplotlib inline     

import tigramite
from tigramite import data_processing as pp
from tigramite.toymodels import structural_causal_processes as toys

from tigramite.toymodels import surrogate_generator 

from tigramite import plotting as tp
from tigramite.pcmci import PCMCI
from tigramite.independence_tests.parcorr import ParCorr
from tigramite.models import Models, Prediction

import math
import sklearn
from sklearn.linear_model import LinearRegression

0、生成一些具有同时依赖关系的示例过程

#正常显示坐标轴负号
matplotlib.rcParams['axes.unicode_minus'] = False
np.random.seed(14)     # Fix random seed
lin_f = lambda x: x
links_coeffs = {0: [((0, -1), 0.7, lin_f)],
                1: [((1, -1), 0.8, lin_f), ((0, -1), 0.3, lin_f)],
                2: [((2, -1), 0.5, lin_f), ((0, -2), -0.5, lin_f)],
                3: [((3, -1), 0., lin_f)], #, ((4, -1), 0.4, lin_f)],
                4: [((4, -1), 0., lin_f), ((3, 0), 0.5, lin_f)], #, ((3, -1), 0.3, lin_f)],
                }
T = 200     # time series length
# Make some noise with different variance, alternatively just noises=None
noises = np.array([(1. + 0.2*float(j))*np.random.randn((T + int(math.floor(0.2*T)))) 
                   for j in range(len(links_coeffs))]).T

data, _ = toys.structural_causal_process(links_coeffs, T=T, noises=noises, seed=14)
T, N = data.shape

# For generality, we include some masking
# mask = np.zeros(data.shape, dtype='int')
# mask[:int(T/2)] = True
mask=None

# Initialize dataframe object, specify time axis and variable names
var_names = [r'$X^0$', r'$X^1$', r'$X^2$', r'$X^3$', r'$X^4$']
dataframe = pp.DataFrame(data, 
                         mask=mask,
                         datatime = {0:np.arange(len(data))}, 
                         var_names=var_names)
tp.plot_timeseries(dataframe=dataframe); plt.show()

现在让我们看一下滞后相关性矩阵。我们希望找到一个解释它的结构性因果模型。

tau_max = 10
parcorr = ParCorr(significance='analytic', 
#                   mask_type='y'
                 )
pcmci = PCMCI(
    dataframe=dataframe, 
    cond_ind_test=parcorr,
    verbosity=0)
original_correlations = pcmci.get_lagged_dependencies(tau_max=tau_max, val_only=True)['val_matrix']
lag_func_matrix = tp.plot_lagfuncs(val_matrix=original_correlations, setup_args={'var_names':var_names, 
                                    'x_base':5, 'y_base':.5}); plt.show()

1、估计（马尔可夫等价类的）因果图

在这里，我们使用一些合理的 tau_max 和 pc_alpha 运行 PCMCI+。当然，估计的模型会略微依赖于这些参数。

results = pcmci.run_pcmciplus(tau_max=tau_max, pc_alpha=0.01)
tp.plot_graph(results['graph']); plt.show()

2、如果马尔可夫等价类有多个成员（存在未定向的边），选择类的一个成员，这可以自动完成

PCMCI+ 的结果是一个 CPDAG（存在未定向的边）。我们使用一些内部函数选择此类别的一个 DAG 成员。

# First create order that is based on some feature of the variables
# to avoid order-dependence of DAG, i.e., it should not matter
# in which order the variables appear in dataframe
# Here we use the sum of absolute val_matrix values incident at j
val_matrix = results['val_matrix']
variable_order = np.argsort(
                    np.abs(val_matrix).sum(axis=(0,2)))[::-1]

# Transform conflicting links to unoriented links as a hack, might not work...
graph = results['graph']
graph[graph=='x-x'] = 'o-o'
dag = pcmci._get_dag_from_cpdag(
                cpdag_graph=graph,
                variable_order=variable_order)

tp.plot_graph(dag); plt.show()

以下步骤都包含在 toymodels.surrogate_generator 模块的 generate_linear_model_from_data 函数中。

3、对从图中提取的因果父节点进行线性结构因果模型拟合

使用调用 Models() 的 Prediction 类。

4、估计残差的噪声协方差矩阵

使用 Prediction 类。

5、利用这种噪声结构构建线性高斯结构因果模型

parents = toys.dag_to_links(dag)
print(parents)

{0: [(0, -1)], 1: [(0, -1), (1, -1)], 2: [(0, -2), (2, -1)], 3: [], 4: [(3, 0)]}

现在使用与数据相同的样本大小生成实际数据。

realizations = 100

generator = surrogate_generator.generate_linear_model_from_data(dataframe, parents, tau_max, realizations=realizations, 
                generate_noise_from='covariance')
datasets = {}
for r in range(realizations):
    datasets[r] = next(generator)

6、展示原始数据的滞后相关性与生成数据的（集合均值和置信区间）滞后相关性。

correlations = np.zeros((realizations, N, N, tau_max + 1))
for r in range(realizations):
    pcmci = PCMCI(
        dataframe=pp.DataFrame(datasets[r]), 
        cond_ind_test=ParCorr(),
        verbosity=0)
    correlations[r] = pcmci.get_lagged_dependencies(tau_max=tau_max, val_only=True)['val_matrix']
# Get mean and 5th and 95th quantile
correlation_mean = correlations.mean(axis=0)
correlation_interval = np.percentile(correlations, q=[5, 95], axis=0)
# Plot lag functions of mean and 5th and 95th quantile together with original correlation in one plot
lag_func_matrix = tp.setup_matrix(N=N, tau_max=tau_max, x_base=5, figsize=(10, 10), var_names=var_names)
lag_func_matrix.add_lagfuncs(val_matrix=correlation_mean, color='black')
lag_func_matrix.add_lagfuncs(val_matrix=correlation_interval[0], color='grey')
lag_func_matrix.add_lagfuncs(val_matrix=correlation_interval[1], color='grey')
lag_func_matrix.add_lagfuncs(val_matrix=original_correlations, color='red')
lag_func_matrix.savefig(name=None)

现在，原始数据中的滞后函数（红色）应该大致在生成数据的相关性范围（灰色）的90%之内。这意味着，具有与重建的因果图相同链接（以及估计系数和噪声结构）的线性高斯模型能够很好地解释原始数据的整体滞后相关性结构。因此，这是一个简洁的因果模型，可以解释数据。

最后，我们还可以展示（混合的）相关性图。

pcmci = PCMCI(
    dataframe=dataframe, 
    cond_ind_test=parcorr,
    verbosity=0)
original_correlations_pvals = pcmci.get_lagged_dependencies(tau_max=tau_max)['p_matrix'] 
tp.plot_graph(graph=original_correlations_pvals<=0.01); plt.show()