如何使用 PyMC (v4.0.1) 和 statsmodels 来估计状态空间模型？

Question

我正在尝试将 PyMC (v4.0.1) 与 statsmodels 一起使用来估计状态空间模型。我正在关注这个使用 PyMC3 的示例：

https://www.statsmodels.org/v0.12.0/examples/notebooks/generated/statespace_sarimax_pymc3.html

该示例使用 pm.DensityDist() 函数，显然 API 已更改。 PyMC 使用 Aesara 而不是 Theano，我不知道这是否重要。

作为一个工作示例，这是我尝试模拟和估计 AR(1) 过程的尝试：

import numpy as np
import statsmodels.api as sm
import pymc as pm
import aesara.tensor as at
from scipy.signal import lfilter

# Generate artificial data
nobs = int(1e3)
true_phi = np.r_[0.5]
true_sigma = 0.5**0.5

np.random.seed(1234)
disturbances = np.random.normal(0, true_sigma, size=(nobs,))
endog = lfilter([1], np.r_[1, -true_phi], disturbances)

# Initialize model
mod = sm.tsa.statespace.SARIMAX(endog, order=(1, 0, 0))


# Helper functions copied. Do not know how they work
class Loglike(at.Op):

    itypes = [at.dvector] # expects a vector of parameter values when called
    otypes = [at.dscalar] # outputs a single scalar value (the log likelihood)

    def __init__(self, model):
        self.model = model
        self.score = Score(self.model)

    def perform(self, node, inputs, outputs):
        theta, = inputs  # contains the vector of parameters
        llf = self.model.loglike(theta)
        outputs[0][0] = np.array(llf) # output the log-likelihood

    def grad(self, inputs, g):
        # the method that calculates the gradients - it actually returns the
        # vector-Jacobian product - g[0] is a vector of parameter values
        theta, = inputs  # our parameters
        out = [g[0] * self.score(theta)]
        return out


class Score(at.Op):
    itypes = [at.dvector]
    otypes = [at.dvector]

    def __init__(self, model):
        self.model = model

    def perform(self, node, inputs, outputs):
        theta, = inputs
        outputs[0][0] = self.model.score(theta)
        
        
loglike = Loglike(mod)

# Set sampling params
ndraws = 3000  # number of draws from the distribution
nburn = 600   # number of "burn-in points" (which will be discarded)

# Sample from posterior
with pm.Model():
    # Priors
    arL1 = pm.Uniform('ar.L1', -0.99, 0.99)
    sigma2 = pm.InverseGamma('sigma2', 2, 4)

    # convert variables to tensor vectors
    theta = at.as_tensor_variable([arL1, sigma2])

    # use a DensityDist (use a lamdba function to "call" the Op)
    pm.DensityDist('likelihood', theta, logp = lambda v: loglike(v))

    # Draw samples
    trace = pm.sample(ndraws, tune=nburn, discard_tuned_samples=True, cores=1)

错误出在对pm.sample() 的调用中。

---> 74     trace = pm.sample(ndraws, tune=nburn, discard_tuned_samples=True, cores=1)

但错误消息表明该问题与似然函数有关：

TypeError: <lambda>() takes 1 positional argument but 2 were given

我已经尝试了很多东西，但我很茫然。我真的很感激任何建议。

Answer 1

通过替换解决了问题：

pm.DensityDist('likelihood', theta, logp = lambda v: loglike(v))

和

pm.Potential('likelihood', loglike(theta))

这是完整的工作代码：

import numpy as np
import statsmodels.api as sm
import pymc as pm
import aesara.tensor as at
from scipy.signal import lfilter

# Generate artificial data
nobs = int(1e3)
true_phi = np.r_[0.5]
true_sigma = 0.5**0.5

np.random.seed(1234)
disturbances = np.random.normal(0, true_sigma, size=(nobs,))
endog = lfilter([1], np.r_[1, -true_phi], disturbances)

# Initialize model
mod = sm.tsa.statespace.SARIMAX(endog, order=(1, 0, 0))


# Helper functions copied. Do not know how they work
class Loglike(at.Op):

    itypes = [at.dvector] # expects a vector of parameter values when called
    otypes = [at.dscalar] # outputs a single scalar value (the log likelihood)

    def __init__(self, model):
        self.model = model
        self.score = Score(self.model)

    def perform(self, node, inputs, outputs):
        theta, = inputs  # contains the vector of parameters
        llf = self.model.loglike(theta)
        outputs[0][0] = np.array(llf) # output the log-likelihood

    def grad(self, inputs, g):
        # the method that calculates the gradients - it actually returns the
        # vector-Jacobian product - g[0] is a vector of parameter values
        theta, = inputs  # our parameters
        out = [g[0] * self.score(theta)]
        return out


class Score(at.Op):
    itypes = [at.dvector]
    otypes = [at.dvector]

    def __init__(self, model):
        self.model = model

    def perform(self, node, inputs, outputs):
        theta, = inputs
        outputs[0][0] = self.model.score(theta)
        
        
loglike = Loglike(mod)

# Set sampling params
ndraws = 3000  # number of draws from the distribution
nburn = 600   # number of "burn-in points" (which will be discarded)

# Sample from posterior
with pm.Model():
    # Priors
    arL1 = pm.Uniform('ar.L1', -0.99, 0.99)
    sigma2 = pm.InverseGamma('sigma2', 2, 4)

    # convert variables to tensor vectors
    theta = at.as_tensor_variable([arL1, sigma2])

    # use a DensityDist (use a lamdba function to "call" the Op)
    pm.Potential('likelihood', loglike(theta))

    # Draw samples
    trace = pm.sample(ndraws, tune=nburn, discard_tuned_samples=True, cores=1)