PyMC3估计GIG分布参数：scipy.special函数接收RV输入问题

我尝试用PyMC3估计广义逆高斯分布（GIG）的参数，其中用到了scipy.special里的贝塞尔函数。但贝塞尔函数要求输入是数组，而alpha、beta、gamma都是PyMC3的随机变量类，不知道怎么让Scipy的函数接收PyMC3随机变量作为输入。运行代码时出现如下错误：

运行代码 #

import pymc3 as pm
from scipy.special import hankel
import numpy as np


def gig(x, a, b, p):
    # c = p, is the order
    kp = special.hankel1e(p, x)
    y1 = ((a / b) ** (p / 2)) / (2 * kp * np.sqrt(a * b))
    y2 = (x ** (p - 1)) * np.exp(-(a * x + b / x) / 2)
    y = y1 * y2
    return y

with pm.Model() as gig_model:
    alpha = pm.Gamma('alpha', alpha=0.5, beta=2)
    beta = pm.Gamma('beta', alpha=0.5, beta=2)
    gamma = pm.Gamma('gamma', alpha=0.5, beta=2)
    
    def giglogp(x):
        lp = np.log(GIG(x, alpha, beta, gamma))
        return lp
        
    # likelihood
    Like = pm.DensityDist('likelihood', giglogp, observed=dt)

报错信息 #

TypeError: ufunc 'hankel1e' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''

解决方案 #

核心原因 #

Scipy的特殊函数（如hankel1e）仅支持NumPy数组输入，无法直接处理PyMC3的Theano张量（PyMC3底层依赖Theano构建计算图）。必须将Scipy函数封装为Theano可识别的操作，同时把所有数值计算替换为Theano张量操作。

修正步骤及代码 #

1. 自定义Theano操作（Op）封装Scipy的hankel1e，让其支持Theano张量输入，并提供梯度计算（MCMC采样需要梯度）。
2. 将原代码中的NumPy函数（如np.sqrt、np.exp）替换为Theano张量函数（tt.sqrt、tt.exp），确保整个计算图兼容PyMC3。
3. 修正原代码中的导入错误和函数调用错误。

修改后的完整代码：

import pymc3 as pm
import numpy as np
from scipy.special import hankel1e
import theano.tensor as tt
from theano.gof import Op, Apply

# 自定义TheanoOp封装hankel1e函数
class Hankel1eOp(Op):
    itypes = [tt.dscalar, tt.dscalar]  # 输入：贝塞尔阶数p、自变量x
    otypes = [tt.dscalar]

    def perform(self, node, inputs, outputs):
        # 执行Scipy的hankel1e计算
        p, x = inputs
        outputs[0][0] = np.array(hankel1e(p, x))

    def grad(self, inputs, g):
        # 数值梯度实现（若有解析导数可替换为精确形式）
        p, x = inputs
        eps = 1e-6
        # 对x求导
        d_hankel_dx = (hankel1e(p, x + eps) - hankel1e(p, x - eps)) / (2 * eps)
        # 对p求导
        d_hankel_dp = (hankel1e(p + eps, x) - hankel1e(p - eps, x)) / (2 * eps)
        return [g[0] * d_hankel_dp, g[0] * d_hankel_dx]

# 实例化自定义操作
hankel1e_tt = Hankel1eOp()

def gig(x, a, b, p):
    kp = hankel1e_tt(p, x)
    y1 = ((a / b) ** (p / 2)) / (2 * kp * tt.sqrt(a * b))
    y2 = (x ** (p - 1)) * tt.exp(-(a * x + b / x) / 2)
    y = y1 * y2
    return y

with pm.Model() as gig_model:
    alpha = pm.Gamma('alpha', alpha=0.5, beta=2)
    beta = pm.Gamma('beta', alpha=0.5, beta=2)
    gamma = pm.Gamma('gamma', alpha=0.5, beta=2)
    
    def giglogp(x):
        lp = tt.log(gig(x, alpha, beta, gamma))
        return lp
        
    # 替换为你的实际观测数据
    dt = np.random.rand(100)  # 示例数据
    Like = pm.DensityDist('likelihood', giglogp, observed=dt)

注意事项 #

• 如果使用MAP估计而非MCMC采样，可省略grad方法的实现，但MCMC（如NUTS）必须有梯度计算才能正常运行。
• 数值梯度的精度和效率略低于解析梯度，若能找到hankel1e的解析导数公式，建议替换grad方法中的数值计算。

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