frequentism-and-bayesianism-chs-iv

frequentism-and-bayesianism-chs-iv

频率主义与贝叶斯主义 IV：Python的贝叶斯工具

 

这个notebook出自Pythonic Perambulations的博文。The content is BSD licensed.

 

这个系列共4个部分：中文版Part I Part II Part III Part IV，英文版Part I Part II Part III Part IV

 

我之前花了一堆时间来分享这两种思想。

• 频率主义与贝叶斯主义 I: 实用介绍 ，论述了频率主义和贝叶斯主义在科学计算方面的差异，在简单问题的处理方面两者基本等价，会得到相同的点估计(point estimates)结果。
• 频率主义 vs 贝叶斯主义 II：当结果不同时，进一步介绍两者的差异，尤其是处理冗余参数方面的差异。
• 频率主义 vs 贝叶斯主义 III：置信与可信，频率主义与科学，不能混为一谈，介绍了频率主义置信区间与贝叶斯主义可信范围的差异，指出在科学研究中频率主义只会回答错误的问题。

这一篇文章，我打算撇开哪些争论，谈谈Python实现贝叶斯研究的工具。 现代贝叶斯主义的核心是马尔可夫链蒙特卡罗算法，MCMC，样本后验分布(sample posterior distributions)的高效算法。

下面我就用Python的三个包来演示MCMC:

• emcee: the MCMC Hammer
• pymc: Bayesian Statistical Modeling in Python
• pystan: The Python Interface to Stan

我不会太关心它们的性能，也不会比较各自的API。这篇文章也不是教程；这三个包都有很完整的文档和教程。我要做的是比较一下三者的用法。用一个相对简单的例子，演示三个包的解法，希望对你有所帮助。

 

测试案例：线性拟合最优解

 

为了解决问题，我们做一个三参数模型来找数据的线性拟合解。参数是斜率，截距和相关性(scatter about the line)；这个相关性可以当作冗余参数

 

数据

 

先来定义一些数据

In [1]:

%matplotlib inline
import matplotlib.pyplot as plt import numpy as np np.random.seed(42) theta_true = (25, 0.5) xdata = 100 * np.random.random(20) ydata = theta_true[0] + theta_true[1] * xdata # add scatter to points xdata = np.random.normal(xdata, 10) ydata = np.random.normal(ydata, 10) 

In [2]:

plt.plot(xdata, ydata, 'ok') plt.xlabel('x') plt.ylabel('y'); 

 

 

数据明显具有相关性，我们假设我们不知道误差。让我们构建一条直线来拟合。

 

模型：斜率，截距和未知相关系数

 

由贝叶斯定义可得

P(θ | D)∝P(D | θ)P(θ)

其中，D表示观察数据，θ表示模型

 

假设线性模型有斜率β，y轴截距α：

y^(xi | α,β)=α+βxi

假设y轴坐标具有正态分布误差， 那么模型下的任何数据点的概率为：

P(xi,yi | α,β,σ)=12πσ2−−−−√exp[−[yi−y^(xi | α,β)]22σ2]

其中，σ是未知测量误差，为冗余参数。

所以i的似然估计的乘积为：

P({xi},{yi} | α,β,σ)∝(2πσ2)−N/2exp[−12σ2∑i−1N[yi−y^(xi | α,β)]2]

 

先验条件

 

这里我们比之前的更仔细地选择先验条件。我们可以简单的认为α，β，和σ都是flat priors(扁平先验)，但是我们必须记住扁平先验并非都无信息先验(uninformative priors)!Jeffreys先验可能是更好的选择，用对称的最大熵(symmetry and/or maximum entropy)来选择最大化的无信息先验(maximally noninformative priors)。虽然这种做法被频率论认为太主观，但是这么做是可以从信息理论中找到理论依据的。

为什么flat prior可能产生坏的选择？看看斜率就明白了。让我们来演示一下斜率从0-10的直线变化情况；

In [3]:

fig, ax = plt.subplots(subplot_kw=dict(aspect='equal')) x = np.linspace(-1, 1) for slope in np.arange(0, 10, 0.1): plt.plot(x, slope * x, '-k') ax.axis([-1, 1, -1, 1], aspect='equal'); 

 

 

这些线按0.1的斜率间隔进行分布，高斜率区域的线十分密集。因为是扁平先验，你肯定会认为这些斜率彼此无差异。由上面的密集区域显然可以看出，斜率的扁平先验满足那些很陡的斜率。斜率的扁平先验不是一个最小化的无信息先验(minimally informative prior)，可能是最终使你的结果有偏差(尽管有足够多的数据，结果可能几乎是0)。

你可能想动手做出一个更好的方案(可能在直线与x轴的夹角θ上用扁平先验)，但是我们有更严谨的方案。这个问题在贝叶斯文献中已经有了答案；我找到的最佳方案来自Jaynes的论文直线拟合：贝叶斯方法(Straight Line Fitting: A Bayesian Solution) (pdf)

 

斜率和截距的先验

 

假设模型为

y=α+βx

那么构建参数空间(parameter-space)，其概率元素为P(α,β) dα dβ

因为x和y是对称的，我们就可以进行参数交互

x=α′+β′y

其概率元素为Q(α′,β′)dα′dβ′，由此可得

(α′, β′)=(−β−1α, β−1).

通过雅可比变换(the Jacobian of the transformation)，可得

Q(α′,β′)=β3P(α,β).

为了保持对称性，需要保证变量的改变不会影响先验概率，因此有：

β3P(α,β)=P(−β−1α,β−1).

此函数满足：

P(α,β)∝(1+β2)−3/2.

即α服从均匀分布，β服从参数为sinθ的均匀分布，其中，θ=tan−1β。

你可能会奇怪，斜率的分布服从参数为sinθ的均匀分布，而非θ。sinθ可以被认为是源自截距。如果把变量α改为α⊥=αcosθ，那么均匀分布就变为(α⊥, θ)。在用PyStan求解时我们会用这个。

 

σ的先验

 

同理，我们希望σ的先验与问题的规模变化无关(如，改变点的数据)。因此其概率必须满足

P(σ)dσ=P(σ/c)dσ/c.

方程等价于

P(σ)∝1/σ.

这就是Harold Jeffreys提出的Jeffreys先验。

 

把先验放一起

 

把它们放到一起，我们用这些对称性参数推导出模型的最小无信息先验：

P(α,β,σ)∝1σ(1+β2)−3/2

综上所述，你可能用扁平先验做参数变换(α,β,σ)→(α⊥,θ,logσ)，来解决这个问题， 但我认为对称/最大熵(symmetry/maximum entropy)方法更简洁明了——而且让我们能看到三个Python包演示这个先验的内涵。

 

用Python解决问题

 

现在有了数据，似然估计和先验，让我们用Python的emcee，PyMC和PyStan来演示。首先，让我们做一些辅助工作来可视化数据：

In [4]:

# Create some convenience routines for plotting

def compute_sigma_level(trace1, trace2, nbins=20): """From a set of traces, bin by number of standard deviations""" L, xbins, ybins = np.histogram2d(trace1, trace2, nbins) L[L == 0] = 1E-16 logL = np.log(L) shape = L.shape L = L.ravel() # obtain the indices to sort and unsort the flattened array i_sort = np.argsort(L)[::-1] i_unsort = np.argsort(i_sort) L_cumsum = L[i_sort].cumsum() L_cumsum /= L_cumsum[-1] xbins = 0.5 * (xbins[1:] + xbins[:-1]) ybins = 0.5 * (ybins[1:] + ybins[:-1]) return xbins, ybins, L_cumsum[i_unsort].reshape(shape) def plot_MCMC_trace(ax, xdata, ydata, trace, scatter=False, **kwargs): """Plot traces and contours""" xbins, ybins, sigma = compute_sigma_level(trace[0], trace[1]) ax.contour(xbins, ybins, sigma.T, levels=[0.683, 0.955], **kwargs) if scatter: ax.plot(trace[0], trace[1], ',k', alpha=0.1) ax.set_xlabel(r'$\alpha$') ax.set_ylabel(r'$\beta$') def plot_MCMC_model(ax, xdata, ydata, trace):