PyMC3自定义似然函数及MAP估计遇数值集成报错求助

我看到你在PyMC3中尝试自定义似然函数并求解MAP估计时遇到了不少问题，包括Theano Op的测试值错误，以及似然函数的实现逻辑不符合PyMC3的张量计算范式。让我们一步步拆解问题并修正你的代码：

关键问题分析 #

你的代码存在几个核心问题：

1. 变量覆盖与作用域混乱：在模型内部重新定义了beta = tt.scalar('beta')和theta = tt.scalar('theta')，这直接覆盖了之前定义的PyMC3随机变量beta=pm.Uniform(...)，导致模型无法识别参数变量。
2. 似然函数实现错误：你预计算了所有beta和M组合的arr，然后在logp_func中返回min(arr)，这完全不是似然函数的正确写法——似然函数应该基于当前模型参数（M和beta）计算观测数据的对数概率，而不是预计算所有可能值的最小值。
3. 自定义Op缺少测试值：你的integrateOut Op中的theta没有设置默认值，Theano需要测试值来验证Op的计算流程，因此抛出了缺失默认值的错误。
4. 数据扰动的位置错误：在模型内部循环生成vnew和rnew的随机扰动会导致每次模型编译时数据都变化，应该提前在模型外处理好带扰动的数据。

修正后的代码实现 #

下面是修复后的完整代码，我会在关键部分添加注释：

import numpy as np
import pandas as pd
import pymc3 as pm
import theano.tensor as tt
from scipy.integrate import quad
import random
import decimal as dm

# ----------------------
# 1. 预处理数据（包括扰动）
# ----------------------
data_ = np.array(pd.read_excel('aaa.xlsx', header=None))
gamma = 3.77
G = 4.302 * 10**-6
rmin = 3.0
R = 95.7
vr = data_[:,1]
r = data_[:,0]

# 提前生成带扰动的数据（模型外处理，避免每次编译模型都重新生成）
np.random.seed(1)  # 设置随机种子保证可复现
vnew = vr + (0.05 * vr * np.random.uniform(0.1, 2.0, size=len(vr)))
rnew = r + (0.05 * r * np.random.uniform(0.1, 2.0, size=len(r)))
vn = np.array(vnew)
rn = np.array(rnew)

# ----------------------
# 2. 修复自定义积分Op（添加测试值）
# ----------------------
class integrateOut(theano.Op):
    def __init__(self, f, t, t0, tf, *args, **kwargs):
        super(integrateOut, self).__init__()
        self.f = f
        self.t = t
        self.t0 = t0
        self.tf = tf
        # 给积分变量t设置默认测试值，解决Theano的测试值错误
        self.t.default = 0.0

    def make_node(self, *inputs):
        self.fvars = list(inputs)
        # 确保输入都是Theano张量类型
        inputs = [tt.as_tensor_variable(inp) for inp in inputs]
        try:
            self.gradF = tt.grad(self.f, self.fvars)
        except:
            self.gradF = None
        return theano.Apply(self, self.fvars, [tt.dscalar().type()])

    def perform(self, node, inputs, output_storage):
        args = tuple(inputs)
        # 编译Theano函数，计算被积函数
        f = theano.function([self.t] + self.fvars, self.f)
        # 调用scipy的quad进行数值积分
        output_storage[0][0] = quad(f, self.t0, self.tf, args=args)[0]

    def grad(self, inputs, grads):
        # 实现梯度的积分（利用积分的导数等于导数的积分）
        if self.gradF is None:
            return [tt.zeros_like(inp) for inp in inputs]
        return [integrateOut(g, self.t, self.t0, self.tf)(*inputs) * grads[0] 
                for g in self.gradF]

# ----------------------
# 3. 定义PyMC3模型与正确的自定义似然
# ----------------------
basic_model = pm.Model()

with basic_model:
    # 定义模型参数（注意：不要在模型内部重新定义beta为tt.scalar！）
    M = pm.Uniform('M', lower=0.5*10**12, upper=3.50*10**12, transform='interval')
    beta = pm.Uniform('beta', lower=2.001, upper=2.999, transform='interval')

    # 定义被积函数（使用Theano张量操作，这样PyMC3能自动求导）
    theta = tt.scalar('theta')
    z = tt.cos(theta)**(2 * ((gamma / (beta - 2)) - 3/2) + 3)
    # 创建积分Op实例，传入beta作为参数
    intZ = integrateOut(z, theta, -np.pi/2, np.pi/2)(beta)

    # 基于参数M和beta计算似然的对数概率
    q = (gamma / (beta - 2)) - 3/2
    B = (G * M * 10**12) / ((beta - 2) * (R**(3 - beta)))
    K = (gamma - 3) / ((rmin**(3 - gamma)) * intZ * tt.sqrt(2 * B))
    
    # 计算每个数据点的对数概率，然后求和
    log_likelihood = -tt.log(K * ((1 - (1/(2*B)) * ((vn**2) * (rn**(beta - 2))))**(q + 1)) * (rn**(1 - gamma + (beta/2))))
    total_logp = tt.sum(log_likelihood)

    # 使用DensityDist定义自定义似然，传入观测数据
    pm.DensityDist("custom_likelihood", lambda v, r: total_logp, observed={"v": vn, "r": rn})

    # 计算MAP估计
    map_estimate = pm.find_MAP(model=basic_model)
    print("MAP Estimates:")
    print(map_estimate)

    # （可选）进行MCMC采样验证
    step = pm.Metropolis()
    trace = pm.sample(10000, step=step, start=map_estimate, random_seed=1, progressbar=True)
    print("\nMCMC Summary:")
    print(pm.summary(trace))

关键修改说明 #

1. 数据预处理：把随机扰动移到模型外，并设置随机种子保证结果可复现，避免模型编译时数据变化。
2. 自定义Op修复：给theta添加default=0.0，解决Theano的测试值缺失错误；同时在make_node中确保输入都是Theano张量类型。
3. 变量作用域修正：不再在模型内部重新定义beta为Theano标量，直接使用PyMC3的随机变量beta参与张量计算。
4. 似然函数重构：完全用Theano张量操作定义似然，直接基于模型参数M和beta计算对数概率总和，这样PyMC3能自动处理求导和参数优化。
5. DensityDist的正确使用：DensityDist的logp函数需要接收观测数据和参数，这里我们直接返回基于参数计算的total_logp，确保PyMC3能关联参数与观测数据。

这样修改后，你的模型应该能正常计算MAP估计，并且可以进行MCMC采样验证结果。

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