无法让 DiscreteMarkovChain 在 PyMC 中工作

Question

我正在尝试使用 PyMC 拟合两国政权切换模型。

完整模型如下。 模型的具体细节并不重要（如果你还是好奇的话，我正在尝试拟合一个几何布朗运动，其系数 μ 和 σ 是该体系的函数）。

我收到一个无法解释的模糊错误（也在下面提供）。 我的问题如下：我使用的是

DiscreteMarkovChain

还是有bug？有没有更好/更惯用的方式来表达这个模型？

with pm.Model() as model
    π0 = pm.Uniform("π01", lower=0., upper=10.)
    π1 = pm.Uniform("π10", lower=0., upper=10.)
    P = pm.math.stack([[1. - π0 * Δt, π0 * Δt], [π1 * Δt, 1. - π1 * Δt]])
    state = pymc_experimental.distributions.timeseries.DiscreteMarkovChain("state", P=P, steps=Δx.size - 1)

    μ0 = pm.Normal("μ0", mu=0, sigma=0.1)
    μ1 = pm.Normal("μ1", mu=0, sigma=0.1)
    
    σ0 = pm.HalfNormal("σ0", sigma=0.1)
    σ1 = pm.HalfNormal("σ1", sigma=0.1)
    
    μ = pm.math.switch(state, μ0, μ1)
    σ = pm.math.switch(state, σ0, σ1)
 
    ΔX = pm.Normal("ΔX", mu=(μ - σ**2 / 2) * Δt, sigma=σ * math.sqrt(Δt), observed=Δx)

Traceback (most recent call last):
  File "/home/parsiad/venv/lib/python3.12/site-packages/pymc/sampling/parallel.py", line 128, in run
    self._start_loop()
  File "/home/parsiad/venv/lib/python3.12/site-packages/pymc/sampling/parallel.py", line 180, in _start_loop
    point, stats = self._step_method.step(self._point)
                   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/parsiad/venv/lib/python3.12/site-packages/pymc/step_methods/compound.py", line 232, in step
    point, sts = method.step(point)
                 ^^^^^^^^^^^^^^^^^^
  File "/home/parsiad/venv/lib/python3.12/site-packages/pymc/step_methods/arraystep.py", line 101, in step
    apoint, stats = self.astep(q)
                    ^^^^^^^^^^^^^
  File "/home/parsiad/venv/lib/python3.12/site-packages/pymc/step_methods/metropolis.py", line 274, in astep
    accept_rate = self.delta_logp(q, q0d)
                  ^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/parsiad/venv/lib/python3.12/site-packages/pytensor/compile/function/types.py", line 972, in __call__
    raise_with_op(
  File "/home/parsiad/venv/lib/python3.12/site-packages/pytensor/link/utils.py", line 528, in raise_with_op
    raise exc_value.with_traceback(exc_trace)
  File "/home/parsiad/venv/lib/python3.12/site-packages/pytensor/compile/function/types.py", line 959, in __call__
    self.vm()
  File "/home/parsiad/venv/lib/python3.12/site-packages/pytensor/graph/op.py", line 524, in rval
    r = p(n, [x[0] for x in i], o)
        ^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/parsiad/venv/lib/python3.12/site-packages/pytensor/tensor/subtensor.py", line 2754, in perform
    rval = inputs[0].__getitem__(tuple(inputs[1:]))
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
IndexError: index 2 is out of bounds for axis 0 with size 2
Apply node that caused the error: AdvancedSubtensor(Join.0, Subtensor{start:stop}.0, Subtensor{start:}.0)
Toposort index: 32
Inputs types: [TensorType(float64, shape=(2, 2)), TensorType(int64, shape=(24113,)), TensorType(int64, shape=(24113,))]
Inputs shapes: [(2, 2), (24113,), (24113,)]
Inputs strides: [(16, 8), (8,), (8,)]
Inputs values: [array([[0.98880855, 0.01119145],
       [0.02356341, 0.97643659]]), 'not shown', 'not shown']
Inputs type_num: [12, 7, 7]
Outputs clients: [[CAReduce{Composite{(i0 + log(i1))}, axes=None}(AdvancedSubtensor.0)]]

我在 PyMC 文档中搜索了

DiscreteMarkovChain

的替代品。我还尝试调试我当前的代码。

Answer 1

我相信你应该这样做。主要思想是使用

state

对

BinaryGibbsMetropolis

变量进行采样。

import pymc as pm
from pymc_experimental.distributions import DiscreteMarkovChain
import math 
import numpy as np



with pm.Model() as model:

    Δx = np.random.randn(20)
    Δt =1

    π0 = pm.Uniform("π01", lower=0., upper=1.)
    π1 = pm.Uniform("π10", lower=0., upper=1.)
    s0 = pm.Bernoulli.dist(p=0.5)
    P = pm.Dirichlet("P", a=[π0 * Δt, π1 * Δt], size=(2,))
    state = DiscreteMarkovChain("state", P=P,init_dist=s0, steps=Δx.size - 1)

    μ0 = pm.Normal("μ0", mu=0, sigma=0.1)
    μ1 = pm.Normal("μ1", mu=0, sigma=0.1)
    
    σ0 = pm.HalfNormal("σ0", sigma=0.1)
    σ1 = pm.HalfNormal("σ1", sigma=0.1)
    
    μ = pm.math.switch(state[-1], μ0, μ1)
    σ = pm.math.switch(state[-1], σ0, σ1)
 
    ΔX = pm.Normal("ΔX", mu=(μ - σ**2 / 2) * Δt, sigma=σ * math.sqrt(Δt), observed=Δx)


    idata = pm.sample(
        step=[
            pm.BinaryGibbsMetropolis([state]),
        ],
        draws=5_000,
    )