用Python模拟贝叶斯信念网络
文章目录
简介问题描述用Python构建贝叶斯信念网络模型识别变量识别变量的可能取值创建网络的心智图创建Python函数完成函数细节组合在一起用BBN网络进行推断简介
本文将向你展示如何利用Python构建简单的贝叶斯信念网络，并用它来进行严格的推理。
我们要建模的问题是著名的 蒙提霍尔问题 (也叫三门问题)。选择蒙提霍尔问题是因为它具有如下3个特征，这3个特征使得它非常适合作为教程举例：
它很简单，只有3个变量；它的结论非常有趣且反经验；它的条件概率很简单，只需要几行python代码即可实现。
问题描述
蒙提霍尔问题来自一个电视游戏节目：主持人蒙提霍尔邀请一位嘉宾选择3扇门中的其中一扇，并告知嘉宾三扇门中有一扇后面放的是汽车，其余两扇后面放的是山羊。主持人知道每扇门后面放的是什么但是他不会告诉嘉宾。
备注
我们假设嘉宾对门后的物品没有任何先验知识。同时我们假设嘉宾都想获得汽车而不是山羊。
当嘉宾做出了初次选择后，为了让节目更加有趣，主持人会先打开另外2扇门中的一扇。
备注
因为主持人知道哪扇门后有汽车，他绝不会打开有汽车的门，而是打开那个有山羊的门。
当打开了其中一扇没有汽车的门后，主持人会给嘉宾以此选择的机会–是否换一扇门。
我们要解答的问题是：嘉宾应该坚持最初的选择还是应该换一扇门？
很多人认为换不换都一样。然而事实真的如此吗？接下来我们用贝叶斯信念网络来推理一下这个问题。
用Python构建贝叶斯信念网络模型
构建贝叶斯信念网络非常简单直接，只要按照规则和约定一步一步走即可。
识别变量
首先，对于任何要建模的问题，识别变量（随机变量）是第一步。对蒙提霍尔问题来说，识别过程中的事件有助于我们找到随机变量。
蒙提霍尔问题中的事件按照时间发生顺序如下：
汽车事先藏在3扇门中的一扇；嘉宾做出初次选择；主持人打开另外两扇门中的一扇；
上面每一个事件都是模型中的一个变量，我们将这些变量命名为：
prize_doorguest_doormonty_door
识别变量的可能取值
当创建贝叶斯网络时，我们不仅需要识别变量还需要知道变量可能的取值。
本案例比较简单：有3扇门，我们将其命名为 A , B , C A, B, C A,B,C。因此三个变量可以取 [ A , B , C ] [A, B, C] [A,B,C] 中的任意值。
注意
我们称变量可取值的集合为变量的定义域。不是所有变量总有相同的定义域。
创建网络的心智图
现在我们可以开始构建贝叶斯信念网络(BBN)。贝叶斯信念网络是个有向无环图（DAG），这意味着我们需要构建一个包含节点（也叫顶点）和有向边的图。
每一个变量将作为网络中的节点。我们可以在图中放置3个节点表示每一个变量：
接下来我们需要添加边。边在概率图模型中表示边连接的两个节点其中一个影响另外一个的事实。在贝叶斯信念网络中，由于边是有向的，这意味节点关系是"父子"关系，即子节点条件依赖父节点。这意味着子节点的取值依赖父节点的取值。考虑蒙提霍尔问题中的事件，当支持人打开后面是样的门时，他的决策依赖车藏在哪个门后（因为他永远不能暴露奖品），当然主持人也不可能选择嘉宾选择的门，这有悖游戏规则。于是我们需要从prize_door节点连一条边到monty_door，同理我们需要从guest_door节点连一条边到monty_door。完成后的贝叶斯信念网络的结构如下图：
创建Python函数
现在我们有了网络结构图，接下来就可以开始编程了。图中所有节点用一个Python函数来表示。函数的参数表示模型中的变量。以prize_door节点为例，它的python代码为:
def node_prize_door(prize_door):pass
同理，guest_door节点代码为：
def node_guest_door(prize_door):pass
对monty_door节点来说，python代码稍有不同，因为这个节点有2条入边，我们需要将prize_door和guest_door两个变量加入到函数的参数列表中：
def node_monty_door(prize_door, guest_door, monty_door):pass
对于上面创建的三个函数，有几点需要说明：
函数的参数名与模型中的变量名一致；函数代表图中的节点；为了区分节点名和变量名，我在函数名加了node_前缀；如果节点有父节点，则需要将父节点的变量添加到参数列表中；三个函数对应三个变量，三个变量代表模型中的按个事件；每个函数引入一个新变量到模型中。
完成函数细节
贝叶斯信念网络不仅需要图，还需要分布来控制变量。接下来我们将加入每个变量的先验和条件概率。要完成这一步可以用下面的方法思考函数：
每一个函数必须返回[0, 1]之间的一个实数，代表这个变量取其定义域上某个值的概率。
以node_prize_door为例，我们想返回奖品藏在A，B，C某扇门后的概率。假设奖品是随机藏在门后的，那么奖品出现在每扇门后的概率是1/3。我们需要将node_prize_door函数的返回值修改为1/3:
def node_prize_door(prize_door):return 0.33333333
同理，由于嘉宾没有先验的知识，他的选择也是随机的，所以嘉宾选择某扇门的概率也是1/3:
def node_guest_door(prize_door):return 0.33333333
node_monty_door要有趣一些，我们需要计算在给定父节点prize_door和guest_door的情况下monty_door变量取A, B, C三个值中每一个的似然值。如果嘉宾有幸选对了，支持人可以随意选择剩下的2扇门；如果嘉宾选错了，此时支持人既不能选择嘉宾选择的门也不能选择有奖品的门，那么他只有一个选择：两扇门中有山羊的门。python代码如下：
def node_monty_door(prize_door, guest_door, monty_door):if guest_door == prize_door:# 嘉宾选对了if monty_door == prize_door:return 0# 主持人绝不会暴露奖品else:return 0.5# 随意选择其余2扇门中的一扇elif monty_door == prize_door:return 0# 主持人绝不会暴露奖品elif monty_door == guest_door:return 0# 支持人也绝不会选择嘉宾选择的门else:return 1# 嘉宾没有猜对，支持人只能选择剩下两扇门中有山羊的那扇
组合在一起
接下来我们要完成整个代码并在其之上运行一些推断。要创建图，我们需要调用build_bnn模块。在程序顶部加入下面一行代码经模块引入程序中：
from bayesian.bbn import build_bbn
接下来加入主函数
if __name__ == "__main__":g = build_bnn(node_prize_door, node_guest_door, node_monty_door,domains=dict(prize_door=['A', 'B', 'C'],guest_door=['A', 'B', 'C'],monty_door=['A', 'B', 'C']))
上面的代码创建了一个贝叶斯信念网络的实例。工厂方法build_bbn接受所有代表图中节点的函数作为参数，附加一个可选参数domains，domains接受一个字典用于指明每个变量的取值范围。注意：图的结构是用过函数和参数推理得到的。
将整个python代码保存为monty_hall.py文件
用BBN网络进行推断
执行下面命令
$ python -i monty_hall.py
-i参数指示python解释器在执行完给定python脚本后运行在交互模式下。这样我们接下来访问模型会更方便。BBN类主要通过query方法来查询。库中提供了一个用户体验友好的查询方法–名叫q。我们可以不带任何参数尝试调用一下q函数，看会输出什么：
>>> g.q()+------------+-------+----------+| Node | Value | Marginal |+------------+-------+----------+| guest_door | A| 0.333333 || guest_door | B| 0.333333 || guest_door | C| 0.333333 || monty_door | A| 0.333333 || monty_door | B| 0.333333 || monty_door | C| 0.333333 || prize_door | A| 0.333333 || prize_door | B| 0.333333 || prize_door | C| 0.333333 |+------------+-------+----------+>>>
要如何解读上面的输出呢？q方法基本上就是用相同的参数调用query方法，然后将query方法返回的结果格式化成易读的表格。表格的列为节点(Node)、值(Value)和边缘概率(Marginal)。你会注意到表格显示的是每一个节点的每一个可取值的边缘概率(Marginal)。表格中所有的边缘概率都是0.333333。这是因为我们没有为模型提供任何额外的证据。在缺乏证据的情况下，模型估计每个门的概率是1/3。
注意
在BBN的术语中我们称给变量赋值为证据。我们称任意为0个或多个变量赋值的集合为图的配置。
接下来让我们为BBN提供一些证据，然后再查询一遍。假设我们观察到嘉宾选择了A门，输入下面的python交互指令：
>>> g.q(guest_door='A')+-------------+-------+----------+| Node | Value | Marginal |+-------------+-------+----------+| guest_door | B| 0.000000 || guest_door | C| 0.000000 || guest_door* | A* | 1.000000 || monty_door | A| 0.000000 || monty_door | B| 0.500000 || monty_door | C| 0.500000 || prize_door | A| 0.333333 || prize_door | B| 0.333333 || prize_door | C| 0.333333 |+-------------+-------+----------+>>>
注意当我们观察到嘉宾选择了A门后边缘概率的变化。变量guest_door取值A的概率变为1，这是理所当然的，因为我们本来就观察到嘉宾选择了A门。同理变量guest_door取B和C的概率变为0，因为我们知道嘉宾一定不会选B门和C门。同时我们注意到作为证据的变量在表格里标记上了星号，以提示我们为本次查询提供了哪些证据。接下来看monty_door变量。monty_door变量取值A的概率变为0，这是因为游戏规则设定为主持人不能选择嘉宾选择的门。monty_door变量取B和取C的概率一半一半，表示在没有其他额外证据的情况下主持人选择打开剩下的两扇门的中哪一扇的概率是相等的。
如果此时我们观察到支持人打开了门B，并询问嘉宾要不要换门。将新的证据加入后再查询得到：
>>> g.q(guest_door='A', monty_door='B')+-------------+-------+----------+| Node | Value | Marginal |+-------------+-------+----------+| guest_door | B| 0.000000 || guest_door | C| 0.000000 || guest_door* | A* | 1.000000 || monty_door | A| 0.000000 || monty_door | C| 0.000000 || monty_door* | B* | 1.000000 || prize_door | A| 0.333333 || prize_door | B| 0.000000 || prize_door | C| 0.666667 |+-------------+-------+----------+>>>
观察prize_door的边缘概率，里面包含了我们最初提出问题的答案。显而易见，奖品在C门后的概率是A门后的2倍。
