PyMC3中自定义多元Dirichlet先验的实现问题

让我来帮你解决在PyMC3中实现Dirichlet线性变换自定义分布的问题——你遇到的错误主要是因为破坏了Theano的符号计算逻辑，还有对DensityDist的用法理解不对。咱们一步步来修正：

原代码的核心问题 #

1. 误用.eval()：你在log_trunc_dir里调用了.eval()，这会尝试在模型构建阶段直接计算符号张量的数值，但此时x是PyMC的潜在随机变量，还没有采样值，Theano无法完成这个计算，自然会抛出数据类型相关的错误。PyMC需要符号化的对数概率函数来构建计算图，才能进行后续的采样和推断。
2. observed参数用错：DensityDist的observed应该传入实际的观测数据数组，而不是模型中的另一个随机变量x。如果你的目标是定义由x变换得到的新分布，不需要设置observed；如果是基于观测数据推断x，应该把观测数据传给observed。

修正后的实现方案 #

要定义Dirichlet变量经过线性变换后的自定义分布，需要考虑概率密度变换的雅可比行列式——因为当你对随机变量做可逆变换时，概率密度会跟着变换，必须加入雅可比项来修正对数概率。

完整可运行代码

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

# 数据设置
n = 5
# Dirichlet的a参数对应n个类别，所以设置为n维
prior_params = np.ones(n)
mx = np.array([[0.25 , 0.5 , 0.75 , 1. ],
               [0.25 , 0.333, 0.25 , 0. ],
               [0.25 , 0.167, 0. , 0. ],
               [0.25 , 0. , 0. , 0. ]])
# 确保mx可逆，计算逆矩阵（如果mx不可逆，这个方法不适用，需要调整逻辑）
mx_inv = np.linalg.inv(mx)

# 自定义对数概率函数：返回Theano符号张量，不能用eval()
def log_transformed_dirichlet(q, mx_inv, prior_params):
    # 逆变换：从q还原出原始Dirichlet变量x
    x = tt.dot(mx_inv, q)
    # 计算x的Dirichlet对数概率
    dirichlet_logp = pm.Dirichlet.dist(a=prior_params).logp(x)
    # 计算线性变换的雅可比行列式对数绝对值（修正概率密度）
    jacobian_logdet = tt.log(tt.abs_(tt.nlinalg.det(mx)))
    # 概率密度变换公式：logp(q) = logp(x) - log|det(雅可比矩阵)|
    return dirichlet_logp - jacobian_logdet

# 构建模型
with pm.Model() as simple_model:
    # 定义自定义分布的潜在变量q，shape为n-1维（对应n类的单纯形）
    q = pm.DensityDist('q', log_transformed_dirichlet,
                       params={'mx_inv': mx_inv, 'prior_params': prior_params},
                       shape=(n-1,))
    
    # 用NUTS采样
    trace = pm.sample(2000, tune=1000, cores=2)

# 查看采样结果
pm.summary(trace)
pm.traceplot(trace)

针对不同需求的调整 #

如果你是想基于观测到的q数据推断潜在变量x，只需要修改模型部分：

# 假设有观测到的q数据
q_data = np.array([0.4, 0.3, 0.2, 0.1])

with pm.Model() as inference_model:
    # 把观测数据传给observed，shape自动匹配数据维度
    pm.DensityDist('q_obs', log_transformed_dirichlet,
                   observed=q_data,
                   params={'mx_inv': mx_inv, 'prior_params': prior_params})
    
    # 采样推断
    trace = pm.sample(2000, tune=1000, cores=2)

关键注意事项 #

• 必须确保变换矩阵mx是可逆的，否则无法通过逆变换还原出原始的Dirichlet变量x。如果mx不可逆，你需要额外添加约束条件来保证x落在单位单纯形内。
• 所有对数概率的计算都要保持符号化（用Theano张量操作），绝对不能在模型构建阶段调用.eval()或.numpy()这类方法强行转换为数值。

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