如何在pymcmcstat中使用自定义对数似然函数？参数解包错误排查

我想用pymcmcstat做参数估计，但不知道怎么让平方和函数适配我的问题。我已经算出两个对数正态分布的对数似然值，想把它用到模型里。当前代码如下：

from pymcmcstat.MCMC import MCMC

mcstat = MCMC()  # 初始化MCMC
x = t
y = I*100000
mcstat.data.add_data_set(x,y) # 用种群数量数据创建数据集
def ssfun(q, data): # 不确定该用平方和函数还是对数似然函数？
    m, b = q # 斜率和截距
    x = data.xdata[0]
    y = data.ydata[0]
    # 计算模型预测值
    ymodel = m*x + b
    res = ymodel - y
    return (res ** 2).sum(axis=0)

mcstat.model_settings.define_model_settings(sos_function=ssfun) # 将函数传入模型
mcstat.simulation_options.define_simulation_options(nsimu=10.0e3) # 设置模拟次数
mcstat.parameters.add_model_parameter(name='beta', theta0=8, minimum = 0, maximum = 20) # 先验参数
mcstat.parameters.add_model_parameter(name='gamma',theta0=0.4, minimum = 0, maximum=2) # 先验参数
mcstat.run_simulation() # 运行MCMC

其中t和I都是长度为500的数组，分别存储时间和对应种群数量的数值。

运行后出现如下错误：

File "<ipython-input-36-b8790f50839a>", line 2, in ssfun
    m, b = q # slope and offset

ValueError: too many values to unpack (expected 2)

我猜测是平方和函数的问题，但不知道怎么在这个场景下使用自定义的对数似然函数？

--- 补充内容 ---
我构建的对数似然函数是将两个对数正态分布的概率密度函数的对数相加得到的：

muPrev, sigmaPrev = 400, 40 # 对应种群数量I的参数
muInc, sigmaInc = 300, 30 # 对应cInc的参数
n = 10000

# 生成对数正态分布样本（已注释）
# logPrev = np.random.lognormal(np.log((muPrev**2) / (muPrev**2 + sigmaPrev**2)**0.5), (np.log(1 + (sigmaPrev**2 / muPrev**2)))**0.5, n)
# logInc = np.random.lognormal(np.log((muInc**2) / (muInc**2 + sigmaInc**2)**0.5), (np.log(1 + (sigmaInc**2 / muInc**2)))**0.5, n)

xPrev = 400 # 计算对数PDF时的输入值（对应I）
xInc = 300 # 计算对数PDF时的输入值（对应cInc）

# 计算对数正态分布在对数空间的均值
logmuPrev = np.log((muPrev**2) / (muPrev**2 + sigmaPrev**2)**0.5)
# 计算对数正态分布在对数空间的标准差
logsdPrev = (np.log(1 + (sigmaPrev**2 / muPrev**2)))**0.5

logmuInc = np.log((muInc**2) / (muInc**2 + sigmaInc**2)**0.5)
logsdInc = (np.log(1 + (sigmaInc**2 / muInc**2)))**0.5

# 计算两个分布的对数PDF值
L_prev = -0.5*((np.log(400) - logmuPrev) / logsdPrev)**2 - np.log(400* logsdPrev * (2*math.pi)**0.5)
L_inc = -0.5*((np.log(300) - logmuInc) / logsdInc)**2 - np.log(300* logsdInc * (2*math.pi)**0.5)

logsum = L_prev + L_inc # 对数似然求和
np.exp(logsum) # 指数转换得到似然值

1. 错误核心原因 #

你定义了2个模型参数（beta和gamma），但原代码存在两个关键矛盾：

• 你想用对数似然做参数估计，却错误绑定了平方和函数（sos_function）
• 原平方和函数里的线性模型（m*x + b）和你实际的对数正态似然模型完全无关，属于模型逻辑不匹配

2. 修正步骤与代码 #

要在pymcmcstat中使用自定义对数似然，需替换sos_function为loglikelihood_function，且函数需返回负对数似然值（因为库内部通过最小化该值实现似然最大化）。

修正后的完整代码

import numpy as np
import math
from pymcmcstat.MCMC import MCMC

mcstat = MCMC()
# 传入时间和种群数量数据
x = t
y = I * 100000
mcstat.data.add_data_set(x, y)

# 定义负对数似然函数
def custom_loglikelihood(q, data):
    # 解包参数：对应你定义的beta和gamma
    beta, gamma = q
    # 获取数据集
    time_data = data.xdata[0]
    pop_data = data.ydata[0]
    
    # --- 关键：用beta和gamma推导对数正态分布的参数 ---
    # 这里需要替换成你实际的模型逻辑，示例仅做演示
    # 假设muPrev由beta和时间决定，sigmaPrev由gamma决定
    muPrev = beta * np.exp(-0.01 * time_data)  # 示例模型
    sigmaPrev = gamma * 20
    # 同理处理第二个对数正态分布的参数（对应cInc）
    muInc = (beta / 2) * np.exp(-0.02 * time_data)
    sigmaInc = gamma * 15
    
    # 转换为对数空间的均值和标准差
    logmu_prev = np.log((muPrev**2) / np.sqrt(muPrev**2 + sigmaPrev**2))
    logsd_prev = np.sqrt(np.log(1 + (sigmaPrev**2 / muPrev**2)))
    logmu_inc = np.log((muInc**2) / np.sqrt(muInc**2 + sigmaInc**2))
    logsd_inc = np.sqrt(np.log(1 + (sigmaInc**2 / muInc**2)))
    
    # 对每个数据点计算对数似然并求和
    # 第一个分布的对数似然（对应种群数量数据）
    ll_prev = -0.5 * ((np.log(pop_data) - logmu_prev) / logsd_prev)**2 - np.log(pop_data * logsd_prev * np.sqrt(2 * math.pi))
    # 第二个分布的对数似然（需对应你的cInc数据，这里假设你有对应数据，需自行调整）
    # 示例：如果cInc是另一个数据集，需在add_data_set时传入，然后通过data.ydata[1]获取
    ll_inc = -0.5 * ((np.log(...) - logmu_inc) / logsd_inc)**2 - np.log(... * logsd_inc * np.sqrt(2 * math.pi))
    
    # 总对数似然
    total_ll = np.sum(ll_prev + ll_inc)
    # 返回负对数似然
    return -total_ll

# 绑定对数似然函数到模型
mcstat.model_settings.define_model_settings(loglikelihood_function=custom_loglikelihood)

# 设置模拟参数
mcstat.simulation_options.define_simulation_options(nsimu=10000)

# 添加参数先验
mcstat.parameters.add_model_parameter(name='beta', theta0=8, minimum=0, maximum=20)
mcstat.parameters.add_model_parameter(name='gamma', theta0=0.4, minimum=0, maximum=2)

# 运行MCMC
mcstat.run_simulation()

3. 重要注意事项 #

• 模型逻辑适配：你需要根据实际业务逻辑，补充beta/gamma到对数正态参数的映射关系，示例中的推导仅为演示。
• 多数据集处理：如果第二个对数正态分布对应另一组数据（如cInc），需通过mcstat.data.add_data_set添加第二个数据集，然后在似然函数中通过data.ydata[1]获取。
• 负对数似然要求：pymcmcstat内部使用优化算法最小化目标函数，因此必须返回负的对数似然值，才能等价于最大化似然。

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