估计置信区间#

我们训练的生成式因果模型和特定的样本集可能无法准确表示真实情况，并且我们的算法中通常存在数值近似。当我们在不计算置信区间的情况下回答因果查询时，我们实际上获得了点估计。然而，在评估我们结果的置信度时，这些点估计的用处不大。置信区间帮助我们量化由于建模不准确和数值近似引入的结果中的不确定性。

因此，我们建议将置信区间作为所有因果查询的默认选项。这可以避免由于数据倾斜或模型拟合不佳而对特定值产生错误置信。

为了计算置信区间，我们使用一种称为自举法（bootstrapping）的方法，该方法已在函数dowhy.gcm.confidence_intervals()中实现。

如何使用#

假设我们有一个函数，它通过蒙特卡洛方法计算圆周率 Pi，并进行 1000 次试验

>>> import numpy as np
>>>
>>> def compute_pi_monte_carlo():
>>>     trials = 1000
>>>     return 4*(np.random.default_rng().uniform(-1, 1, (trials,))**2+
>>>               np.random.default_rng().uniform(-1, 1, (trials,))**2 <= 1).sum() / trials

执行此函数将得到一个结果

>>> compute_pi_monte_carlo()
3.188

这看起来是一个合理的值，但考虑到它是使用随机方法计算的，我们无法判断该结果有多可靠。现在，让我们使用dowhy.gcm.confidence_intervals()来计算置信区间

>>> from dowhy import gcm
>>>
>>> median, intervals = gcm.confidence_intervals(compute_pi_monte_carlo,
>>>                                              num_bootstrap_resamples=1000)
>>> median, intervals
(array([3.136]), array([[3.0518, 3.224 ]]))

结果中的第一项是所有计算结果的中位数。第二项是我们有 95% 的信心计算值将落入的区间。通过增加trials的数量，我们可以提高compute_pi_monte_carlo的可靠性，从而缩小置信区间。例如，当trials = 10_000时，我们得到一个区间 [3.1132 , 3.16602]；当trials = 100_000时，我们得到 [3.13256 , 3.150016]。

计算因果查询的置信区间#

有了这个理解，我们现在可以将confidence_intervals应用于我们的一个因果查询，例如distribution_change。首先生成一些数据

>>> import pandas as pd
>>> from scipy.stats import halfnorm
>>>
>>> X = halfnorm.rvs(size=1000, loc=0.5, scale=0.2)
>>> Y = halfnorm.rvs(size=1000, loc=1.0, scale=0.2)
>>> Z = np.maximum(X, Y) + np.random.normal(loc=0, scale=1, size=1000)
>>> data_old = pd.DataFrame(data=dict(X=X, Y=Y, Z=Z))
>>>
>>> X = halfnorm.rvs(size=1000, loc=0.5, scale=0.2)
>>> Y = halfnorm.rvs(size=1000, loc=1.0, scale=0.2)
>>> Z = X + Y + np.random.normal(loc=0, scale=1, size=1000)
>>> data_new = pd.DataFrame(data=dict(X=X, Y=Y, Z=Z))

这里，Z 的因果机制在旧数据和新数据之间发生变化。和往常一样，接下来我们将因果关系建模为 SCM 并分配因果机制

>>> import networkx as nx
>>> causal_model = gcm.StructuralCausalModel(nx.DiGraph([('X', 'Z'), ('Y', 'Z')]))  # X -> Z <- Y
>>> gcm.auto.assign_causal_mechanisms(causal_model, data_old)

现在，不再直接调用distribution_change，而是使用confidence_intervals调用它。然而，由于confidence_intervals接受一个不带参数的函数，我们必须首先定义一个不带参数的函数

>>> def f():
>>>     return gcm.distribution_change(causal_model, data_old, data_new, 'Z')

然后

>>> gcm.confidence_intervals(f)
({'X': 0.0005804787010800414, 'Y': 0.0030143299612704804, 'Z': 0.1724206687905231},
{'X': array([-0.00427554,  0.00891714]),
'Y': array([-0.00605089,  0.01782684]),
'Z': array([0.15441771, 0.19062329])})

结果再次包含第一项，即贡献分值的中位数。第二项包含每个变量的贡献分值的区间。

为了避免定义新函数，我们也可以通过使用 lambda 来简化该调用

>>> gcm.confidence_intervals(lambda: gcm.distribution_change(causal_model,
>>>                                                          data_old, data_new,
>>>                                                          target_node='Z'))

便捷地在训练数据的随机子集上进行图训练自举#

GCM 软件包中的许多因果查询都需要一个训练好的因果图作为第一个参数。为了计算这些方法的置信区间，我们需要使用不同的随机数据子集多次显式地重新训练我们的因果图，并使用每个新训练的图运行因果查询。为了方便地做到这一点，GCM 软件包提供了一个函数fit_and_compute。假设我们有data和一个因果图

>>> Z = np.random.normal(loc=0, scale=1, size=1000)
>>> X = 2*Z + np.random.normal(loc=0, scale=1, size=1000)
>>> Y = 3*X + 4*Z + np.random.normal(loc=0, scale=1, size=1000)
>>> data = pd.DataFrame(dict(X=X, Y=Y, Z=Z))
>>>
>>> causal_model = gcm.StructuralCausalModel(nx.DiGraph([('Z', 'Y'), ('Z', 'X'), ('X', 'Y')]))
>>> gcm.auto.assign_causal_mechanisms(causal_model, data)

我们现在可以使用fit_and_compute，如下所示

>>> strength_median, strength_intervals = gcm.confidence_intervals(
>>>     gcm.fit_and_compute(gcm.arrow_strength,
>>>                         causal_model,
>>>                         bootstrap_training_data=data,
>>>                         target_node='Y'))
>>> strength_median, strength_intervals
({('X', 'Y'): 45.90886398636573, ('Z', 'Y'): 15.47129383737619},
{('X', 'Y'): array([42.88319632, 50.43890079]), ('Z', 'Y'): array([13.44202416, 17.74266107])})

运行时成本与置信度#

在某些场景下，多次重新训练因果图的成本非常高。例如，当使用 AutoGluon 作为加性噪声模型的预测模型时，单个fit执行可能很快需要数小时。对 20 个重采样进行自举则会很快需要数天，具体取决于并行化程度。

因此，有时需要权衡的是只在因果查询上进行自举，而不是在训练上。为了使其与我们上面使用的方法类似，有一个函数bootstrap_sampling()。此函数假定因果图已经过训练。然后可以按如下方式使用它

>>> gcm.fit(causal_model, data)
>>>
>>> strength_median, strength_intervals = gcm.confidence_intervals(
>>>     gcm.bootstrap_sampling(gcm.arrow_strength,
>>>                            causal_model,
>>>                            target_node='Y'))
>>>
>>> strength_median, strength_intervals
({('X', 'Y'): 46.07299374572871, ('Z', 'Y'): 15.358850280972195},
{('X', 'Y'): array([44.95914495, 47.63918151]), ('Z', 'Y'): array([15.04323069, 15.72570547])})

dowhy.gcm.bootstrap_sampling()实现了与我们之前看到的 lambda 相同的功能。然而，使用bootstrap_sampling更具表现力。因此，对于所有将因果图作为第一个参数的因果查询，我们建议使用它而不是 lambda。

理解置信区间的必要性#

正如开头所解释的，模型和特定的样本集可能无法准确表示真实情况，并且我们的算法中通常存在数值近似。误差主要来自三个来源

  • 模型“最优”参数的方差，即使用相同/略有不同的数据运行fit两次可以得到两个不同的模型（当预测模型的拟合过程中存在随机因素时，情况会更糟）。例如，即使使用完全相同的数据训练 AutoGluon 模型两次，也会返回两个性能略有不同的模型。

  • 算法中的近似。例如，当估计具有 6 个或更多上游节点的分布变化时，我们将只近似 Shapley 值（该近似具有随机因素）。因此，使用完全相同的输入和生成的样本运行分布变化两次将返回两个不同的结果。

  • 即使我们不近似 Shapley 值，算法通常也基于生成模型中的一些特定样本，即方差也可能来自两组抽取样本之间的方差。因此，即使算法本身是确定性的（例如评估 Shapley 值的​​所有可能子集），我们也会使用随机生成的样本，这在两次调用同一算法时会产生不同的结果。