frequentism-and-bayesianism-chs-iv
frequentism-and-bayesianism-chs-iv
频率主义与贝叶斯主义 IV：Python的贝叶斯工具
这个notebook出自Pythonic Perambulations的博文。The content is BSD licensed.
我之前花了一堆时间来分享这两种思想。
- 频率主义与贝叶斯主义 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)；这个相关性可以当作冗余参数
数据
先来定义一些数据
%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)
plt.plot(xdata, ydata, 'ok')
plt.xlabel('x')
plt.ylabel('y');
数据明显具有相关性，我们假设我们不知道误差。让我们构建一条直线来拟合。
模型：斜率，截距和未知相关系数
由贝叶斯定义可得
其中，D表示观察数据，θ表示模型
假设线性模型有斜率β，y轴截距α：
假设y轴坐标具有正态分布误差， 那么模型下的任何数据点的概率为：
其中，σ是未知测量误差，为冗余参数。
所以i的似然估计的乘积为：
先验条件
这里我们比之前的更仔细地选择先验条件。我们可以简单的认为α，β，和σ都是flat priors(扁平先验)，但是我们必须记住扁平先验并非都无信息先验(uninformative priors)!Jeffreys先验可能是更好的选择，用对称的最大熵(symmetry and/or maximum entropy)来选择最大化的无信息先验(maximally noninformative priors)。虽然这种做法被频率论认为太主观，但是这么做是可以从信息理论中找到理论依据的。
为什么flat prior可能产生坏的选择？看看斜率就明白了。让我们来演示一下斜率从0-10的直线变化情况；
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)
斜率和截距的先验
假设模型为
那么构建参数空间(parameter-space)，其概率元素为P(α,β) dα dβ
因为x和y是对称的，我们就可以进行参数交互
其概率元素为Q(α′,β′)dα′dβ′，由此可得
通过雅可比变换(the Jacobian of the transformation)，可得
为了保持对称性，需要保证变量的改变不会影响先验概率，因此有：
此函数满足：
即α服从均匀分布，β服从参数为sinθ的均匀分布，其中，θ=tan−1β。
你可能会奇怪，斜率的分布服从参数为sinθ的均匀分布，而非θ。sinθ可以被认为是源自截距。如果把变量α改为α⊥=αcosθ，那么均匀分布就变为(α⊥, θ)。在用PyStan求解时我们会用这个。
σ的先验
同理，我们希望σ的先验与问题的规模变化无关(如，改变点的数据)。因此其概率必须满足
方程等价于
这就是Harold Jeffreys提出的Jeffreys先验。
把先验放一起
把它们放到一起，我们用这些对称性参数推导出模型的最小无信息先验：
综上所述，你可能用扁平先验做参数变换(α,β,σ)→(α⊥,θ,logσ)，来解决这个问题， 但我认为对称/最大熵(symmetry/maximum entropy)方法更简洁明了——而且让我们能看到三个Python包演示这个先验的内涵。
用Python解决问题
现在有了数据，似然估计和先验，让我们用Python的emcee，PyMC和PyStan来演示。首先，让我们做一些辅助工作来可视化数据：
# 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