使用pgmpy构建贝叶斯网络并且进行参数学习
pgmpy - Python中的概率图模型库
PGMPY是Python中的一个概率图模型库，用于创建、操作和推断概率图模型。它提供了一组工具，用于参数学习和结构学习，以帮助我们理解变量之间的依赖关系。以下是一些主要概念：
-
参数学习：当我们有一组数据样本和一个有向无环图（DAG）来表示变量之间的依赖关系时，参数学习旨在估计单个变量的（条件）概率分布。
-
结构学习：当我们有一组数据样本时，结构学习旨在估计一个DAG，该DAG 捕获变量之间的依赖关系。
注意：在pgmpy的最新版本0.1.23中，
BayesianModel已经更名为BayesianNetwork。
你可以在官方链接中找到更多信息和示例代码：pgmpy官方链接
这个库使得处理概率图模型变得更加容易，可以帮助你分析和推断各种依赖关系。下面是一个构建贝叶斯网络并且进行学习的示例代码
# 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）中估计单个变量的概率分布。虽然具体的代码被省略，但通常可以使用统计方法或概率模型来进行参数学习。
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 进行文本的自动分行。这有助于我们清晰地呈现模型的结构和依赖关系。
%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构建贝叶斯网络并且进行参数学习
https://blog.fulequn.online/2023/10/Article202310301/