使用Python实现马尔科夫随机场、蒙特卡洛采样等随机过程算法的前提，就是用Python实现概率的计算。并不只是数值计算，而是能够将随机模拟中常用的各种概率相关的操作，都能用计算机的数据结构来表达，其关键在于对【随机变量】的适当定义处理。
因此本文介绍一下概率分布在Python中定义的一种数据结构。
一个概率分布的组成要素包含：随机变量、变量的维度、变量不同取值状态的对应概率值。
在一个有向图中(贝叶斯网络)，常有如下所示的结构：
典型的贝叶斯网络
在这个贝叶斯网络中，有五个随机变量，各自有不同的维度，并且在列表中给出了对应变量取值的概率值。因此我们的目标可以拆解为两步，第一：
用Python精准地定义一个因子对象(记为Factor)，这个对象表达一个联合分布，其中包含随机变量、变量维度、对应概率值这三个信息。
第二：
实现概率相乘、求边缘概率、求条件概率运算。
这次先解决问题一。
将概率分布间的计算用python表达出来，包括对一个联合分布求某个变量的边缘概率分布，对一个联合分布给定一部分变量的取值，求条件概率分布。以及更基础的，求两个因子的乘积-状态维度和状态维度的相乘，所对应的概率值的数量也同样扩大。
最后，分别获取这个贝叶斯网络中每个因子的概率分布，我们要求出整个网络的联合分布：
求联合分布
联合分布概率为每个因子在父节点条件下条件概率的乘积。
具体实现如下：
首先实现因子对象，Factor，其中随机变量用属性var(variables)表示，为python_list结构，变量维度card(cardinality)，同样为list，概率取值val(value)列表list。对因子对象实现获取状态结构函数(get_all_assignments),正态化函数(将val列表中所有概率值正态化，使其和为1)。
注意，val概率值列表的顺序，对应着变量取值的某个特定状态。
实现打印功能，方便查看一个因子的具体情况。
class Factor: def __init__(self, var=None, card=None, val=None): if var is None: var = np.array([], np.int64) if card is None: card = np.array([], np.int64) if val is None: val = np.array([], np.float64) self.var = np.array(var) self.card = np.array(card) self.val = np.array(val) def is_empty(self): """Returns true if the factor is empty (i.e. not initialized)""" return len(self.var) == 0 def get_all_assignments(self): assignments = index_to_assignment(np.arange(int(np.prod(self.card))), self.card) return assignments def __repr__(self): if self.is_empty(): str = 'Empty factor' else: str = 'Factor containing {} variables'.format(len(self.var)) num_states = int(np.prod(self.card)) assigments = index_to_assignment(np.arange(num_states), self.card) # Print header row header = '| ' + ' '.join(['X_{}'.format(i) for i in self.var]) + ' | Probability |' col_width = len(header) line = '-' * col_width + '' str += line + header + '' + line # Assignments for i in range(assigments.shape[0]): lhs = ' '.join(['{}'.format(a) for a in assigments[i]]) row = '| ' + lhs + ' | ' + '{:>11g}'.format(self.val[i]) + ' |' str = str + row str += line + '' return str def normalize(self): """Normalize the probablity such that it sums to 1. Use this function with care since not all factor tables should be normalized. """ if not self.is_empty(): self.val /= np.sum(self.val)比如，定义一个因子:
card=[2,3,2]phi = Factor(var=[1, 2,5], card=card, val=[0.6, 0.3, 0.2, 0.7, 0.2, 0.0,0.1,0.2,0.1,0.1,0.2,0.3])使得该因子含3个随机变量，随机变量维度分别为2，3，2，概率值如val列表所示。
将其打印，结果如下：
print(phi)———————————————Factor containing 2 variables-------------------------| X_1 X_2 X_5 | Probability |-----------------------------| 0 0 0 | 0.6 || 1 0 0 | 0.3 || 0 1 0 | 0.2 || 1 1 0 | 0.7 || 0 2 0 | 0.2 || 1 2 0 | 0 || 0 0 1 | 0.1 || 1 0 1 | 0.2 || 0 1 1 | 0.1 || 1 1 1 | 0.1 || 0 2 1 | 0.2 || 1 2 1 | 0.3 |-----------------------------这样就完成了问题一，实现了概率分布这个数据结构。这样做的意义，在于将随机变量的信息和概率值绑定在一起，为后续的概率计算，乃至各种随机模拟算法的实现提供条件。
为了后续的实现，开发两个因子函数，用以方便地查看概率值，及其所对应的变量状态取值：
def index_to_assignment(index: Union[int, np.ndarray], card: np.ndarray): """Convert index to variable assignment Args: index: Index to convert into assignment. If index is a vector of numbers, the function will return a matrix of assignments, one assignment per row. card: Cardinality of the factor """ if isinstance(index, int): is_scalar = True index = np.array([index]) else: is_scalar = False divisor = np.cumprod(np.concatenate([[1.], card[:-1]])) assignment = np.mod( np.floor(index[:, None] / divisor[None, :]), card[None, :] ).astype(np.int64) if is_scalar: assignment = assignment[0] return assignmentdef assignment_to_index(assignment: np.ndarray, card: np.ndarray): """Convert assignment to index. Args: assignment: Assignment to convert to index card: Cardinality of the factor """ multiplier = np.cumprod(np.concatenate([[1.], card[:-1]])) index = np.sum(assignment * multiplier, axis=-1).astype(np.int64) return index以上是用Python定义概率分布数据结构的方法。
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
