使用pgmpy构建贝叶斯网络并且进行参数学习

pgmpy - Python中的概率图模型库

PGMPY是Python中的一个概率图模型库，用于创建、操作和推断概率图模型。它提供了一组工具，用于参数学习和结构学习，以帮助我们理解变量之间的依赖关系。以下是一些主要概念：

  • 参数学习：当我们有一组数据样本和一个有向无环图（DAG）来表示变量之间的依赖关系时，参数学习旨在估计单个变量的（条件）概率分布。

  • 结构学习：当我们有一组数据样本时，结构学习旨在估计一个DAG，该DAG 捕获变量之间的依赖关系。

注意：在pgmpy的最新版本0.1.23中，BayesianModel 已经更名为 BayesianNetwork。

你可以在官方链接中找到更多信息和示例代码：pgmpy官方链接

这个库使得处理概率图模型变得更加容易，可以帮助你分析和推断各种依赖关系。下面是一个构建贝叶斯网络并且进行学习的示例代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
# encoding: utf-8
# ====================BN模型=========================
# 贝叶斯模型
from pgmpy.models import BayesianNetwork
# ====================参数学习=========================
# 参数估计
from pgmpy.estimators import ParameterEstimator
# MLE参数估计
from pgmpy.estimators import MaximumLikelihoodEstimator
# Bayesian参数估计
from pgmpy.estimators import BayesianEstimator
# ====================结构学习=========================
# ========评分搜索=========================
# 评分
from pgmpy.estimators import BDeuScore, K2Score, BicScore
# 穷举搜索
from pgmpy.estimators import ExhaustiveSearch
# 爬山搜索
from pgmpy.estimators import HillClimbSearch
# ========  混合  =========================
from pgmpy.estimators import MmhcEstimator
# ==================== 通用库 =========================
import pandas as pd
import numpy as np
import random
from sklearn.metrics import accuracy_score


# 设置随机种子以便复现
random.seed(0)

# 生成示例数据
data = {
    'average_weekly_hours': [random.choice([-1, 0, 1]) for _ in range(1000)],
    'unemployment': [random.choice([-1, 1]) for _ in range(1000)],
    'eci_leading_indicator': [random.choice([-1, 1]) for _ in range(1000)],
    'manufacturing_pmi_inflation': [random.choice([-1, 1]) for _ in range(1000)],
    '10yr_treasury_yield': [random.choice([-2, -1, 0, 1, 2]) for _ in range(1000)],
    'inflation_expectation': [random.choice([-2, -1, 0, 1, 2]) for _ in range(1000)],
    'comex_gold_settlement': [random.choice([-2, -1, 0, 1, 2]) for _ in range(1000)]
}

# 创建DataFrame
df = pd.DataFrame(data)
# 打印前几行数据
print(df.head())

# 变量附录
# 平均每周工时变化量：average_weekly_hours
# 新增失业人数变化量：unemployment
# ECI领先指标变化量：eci_leading_indicator
# 制造业PMI物价变化量：manufacturing_pmi_inflation
# 10年期国债收益率变化量：10yr_treasury_yield
# 通胀预期变化量：inflation_expectation
# COMEX黄金结算价变化量：comex_gold_settlement
node_tree = [
    ('average_weekly_hours', '10yr_treasury_yield'), ('average_weekly_hours', 'inflation_expectation'),
    ('unemployment', '10yr_treasury_yield'), ('unemployment', 'inflation_expectation'),
    ('eci_leading_indicator', '10yr_treasury_yield'), ('eci_leading_indicator', 'inflation_expectation'),
    ('manufacturing_pmi_inflation', '10yr_treasury_yield'), ('manufacturing_pmi_inflation', 'inflation_expectation'),
    ('10yr_treasury_yield', 'comex_gold_settlement'), ('inflation_expectation', 'comex_gold_settlement'),
]
model = BayesianNetwork(node_tree)

train_data = df[:800]
predict_data = df[800:]
# 用于统计正确率
actual_values = predict_data['comex_gold_settlement'].values
print(actual_values)

mle = MaximumLikelihoodEstimator(model, df)
# print(mle.estimate_cpd('fruit'))  # unconditional
# print(mle.estimate_cpd('tasty'))  # conditional
model.fit(df, estimator=MaximumLikelihoodEstimator)

predict_data = predict_data.copy()
predict_data.drop('comex_gold_settlement', axis=1, inplace=True)
y_pred = model.predict(predict_data)
print(y_pred)

# 计算准确率
accuracy = accuracy_score(actual_values, y_pred)

print(f"Accuracy: {accuracy * 100:.2f}%")

学习参数

学习参数通常涉及从一组数据样本和一个有向无环图（DAG）中估计单个变量的概率分布。虽然具体的代码被省略，但通常可以使用统计方法或概率模型来进行参数学习。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
def parameterLearning():
    # ================= 数据部分 ======================
    data = pd.DataFrame(data={'fruit': ["banana", "apple", "banana", "apple", "banana", "apple", "banana",
                                        "apple", "apple", "apple", "banana", "banana", "apple", "banana", ],
                              'tasty': ["yes", "no", "yes", "yes", "yes", "yes", "yes",
                                        "yes", "yes", "yes", "yes", "no", "no", "no"],
                              'size': ["large", "large", "large", "small", "large", "large", "large",
                                       "small", "large", "large", "large", "large", "small", "small"]})

    # ================= 构建网络模型结构 ======================
    model = BayesianNetwork([('fruit', 'tasty'), ('size', 'tasty')])  # fruit -> tasty <- size

    # ================= 选择合适的参数估计器 =====================
    # 参数估计（Parameter Estimation）：
    #
    # 参数估计是一种直接从数据集中计算条件概率分布的参数的方法。
    # 它不使用任何先验信息，仅根据观察的数据来估计参数。
    # 参数估计器通常用于有限数据集和在没有领域专业知识的情况下。
    # 最大似然估计（Maximum Likelihood Estimation, MLE）：
    #
    # MLE是一种参数估计方法，旨在最大化观察到的数据发生的可能性。
    # 它基于样本数据计算参数，使观察到的数据在给定模型下的概率最大。
    # MLE通常用于大型数据集和参数空间较大的情况，因为它不需要先验知识。
    # 贝叶斯参数估计（Bayesian Parameter Estimation）：
    #
    # Bayesian参数估计是一种使用贝叶斯方法估计参数的方法，它结合了先验知识和观察数据。
    # 它基于贝叶斯定理，将参数的后验分布计算为先验分布和似然性的乘积。
    # Bayesian参数估计通常用于小型数据集和需要领域专业知识的情况，以指导先验的选择。
    # 比较：
    #
    # 参数估计是一种通用方法，但在某些情况下可能不够准确，因为它不考虑先验信息。
    # MLE是无信息估计，通常用于大型数据集，但对数据不足或参数空间复杂的情况不够稳健。
    # Bayesian参数估计结合了先验知识和观察数据，通常更准确和稳健，但需要领域专业知识来指导先验的选择。
    print("========================================================")
    pe = ParameterEstimator(model, data)
    print("\n", pe.state_counts('fruit'))  # unconditional
    print("\n", pe.state_counts('size'))  # unconditional
    print("\n", pe.state_counts('tasty'))  # conditional on fruit and size
    print("========================================================")
    mle = MaximumLikelihoodEstimator(model, data)
    print(mle.estimate_cpd('fruit'))  # unconditional
    print(mle.estimate_cpd('tasty'))  # conditional

    print("========================================================")
    est = BayesianEstimator(model, data)
    print(est.estimate_cpd('tasty', prior_type='BDeu', equivalent_sample_size=10))
    # Setting equivalent_sample_size to 10 means
    # that for each parent configuration, we add the equivalent of 10 uniform samples
    # (here: +5 small bananas that are tasty and +5 that aren't).

    print("========================================================")
    # Calibrate all CPDs of `model` using MLE:
    model.fit(data, estimator=MaximumLikelihoodEstimator)
    print("========================================================")


def structuralLearning_Score():
    """
        score-based structure learning
        constraint-based structure learning
    The combination of both techniques allows further improvement:
        hybrid structure learning
    """
    print("===================基于评分=================================")
    # create random data sample with 3 variables, where Z is dependent on X, Y:
    data = pd.DataFrame(np.random.randint(0, 4, size=(5000, 2)), columns=list('XY'))
    data['Z'] = data['X'] + data['Y']

    bdeu = BDeuScore(data, equivalent_sample_size=5)
    k2 = K2Score(data)
    bic = BicScore(data)

    model1 = BayesianNetwork([('X', 'Z'), ('Y', 'Z')])  # X -> Z <- Y
    model2 = BayesianNetwork([('X', 'Z'), ('X', 'Y')])  # Y <- X -> Z
    print("==========基于评分===model1===============")
    print(bdeu.score(model1))
    print(k2.score(model1))
    print(bic.score(model1))
    print("==========基于评分===model2===============")
    print(bdeu.score(model2))
    print(k2.score(model2))
    print(bic.score(model2))
    print("==========基于评分===局部评分==============")
    print(bdeu.local_score('Z', parents=[]))
    print(bdeu.local_score('Z', parents=['X']))
    print(bdeu.local_score('Z', parents=['X', 'Y']))
    print("==========基于评分===穷举搜索算法==============")
    # 穷举搜索（计算困难），启发式搜索
    es = ExhaustiveSearch(data, scoring_method=bic)
    # 获取分数最高的分数
    best_model = es.estimate()
    print(best_model.edges())
    print("\n 遍历所有的分数:")
    for score, dag in reversed(es.all_scores()):
        print(score, dag.edges())
    print("==========基于评分===爬山搜索算法==============")
    data = pd.DataFrame(np.random.randint(0, 3, size=(2500, 8)), columns=list('ABCDEFGH'))
    data['A'] += data['B'] + data['C']
    data['H'] = data['G'] - data['A']
    hc = HillClimbSearch(data, scoring_method=BicScore(data))
    best_model = hc.estimate()
    print(best_model.edges())


def structuralLearning_Hybrid():
    """
    MMHC算法[3]结合了基于约束和基于分数的方法。它有两部分:
        1. 使用基于约束的构造过程MMPC学习无向图骨架
        2. 基于分数的优化(BDeu分数+修改爬山)
    """
    print("===================混合方法=================================")
    # 实验数据生成
    data = pd.DataFrame(np.random.randint(0, 3, size=(2500, 8)), columns=list('ABCDEFGH'))
    data['A'] += data['B'] + data['C']
    data['H'] = data['G'] - data['A']
    data['E'] *= data['F']
    # 构建无向图骨架
    mmhc = MmhcEstimator(data)
    skeleton = mmhc.mmpc()
    print("Part 1) Skeleton: ", skeleton.edges())
    # 基于分数优化
    # use hill climb search to orient the edges:
    hc = HillClimbSearch(data, scoring_method=BDeuScore(data))
    model = hc.estimate(tabu_length=10, white_list=skeleton.to_directed().edges())
    print("Part 2) Model:    ", model.edges())
    print("===================两步划为一步=================================")
    # MmhcEstimator.estimate(self, scoring_method=None, tabu_length=10,
    #       significance_level=0.01)
    # mmhc.estimate(scoring_method=BdeuScore(data),tabu_length=10)

绘制概率图

为了可视化概率图模型，我们需要使用工具来绘制图像。在这里，我们使用 daft 中的 PGM，同时使用 textwrap 进行文本的自动分行。这有助于我们清晰地呈现模型的结构和依赖关系。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
%matplotlib inline
import daft
from daft import PGM
import matplotlib.pyplot as plt
import textwrap


# 设置字体为 SimHei（中文黑体）
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号

# 7个节点 9条边
pgm = PGM(shape=[9, 5.4])

def fit_node(text, x, y, min_width=2, min_height=1):
    # 使用 textwrap 模块来自动分行
    wrapped_text = textwrap.fill(text, width=4)  # 设置每行最大宽度
    
    # 计算节点宽度和高度，根据文本行数自动调整
    num_lines = len(wrapped_text.split('\n'))
    node_width = max(len(wrapped_text) // num_lines, min_width)
    node_height = max(0.5 * num_lines, min_height)
    
    pgm.add_node(daft.Node(text, wrapped_text, x, y, aspect=node_width/node_height))
    
fit_node("平均每周工时", 1, 5)
fit_node("新增失业人数", 3, 5)
fit_node("ECRI领先指标", 5, 5)
fit_node("制造业PMI物价", 7.4, 5)
fit_node("名义利率", 2.5, 3)
fit_node("通货预期", 5.5, 3)
fit_node("黄金价格", 4, 1)


pgm.add_edge("平均每周工时", "名义利率")
pgm.add_edge("平均每周工时", "通货预期")
pgm.add_edge("新增失业人数", "名义利率")
pgm.add_edge("新增失业人数", "通货预期")
pgm.add_edge("ECRI领先指标", "名义利率")
pgm.add_edge("制造业PMI物价", "名义利率")
pgm.add_edge("制造业PMI物价", "通货预期")
pgm.add_edge("名义利率", "黄金价格")
pgm.add_edge("通货预期", "黄金价格")

pgm.render()
plt.show()

未来，我们还可以考虑使用 bnlearn 库，这是一个专门用于学习贝叶斯网络的图形结构、参数学习、推理和抽样方法的 Python 包。值得注意的是，bnlearn 建立在 pgmpy 包的基础上，提供更多功能和工具。

如果你对学习更多关于贝叶斯网络和如何使用这些库感兴趣，可以查看以下链接：

参考链接：

Python的贝叶斯网络学习库pgmpy介绍和使用-CSDN博客

贝叶斯网络——pgmpy 教程 - bonelee - 博客园 (cnblogs.com)

pgmpy——高斯贝叶斯网络-CSDN博客

使用pgmpy构建贝叶斯网络并且进行参数学习
https://fulequn.github.io/2023/10/Article202310301/
作者
Fulequn
发布于
2023年10月30日
许可协议