贝叶斯网络

作者: Jacob Schreiber 联系方式: jmschreiber91@gmail.com

贝叶斯网络是一种通用的概率模型，它是pomegranate中展示的所有其他模型的超集。具体来说，贝叶斯网络是通过图结构对联合概率分布进行因子分解的一种方法，其中边的存在表示两个变量之间存在有向依赖关系，而边的缺失则表示条件独立性。重要的是，边的缺失本身并不表示两个变量完全独立，而是表示条件独立性。

%pylab inline
import seaborn; seaborn.set_style('whitegrid')
import torch

%load_ext watermark
%watermark -m -n -p torch,pomegranate

Populating the interactive namespace from numpy and matplotlib
torch      : 1.13.0
pomegranate: 1.0.0

Compiler    : GCC 11.2.0
OS          : Linux
Release     : 4.15.0-208-generic
Machine     : x86_64
Processor   : x86_64
CPU cores   : 8
Architecture: 64bit

初始化与拟合

与隐马尔可夫模型类似，贝叶斯网络由一组分布和连接它们的图结构组成。在这种情况下，该图仅是一系列有向无权重边。大多数贝叶斯网络要求该图是无环的。然而，由于pomegranate使用因子图进行推断，因此并不严格要求这一点。请参阅下面的推断部分。

同样地，与pomegranate中的其他模型类似，贝叶斯网络可以完全从数据中学习得到。然而，精确的结构学习是难以处理的，因此该领域已发展出多种近似方法。更多详情请参阅贝叶斯网络结构学习教程。

现在，让我们实现最简单的贝叶斯网络：一个子节点和一个父节点。

from pomegranate.distributions import Categorical
from pomegranate.distributions import ConditionalCategorical
from pomegranate.bayesian_network import BayesianNetwork

d1 = Categorical([[0.1, 0.9]])
d2 = ConditionalCategorical([[[0.4, 0.6], [0.3, 0.7]]])

model = BayesianNetwork([d1, d2], [(d1, d2)])

或者，可以使用add_distributions和add_edge方法以编程方式创建网络。

model2 = BayesianNetwork()
model2.add_distributions([d1, d2])
model2.add_edge(d1, d2)

一旦这些模型通过结构初始化后，就可以对数据进行拟合。

X = numpy.random.randint(2, size=(10, 2))
X

array([[0, 0],
       [0, 0],
       [0, 1],
       [0, 0],
       [0, 1],
       [0, 1],
       [0, 0],
       [1, 1],
       [1, 0],
       [0, 1]])

model.fit(X)

BayesianNetwork(
  (distributions): ModuleList(
    (0): Categorical()
    (1): ConditionalCategorical(
      (probs): ParameterList(  (0): Parameter containing: [torch.float32 of size 2x2])
      (_w_sum): [tensor([0., 0.])]
      (_xw_sum): [tensor([[0., 0.],
              [0., 0.]])]
      (_log_probs): [tensor([[-0.6931, -0.6931],
              [-0.6931, -0.6931]])]
    )
  )
  (_factor_graph): FactorGraph(
    (factors): ModuleList(
      (0): Categorical()
      (1): JointCategorical()
    )
    (marginals): ModuleList(
      (0): Categorical()
      (1): Categorical()
    )
  )
)

接下来，我们可以检查学习到的参数是否正确。

model.distributions[0].probs, X[:,0].mean()

(Parameter containing:
 tensor([[0.8000, 0.2000]]),
 0.2)

看起来模型学习到的边缘分布中1的概率等于均值，这是正确的。

我们还可以查看条件概率表，并与第一列为0时第二列为1的概率（通过取平均值）进行比较。

model.distributions[1].probs[0], (X[X[:,0] == 0][:,1]).mean()

(Parameter containing:
 tensor([[0.5000, 0.5000],
         [0.5000, 0.5000]]),
 0.5)

看起来也是正确的。

最后，如果我们知道想要学习的贝叶斯网络结构但不了解其参数表的参数，可以将该结构传入初始化函数并调用fit方法。

model3 = BayesianNetwork(structure=[(), (0,)])
model3.fit(X)

model3.distributions[1].probs[0]

Parameter containing:
tensor([[0.5000, 0.5000],
        [0.5000, 0.5000]])

概率与对数概率

与其他生成模型类似，贝叶斯网络可以根据分布和图结构计算示例的概率。由于贝叶斯网络只是联合概率表沿图结构的因式分解，因此示例的概率就是每个变量给定其父变量概率的乘积。

model.probability(X)

tensor([0.4000, 0.4000, 0.4000, 0.4000, 0.4000, 0.4000, 0.4000, 0.1000, 0.1000,
        0.4000])

model.log_probability(X)

tensor([-0.9163, -0.9163, -0.9163, -0.9163, -0.9163, -0.9163, -0.9163, -2.3026,
        -2.3026, -0.9163])

model.distributions[0].log_probability(X[:,:1]) + model.distributions[1].log_probability(X[:, :, None])

tensor([-0.9163, -0.9163, -0.9163, -0.9163, -0.9163, -0.9163, -0.9163, -2.3026,
        -2.3026, -0.9163])

与其他pomegranate模型类似，输入可以是列表、numpy数组或torch张量。

预测

学习型贝叶斯网络最实用的应用之一是对缺失值进行推理的能力。与传统预测问题不同（传统预测具有固定的输入集和一个或多个固定输出），贝叶斯网络推理会利用任何已知值的变量来推断未知值的变量。已知变量集在不同样本间可以变化，因此无需预先确定。

在pomegranate中，这是通过使用环状置信传播算法（有时也称为"和积"算法）实现的。该算法在后端构建的因子图上运行。与常规的联结树推理相比，这种方法的权衡在于：算法速度更快、更易于实现、对树状贝叶斯网络具有精确性，甚至能为循环网络提供估计值；但在其他情况下不能保证推理的精确性，当网络存在循环时甚至可能无法收敛。

pomegranate中预测方法的实现与其他模型略有不同。首先，未观测变量使用torch.masked.MaskedTensor对象表示，该对象包含基础数据和掩码，其中True表示值已观测，False表示未观测。当掩码为False时，基础值是什么并不重要。

X_torch = torch.tensor(X[:4])
mask = torch.tensor([[True, False],
                     [False, True],
                     [True, True],
                     [False, False]])

X_masked = torch.masked.MaskedTensor(X_torch, mask=mask)

/home/jmschr/anaconda3/lib/python3.9/site-packages/torch/masked/maskedtensor/core.py:156: UserWarning: The PyTorch API of MaskedTensors is in prototype stage and will change in the near future. Please open a Github issue for features requests and see our documentation on the torch.masked module for further information about the project.
  warnings.warn(("The PyTorch API of MaskedTensors is in prototype stage "

如果您已经在张量中将某些值设置为缺失值（例如-1），您可以简单地执行torch.masked.MaskedTensor(X, mask=(X != -1))。

创建MaskedTensor后，您可以将其传入贝叶斯网络的predict、predict_proba或predict_log_proba方法。需要注意的是，当提供数据时，这些值将以概率1的形式出现在输出中。概率分布仅针对未观测变量提供。

model.predict_proba(X_masked)

[tensor([[1.0000, 0.0000],
         [0.8000, 0.2000],
         [1.0000, 0.0000],
         [0.8000, 0.2000]]),
 tensor([[0.5000, 0.5000],
         [1.0000, 0.0000],
         [0.0000, 1.0000],
         [0.5000, 0.5000]])]

您可能会注意到，这些函数的输出形状与其他方法不同。由于无法保证所有变量都具有相同数量的类别，pomegranate无法返回一个维度为类别数量的单一张量。相反，pomegranate选择返回一个张量列表，其中列表中的每个元素代表一个变量，该张量的维度为(n_examples, n_categories)，表示该维度的类别数量。理论上，可以返回一个大小为(n_examples, n_dimensions, max_n_categories)的单一张量，其中max_n_categories是所有维度中类别的最大数量，但用户可能仍会选择切掉不必要的类别，而且无法保证不会出现具有大量类别的单个变量，从而大幅增加所需内存量。

model.predict_log_proba(X_masked)

[tensor([[ 0.0000,    -inf],
         [-0.2231, -1.6094],
         [ 0.0000,    -inf],
         [-0.2231, -1.6094]]),
 tensor([[-0.6931, -0.6931],
         [ 0.0000,    -inf],
         [   -inf,  0.0000],
         [-0.6931, -0.6931]])]

与前两种方法不同，predict方法会保留原始数据的形状，但会根据掩码将缺失值替换为predict_proba方法计算出的最可能值。

model.predict(X_masked)

tensor([[0, 0],
        [0, 0],
        [0, 1],
        [0, 0]])

摘要

贝叶斯网络的摘要统计与其他模型的摘要统计工作方式完全相同。当给定完整数据时，摘要统计信息是通过最大似然估计(MLE)得出的，这些统计信息可以跨批次相加。这意味着贝叶斯网络可以通过核外计算的方式进行学习。

重要的是，贝叶斯网络必须已经具备结构，否则它将不知道每批次需要计算哪些统计量。这是因为结构学习尚未实现外存计算方式，否则必须在第一批数据上完成。如果您希望这样做，则应在第一批数据上使用fit方法，并在后续批次上使用summarize方法。

X = torch.randint(2, size=(20, 4))

d1 = Categorical([[0.1, 0.9]])
d2 = Categorical([[0.3, 0.7]])
d3 = ConditionalCategorical([[[0.4, 0.6], [0.3, 0.7]]])
d4 = ConditionalCategorical([[[[0.8, 0.2], [0.1, 0.9]], [[0.1, 0.9], [0.5, 0.5]]]])

model = BayesianNetwork([d1, d2, d3, d4], [(d1, d3), (d2, d4), (d3, d4)])
model.summarize(X[:5])
model.summarize(X[5:])
model.from_summaries()

拟合完成后，我们可以查看学习到的参数是什么。

model.distributions[0].probs, model.distributions[1].probs

(Parameter containing:
 tensor([[0.6000, 0.4000]]),
 Parameter containing:
 tensor([[0.5500, 0.4500]]))

我们可以将它们与数据中的实际平均值进行比较，这应该与前两个维度的学习参数相同。

torch.mean(X.type(torch.float32), dim=0)

tensor([0.4000, 0.4500, 0.4500, 0.6000])

采样

X_sample = model.sample(1000000).type(torch.float32)

model._parents

[(), (), (0,), (1, 2)]

model.distributions[0].probs

Parameter containing:
tensor([[0.6000, 0.4000]])

model.distributions[1].probs

Parameter containing:
tensor([[0.5500, 0.4500]])

X_sample[:, :2].mean(dim=0)

tensor([0.4006, 0.4501])

model.distributions[2].probs[0]

Parameter containing:
tensor([[0.4167, 0.5833],
        [0.7500, 0.2500]])

X_sample[X_sample[:, 0] == 0, 2].mean(), X_sample[X_sample[:, 0] == 1, 2].mean()

(tensor(0.5836), tensor(0.2500))

model.distributions[3].probs[0]

Parameter containing:
tensor([[[0.5000, 0.5000],
         [0.2000, 0.8000]],

        [[0.4000, 0.6000],
         [0.5000, 0.5000]]])

X_sample[(X_sample[:, 1] == 1) & (X_sample[:, 2] == 0), 3].mean()

tensor(0.5998)