問題

我想擬合一個簡單的 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)

然而結果非常奇怪。例如軌跡（第一個總體的分數）趨於 0 而不是收斂到 0.4，並且具有非常強的自相關性：

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

正如您在先前的選擇中看到的，我嘗試使用先前人口比例的 Beta 分佈來“幫助”MCMC 算法. 中心和西格瑪的先驗分佈也很合理（我認為）。

那麼這裡發生了什麼？是我做錯了什麼還是 MCMC 不適合這個問題？

我知道 MCMC 方法會更慢，但平凡的直方圖擬合似乎在解決人口問題方面表現得更好。

該問題是由 PyMC 為該模型抽取樣本的方式引起的。正如 PyMC 文檔的第 5.8.1 節所述，數組變量的所有元素都會一起更新。對於像這樣的小陣列center不是問題，但是對於像這樣的大陣列，category它會導致低接受率。您可以通過以下方式查看接受率

print mcmc.step_method_dict[category][0].ratio

文檔中建議的解決方案是使用標量值變量數組。此外，您需要配置一些提案分佈，因為默認選擇不好。這是對我有用的代碼：

import pymc as pm
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)
alpha  = pm.Beta('alpha', alpha=2, beta=3)
category = pm.Container([pm.Categorical("category%i" % i, [alpha, 1 - alpha]) 
                        for i in range(nsamples)])
observations = pm.Container([pm.Normal('samples_model%i' % i, 
                  mu=centers[category[i]], tau=1/(sigmas[category[i]]**2), 
                  value=samples[i], observed=True) for i in range(nsamples)])
model = pm.Model([observations, category, alpha, sigmas, centers])
mcmc = pm.MCMC(model)
# initialize in a good place to reduce the number of steps required
centers.value = [mu1_true, mu2_true]
# set a custom proposal for centers, since the default is bad
mcmc.use_step_method(pm.Metropolis, centers, proposal_sd=sig1_true/np.sqrt(nsamples))
# set a custom proposal for category, since the default is bad
for i in range(nsamples):
   mcmc.use_step_method(pm.DiscreteMetropolis, category[i], proposal_distribution='Prior')
mcmc.sample(100)  # beware sampling takes much longer now
# check the acceptance rates
print mcmc.step_method_dict[category[0]][0].ratio
print mcmc.step_method_dict[centers][0].ratio
print mcmc.step_method_dict[alpha][0].ratio

和選項在5.7.1 節proposal_sd中解釋。對於中心，我將建議設置為大致匹配後驗的標準差，由於數據量的原因，該標準差遠小於默認值。PyMC 確實會嘗試調整提案的寬度，但這僅在您的接受率足夠高的情況下才有效。對於，默認值會產生較差的結果（我不知道為什麼會這樣，但它肯定聽起來不像是二進制變量的明智提議）。proposal_distribution``category``proposal_distribution = 'Poisson'

引用自：https://stats.stackexchange.com/questions/120028