我是概率编程的新手,正在研究 PyMC3 的高斯混合模型示例笔记本:
import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc3 as pm
import theano.tensor as tt
# simulate data from a known mixture distribution
np.random.seed(12345) # set random seed for reproducibility
k = 3
ndata = 500
spread = 5
centers = np.array([-spread, 0, spread])
# simulate data from mixture distribution
v = np.random.randint(0, k, ndata)
data = centers[v] + np.random.randn(ndata)
plt.hist(data);
# setup model
model = pm.Model()
with model:
# cluster sizes
p = pm.Dirichlet("p", a=np.array([1.0, 1.0, 1.0]), shape=k)
# ensure all clusters have some points
p_min_potential = pm.Potential("p_min_potential", tt.switch(tt.min(p) < 0.1, -np.inf, 0))
# cluster centers
means = pm.Normal("means", mu=[0, 0, 0], sigma=15, shape=k)
# break symmetry
order_means_potential = pm.Potential(
"order_means_potential",
tt.switch(means[1] - means[0] < 0, -np.inf, 0)
+ tt.switch(means[2] - means[1] < 0, -np.inf, 0),
)
# measurement error
sd = pm.Uniform("sd", lower=0, upper=20)
# latent cluster of each observation
category = pm.Categorical("category", p=p, shape=ndata)
# likelihood for each observed value
points = pm.Normal("obs", mu=means[category], sigma=sd, observed=data)
# fit model
with model:
step1 = pm.Metropolis(vars=[p, sd, means])
step2 = pm.ElemwiseCategorical(vars=[category], values=[0, 1, 2])
tr = pm.sample(10000, step=[step1, step2], tune=5000)
我正在为以下表达而苦恼:
# ensure all clusters have some points
p_min_potential = pm.Potential("p_min_potential", tt.switch(tt.min(p) < 0.1, -np.inf, 0))
和
# break symmetry
order_means_potential = pm.Potential(
"order_means_potential",
tt.switch(means[1] - means[0] < 0, -np.inf, 0)
+ tt.switch(means[2] - means[1] < 0, -np.inf, 0),
)
通过查看相关问题以及 PyMC3 和 Theano 的文档,我想我明白 pm.Potential()
是一种在采样期间设置模型中事件的对数似然性而不向其提供观察结果的方法,并且 tt.switch()
检查是否满足特定条件并相应地返回两个值之一。
因此 p_min_potential
通过设置一个事件的对数似然值来确保 p
中的所有值都大于 0.1,其中 p
中的一个值为负无穷大,同样 order_means_potential
确保 means
中的值与另一个不同,并且它们的顺序在采样期间保持不变.
我不明白的是:
这些表达式的结果如何反馈到模型中,因为 p_min_potential
和 order_means_potential
都不会作为任何其他表达式的输入出现?
我对 tt.switch()
的工作方式是否正确,因此如果条件 tt.min(p) < 0.1
被满足 -np.inf
作为事件对数似然返回,而在任何其他情况下返回 0?
任何帮助将不胜感激,我想了解这个例子是如何工作到我能够改变和扩展它的程度。具体来说,我想为两个或多个 beta 分布实现一个混合模型。