《贝叶斯方法：概率编程与贝叶斯推断》一1.6　补充说明

本节书摘来自异步社区《贝叶斯方法：概率编程与贝叶斯推断》一书中的第1章，第1.6节，作者 【加】Cameron Davidson-Pilon（卡梅隆 戴维森-皮隆），更多章节内容可以访问云栖社区“异步社区”公众号查看

1.6　补充说明

1.6.1　从统计学上确定两个λ值是否真的不一样

在短信接收例子中，我们直观地观测了λ1和λ2的先验信息并认为它们是不同的。这很公平，毕竟先验的位置基本离得非常远。但如果这并不是真相，有一部分分布其实是重合的呢？我们怎么才能让上面的结论更加的正式呢？

一种方法就是计算出P(λ1<λ2|data)，即在获得观察数据的前提下，λ1的真实值比λ2小的概率。如果这个概率接近50%，那相当于抛硬币得到的结果，这样我们就不能确定它们是否真的不同。如果这个概率接近100%，那么我们就能确定地说这两个值是不同的。利用后验中的样本，这种计算非常简单——我们计算λ1后验中的样本比λ2后验中的样本小的次数占比：

print (lambda_1_samples < lambda_2_samples) # Boolean array: True if lambda_1 is less than lambda_2. [Output]: 　 [ True True True True ..., True True True True] # How often does this happen? print (lambda_1_samples < lambda_2_samples).sum() 　 # How many samples are there? print lambda_1_samples.shape[0] [Output]: 　 29994 30000 # The ratio is the probability. Or, we can just use .mean: print (lambda_1_samples < lambda_2_samples).mean() [Output]: 　 0.9998

结果很显然，有几乎100%的把握可以说这两个值是不等的。

我们也可以再问详细一点，比如：“两个值之间相差1、2、5、10的概率有多大？”

# The vector abs(lambda_1_samples - lambda_2_samples) > 1 is a boolean, # True if the values are more than 1 apart, False otherwise. # How often does this happen? Use .mean() for d in [1,2,5,10]: v = (abs(lambda_1_samples - lambda_2_samples) >= d).mean() print "What is the probability the difference is larger than %d\ ? %.2f"%(d,v) [Output]: 　 What is the probability the difference is larger than 1? 1.00 What is the probability the difference is larger than 2? 1.00 What is the probability the difference is larger than 5? 0.49 What is the probability the difference is larger than 10? 0.00

1.6.2　扩充至两个转折点

读者们或许会对前面模型中转折点个数的扩充，即如果不止一个转折点会怎么样感兴趣，或者会对只有一个转折点的结论表示怀疑。下面我们把模型扩充至两个转折点（意味着会出现3个λi）。新模型跟之前的比较相像。

我们把这个模型也编译成代码，跟前面的代码看上去差不多。

lambda_1 = pm.Exponential("lambda_1", alpha) lambda_2 = pm.Exponential("lambda_2", alpha) lambda_3 = pm.Exponential("lambda_3", alpha) 　 tau_1 = pm.DiscreteUniform("tau_1", lower=0, upper=n_count_data-1) tau_2 = pm.DiscreteUniform("tau_2", lower=tau_1, upper=n_count_data) 　 @pm.deterministic def lambda_(tau_1=tau_1, tau_2=tau_2, lamb-da_1=lambda_1, lambda_2=lambda_2, lambda_3 = lambda_3): out = np.zeros(n_count_data) # number of data points out[:tau_1] = lambda_1 # lambda before tau is lambda_1 out[tau_1:tau_2] = lambda_2 out[tau_2:] = lambda_3 # lambda after (and including) tau # is lambda_2 return out observa-tion = pm.Poisson("obs", lambda_, value=count_data,observed=True) mod-el = pm.Model([observation, lambda_1, lambda_2, lambda_3, tau_1, tau_2]) mcmc = pm.MCMC(model) mcmc.sample(40000, 10000)

图1.6.1展示了5个未知数的后验。我们可以看到模型的转折点大致在第45天和第47天的时候取得。对此你怎么认为呢？我们的模型是否对数据过拟合呢？

确实，我们都可能对数据中有多少个转折点抱有疑惑的态度。例如，我就认为一个转折点好过两个转折点，两个转折点好过三个转折点，以此类推。这意味着对于应该有多少个转折点可以设置一个先验分布并让模型自己做决定！在对模型进行调整之后，答案是肯定的，一个转折点确实比较适合。代码在本章就不展示了，这里我只是想介绍一种思想：用怀疑数据那样的眼光审视我们的模型。

相关资源：敏捷开发V1.0.pptx