PyMC3拟合双峰分布报TypeError: Unsupported dtype for TensorType: object及收敛问题

初始报错原因 #

初始代码中使用scipy.stats.norm.pdf计算似然，该函数是基于numpy的数值计算接口，无法识别PyMC3内部的Theano张量类型参数，因此触发类型错误。编写自定义对数似然时必须全量使用pm.math下的张量运算函数，不能混合numpy/scipy的数值运算接口。

修改后参数不收敛的核心原因 #

自定义的正态分布概率密度公式存在严重错误：
正确的正态分布PDF公式为：
$$p(y|\mu,\sigma) = \frac{1}{\sigma\sqrt{2\pi}} \exp\left( -\frac{(y-\mu)^2}{2\sigma2} \right)$$
你的代码存在3个问题：

• 指数部分缺少负号，导致似然计算逻辑完全颠倒
• 指数分母缺少系数2
• 归一化项缺少$\sqrt{2\pi}$系数（该问题不影响极值位置，但前两个错误直接导致似然计算完全失效）

方案1：修复自定义对数似然 #

import pymc3 as pm
import numpy as np
import theano.tensor as tt

# 生成模拟数据
mu1_true = 1
sigma1_true = 0.1
mu2_true = 2
sigma2_true = 0.2
w1_true = 2/3
n = 7500
n1 = int(n*w1_true)
n2 = n - n1
y = np.concatenate((np.random.normal(mu1_true, sigma1_true, n1), 
                    np.random.normal(mu2_true, sigma2_true, n2)))

# 修正后的正态分布PDF张量运算
def normpdf(y, mu, sigma):
    return (1 / (sigma * tt.sqrt(2 * np.pi))) * tt.exp(-tt.pow(y - mu, 2) / (2 * tt.pow(sigma, 2)))

def logp(mu, sigma, w, y):
    comp1 = w[0] * normpdf(y, mu[0], sigma[0])
    comp2 = w[1] * normpdf(y, mu[1], sigma[1])
    return tt.sum(tt.log(comp1 + comp2))

with pm.Model() as model:
    # 给mu加有序约束，避免标签交换问题
    mu = pm.Normal("mu", mu=np.array([0.5, 1.5]), sigma=1, shape=2, 
                  transform=pm.distributions.transforms.ordered)
    sigma = pm.HalfNormal("sigma", sigma=1, shape=2)
    w = pm.Dirichlet("w", a=np.array([1, 1]), shape=2)
    
    likelihood = pm.DensityDist('Likelihood', logp, observed=dict(mu=mu, sigma=sigma, w=w, y=y))
    # 优先用MCMC采样，比MAP点估计更稳定，避免陷入局部最优
    trace = pm.sample(2000, tune=1000, cores=2, return_inferencedata=True)

# 查看参数均值
print(trace.posterior.mean())

方案2：使用PyMC3内置Mixture分布（更稳定，无需自定义logp） #

PyMC3已经内置了混合分布实现，无需手动编写似然，避免公式写错的问题：

with pm.Model() as model:
    mu = pm.Normal("mu", mu=np.array([0.5, 1.5]), sigma=1, shape=2,
                  transform=pm.distributions.transforms.ordered)
    sigma = pm.HalfNormal("sigma", sigma=1, shape=2)
    w = pm.Dirichlet("w", a=np.array([1, 1]), shape=2)
    
    # 定义混合分量
    comps = [pm.Normal.dist(mu=mu[i], sigma=sigma[i]) for i in range(2)]
    # 定义混合分布作为似然
    likelihood = pm.Mixture('Likelihood', w=w, comp_dists=comps, observed=y)
    
    trace = pm.sample(2000, tune=1000, cores=2, return_inferencedata=True)

print(trace.posterior.mean())

运行后得到的参数均值会非常接近真实值：mu接近[1,2]，sigma接近[0.1,0.2]，w接近[0.666, 0.333]。

内容的提问来源于stack exchange，提问作者Tarik