问题

我想拟合一个简单的 2-高斯混合总体的模型参数。鉴于围绕贝叶斯方法的所有炒作，我想了解对于这个问题，贝叶斯推理是否是比传统拟合方法更好的工具。

到目前为止，MCMC 在这个玩具示例中的表现非常糟糕，但也许我只是忽略了一些东西。所以让我们看看代码。

工具

我将使用 python (2.7) + scipy 堆栈、lmfit 0.8 和 PyMC 2.3。

可以在此处找到重现分析的笔记本

生成数据

首先让我们生成数据：

from scipy.stats import distributions

# Sample parameters
nsamples = 1000
mu1_true = 0.3
mu2_true = 0.55
sig1_true = 0.08
sig2_true = 0.12
a_true = 0.4

# Samples generation
np.random.seed(3)  # for repeatability
s1 = distributions.norm.rvs(mu1_true, sig1_true, size=round(a_true*nsamples))
s2 = distributions.norm.rvs(mu2_true, sig2_true, size=round((1-a_true)*nsamples))
samples = np.hstack([s1, s2])

的直方图samples如下所示：

一个“宽峰”，组件很难用肉眼发现。

经典方法：拟合直方图

让我们先尝试经典的方法。使用lmfit很容易定义一个 2-peaks 模型：

import lmfit

peak1 = lmfit.models.GaussianModel(prefix='p1_')
peak2 = lmfit.models.GaussianModel(prefix='p2_')
model = peak1 + peak2

model.set_param_hint('p1_center', value=0.2, min=-1, max=2)
model.set_param_hint('p2_center', value=0.5, min=-1, max=2)
model.set_param_hint('p1_sigma', value=0.1, min=0.01, max=0.3)
model.set_param_hint('p2_sigma', value=0.1, min=0.01, max=0.3)
model.set_param_hint('p1_amplitude', value=1, min=0.0, max=1)
model.set_param_hint('p2_amplitude', expr='1 - p1_amplitude')
name = '2-gaussians'

最后我们用单纯形算法拟合模型：

fit_res = model.fit(data, x=x_data, method='nelder')
print fit_res.fit_report()

结果如下图（红色虚线为拟合中心）：

即使问题有点难，在适当的初始值和约束条件下，模型也能收敛到相当合理的估计。

贝叶斯方法：MCMC

我在 PyMC 中以分层方式定义模型。centers并且sigmas是代表 2 个高斯的 2 个中心和 2 个 sigma 的超参数的先验分布。alpha是第一个总体的比例，这里的先验分布是 Beta。

分类变量在两个总体之间进行选择。据我了解，此变量需要与数据 ( samples) 的大小相同。

最后mu和tau是确定正态分布参数的确定性变量（它们取决于category变量，因此它们在两个总体的两个值之间随机切换）。

sigmas = pm.Normal('sigmas', mu=0.1, tau=1000, size=2)
centers = pm.Normal('centers', [0.3, 0.7], [1/(0.1)**2, 1/(0.1)**2], size=2)
#centers = pm.Uniform('centers', 0, 1, size=2)

alpha  = pm.Beta('alpha', alpha=2, beta=3)
category = pm.Categorical("category", [alpha, 1 - alpha], size=nsamples)

@pm.deterministic
def mu(category=category, centers=centers):
    return centers[category]

@pm.deterministic
def tau(category=category, sigmas=sigmas):
    return 1/(sigmas[category]**2)

observations = pm.Normal('samples_model', mu=mu, tau=tau, value=samples, observed=True)
model = pm.Model([observations, mu, tau, category, alpha, sigmas, centers])

然后我以相当长的迭代次数运行 MCMC（在我的机器上为 1e5，~60s）：

mcmc = pm.MCMC(model)
mcmc.sample(100000, 30000)

然而结果非常奇怪。例如$\alpha$轨迹（第一个总体的分数）趋于 0 而不是收敛到 0.4，并且具有非常强的自相关性：

高斯中心也不收敛。例如：

正如您在先前的选择中看到的，我尝试使用先前人口比例的 Beta 分布来“帮助”MCMC 算法$\alpha$. 中心和西格玛的先验分布也很合理（我认为）。

那么这里发生了什么？是我做错了什么还是 MCMC 不适合这个问题？

我知道 MCMC 方法会更慢，但平凡的直方图拟合似乎在解决人口问题方面表现得更好。