互联学术本问题例取自于Ref1-Chapter1.
问题描述：假设有一个制作灯泡的机器。你想知道机器是正常工作还是有问题。为了得到答案你可以测试每一个灯泡，但是灯泡数量很多，每一个都测试在实际生产过程中可能是无法承受的。使用贝叶斯推断，你可以基于少量样本（比如说抽检结果）来估计机器是否在正常地工作（in probabilistic way）。
构建贝叶斯推断时，首先需要两个要素:
(1) 先验分布
(2) 似然率
先验分布是我们关于机器工作状态的初始信念。首先我们确定第一个刻画机器工作状态的随机变量，记为M。这个随机变量有两个工作状态：{working, broken}，以下简写成{w, br}(缩写成br是为了与下面的Bad缩写成b区分开来)。作为初始信念，我们相信机器是好的，是可以正常工作的，定义先验分布如下：
P(M=working) = 0.99
P(M=broken ) = 0.01
这表明我们对于机器正常工作的信念度很高，有99%的概率能够正常工作。
第二个随机变量是L，表示机器生产的灯泡的工作状态。灯泡可能是好，也可能是坏的，包含两个状态：{good, bad}，以下简写成{g,b}，注意br与b的区别。
我们需要基于机器工作状态给出L的先验分布，也就是条件概率P(L|M)，在贝叶斯公式中它代表似然概率(likelihood)。
定义这个似然概率分布（由于M和L各有两种状态，所以一共包含4个条件概率）如下：
P(L=Good|M=w) = 0.99
P(L=Bad |M=w) = 0.01
P(L=Good|M=br ) = 0.6
P(L=Bad |M=br ) = 0.4
现在，我们已经完整地刻画了贝叶斯模型，可以用它来做一些神奇的估计和预测的工作了。
我们的输入是一些灯泡的抽检结果。假设我们抽检了十个灯泡其抽检结果如下：
{bad, good, good, good, good, good, good, good, good, good}
让我们来看看基于贝叶斯推断的我们对于机器工作状态的信念（后验概率）如何变化。
3. 贝叶斯规则
贝叶斯推断规则以贝叶斯公式的形式表示为：
具体映射到本问题中可以表达如下：
贝
叶斯推断的优先在于可以以在线（online）的方式进行，即观测数据可以一个一个地到来，每次受到一个新的观测数据，就进行一次基于贝叶斯公式的后验概率的计算更新，而更新后的后验概率又作为下一贝叶斯推断的先验概率使用。因此在线的贝叶斯推断的基本处理流程如下所示：
4. Bayes engine: scalar implementation
首先，我们以标量运算的方式写一个函数来进行bayes推断处理。
prior以向量的形式存储先验概率分布，prior[0]表示P(M=working)，prior[1]表示P(M=broken)。
likelihood以矩阵的形式方式存储似然概率分布。其中第1行表示P(L/M=working)，第2行表示P(L/M=broken).
在本例中，当输入
时，evidence的计算式（注意evidence是依赖于输入的观测数据的）是:
注意，当我写P(w)其实是表示P(M=w)，而
其实是表示
，余者类推。根据上下文，这些应该不会导致混淆。
第一个函数的代码如下：
def bayes_scalar(prior, likelihood, data): """ Bayesian inference function example. Parameters ---------- prior : float, 1-D vector prior information, P(X). likelihood : float 2-D matrix likelihood function, P(Y|X). data : List of strings. Value: 'Good','Bad' Observed data samples sequence Returns ------- posterior : float P(X,Y), posterior sequence. """posterior = np.zeros((len(data)+1,2)) posterior[0,:] = prior# Not used in computation, just for the later plotting for k,L in enumerate(data): if L == 'good': L_value = 0 else: L_value = 1 #print(L, L_value, likelihood[:,L_value]) evidence= likelihood[0,L_value] * prior[0] + likelihood[1,L_value] * prior[1] LL0_prior_prod= likelihood[0,L_value] * prior[0] posterior[k+1,0]= LL0_prior_prod / evidence LL1_prior_prod= likelihood[1,L_value] * prior[1] posterior[k+1,1]= LL1_prior_prod / evidenceprior = posterior[k+1,:] # Using the calculated posterior at this step as the prior for the next stepreturn posterior
5.Bayes engine: vectorization
我们注意到，evidence的计算可以表示成两个向量的点积，如下所示。
这样就非常方便用numpy来实现了。本例中每个随机变量只有两种取值，在复杂的情况下，每个随机变量有很多种取值时，有效利用向量或矩阵的运算是简洁的运算实现的必不可缺的要素。以上这两个向量的点积可以用numpy.dot()来实现。
另外，likelihood和prior的乘积是分别针对M的两种状态进行计算（注意，我们需要针对M的两种不同状态分别计算posterior），不是用向量的点积进行计算，而是一种element-wise multiplication，可以用numpy.multiply()进行计算。所以在vectorization版本中贝叶斯更新处理削减为两条语句，与上面的scalar版本相比显得非常优雅简洁（好吧，也许这个简单例子中还显不出那么明显的优势，但是随着问题的复杂度的增加，这种优势就会越来越明显了。）
由此我们得到向量化处理的函数如下：
def bayes_vector(prior, likelihood, data): """ Bayesian inference function example. Parameters ---------- prior : float, 1-D vector prior information, P(X). likelihood : float 2-D matrix likelihood function, P(Y|X). data : List of strings. Value: 'Good','Bad' Observed data samples sequence Returns ------- posterior : float P(X,Y), posterior sequence. """posterior = np.zeros((len(data)+1,2)) posterior[0,:] = prior# Not used in computation, just for the later plotting for k,L in enumerate(data): if L == 'good': L_value = 0 else: L_value = 1 #print(L, L_value, likelihood[:,L_value]) evidence= np.dot(likelihood[:,L_value], prior[:]) posterior[k+1,:]= np.multiply(likelihood[:,L_value],prior)/evidence prior = posterior[k+1,:] # Using the calculated posterior at this step as the prior for the next stepreturn posterior
6. 测试
让我们来看看利用以上函数对我们的观测数据进行处理，后验概率将会如何变化。
import numpy as np
import matplotlib.pyplot as plt
prior= np.array([0.99,0.01])
likelihood = np.array([[0.99,0.01],[0.6,0.4]])
data = ['bad','good','good','good','good','good','good','good']
posterior1 = bayes_scalar(prior,likelihood,data)
posterior2 = bayes_vector(prior,likelihood,data)
if np.allclose(posterior1,posterior2):
print('posterior1 and posterior2 are identical!')
fig, ax = plt.subplots()
ax.plot(posterior1[:,0])
ax.plot(posterior1[:,1])
ax.grid()
# fig.suptitle('Poeterior curve vs observed data')
ax.set_title('Posterior curve vs observed data')
plt.show()
运行以上代码可以得到后验概率的变化如下图所示（注意第一个点是prior）：
