用 Python 进行贝叶斯模型建模(2)

444 查看

本系列:

第2节:模型检验

在这一节中,我们将用到两种技术,旨在回答:

  1. 模型和估计的参数是否对潜在数据有很好的模拟?
  2. 给定两个独立的模型,哪个对潜在数据模拟的更好?

模型检验I:后验估计检验

一种检验模型拟合的方法是后验估计检验。这种方法很直观,回顾上节中,我们通过收集 200,000 个 μ 的后验分布样本来对泊松分布的参数 μ进行估计,每个样本都被认为是可信的参数值。

后验预测检验需要从预测模型中产生新的数据。具体来说就是,我们已经估计了 200,000 个可信的泊松分布参数值μ,这意味着我们可以根据这些值建立 200,000 个泊松分布,然后从这些分布中随机采样。用公式表示为:
p(\tilde{y}|y) = \int p(\tilde{y}|\theta) f(\theta|y) d\theta

理论上,如果模型对潜在数据拟合的很好,那么产生的数据应该和原始观测数据近似。PyMC 提供一种方便的方式来从拟合模型中采样。你可能已经注意到了上面模型应用中新的一行代码:

这几乎和 y_est 一样,只不过没给观测数据赋值。PyMC 把它当做随机节点(和观测节点相反),当 MCMC 采样器运行时,它也从 y_est 中采集数据。

然后画出y_pred并和观测数据y_est比较。

选择正确的分布

我对上面的结果不是很满意。理想情况下,我希望后验预测分布和观测数据的分布一定程度上近似。直观上,如果我们正确估计了模型参数,那么我们应该可以从模型中采样得到类似的数据。结果很明显不是这样。

可能泊松分布不适合拟合这些数据。一种可选模型是负二项分布,特点和泊松分布很像,只是有两个参数(μ 和 α),使得方差和均值独立。回顾泊松分布只有一个参数 μ,既是均值,又是方差。

使用之前相同的数据集,继续对负二项分布的参数进行估计。同样地,使用均匀分布来估计 μ 和 α。模型可以表示为:

y_{j} \sim NegBinomial(\mu, \alpha)
\alpha = Exponential(0.2)
\mu = Uniform(0,100)

代码: