pymc3: как реализовать детерминированную c операцию, требующую итерации? - PullRequest
Новая
       105

pymc3: как реализовать детерминированную c операцию, требующую итерации?

1 голос
/ 08 марта 2020

Ради этого вопроса я пытаюсь перенести пример pymc2 из этого сообщения в блоге , интерес представляет: 

## SI model

from pymc import *
from numpy import *

#observed data
T = 10
susceptible_data = array([999,997,996,994,993,992,990,989,986,984])
infected_data = array([1,2,5,6,7,8,9,11,13,15])

# stochastic priors
beta = Uniform('beta', 0., 40., value=1.)
gamma = Uniform('gamma', 0., 1., value=.001)
SI_0 = Uninformative('SI_0', value=[999., 1.])

# deterministic compartmental model
@deterministic
def SI(SI_0=SI_0, beta=beta, gamma=gamma):
    S = zeros(T)
    I = zeros(T)
    S[0] = SI_0[0]
    I[0] = SI_0[1]
    for i in range(1,T):
        S[i] = S[i-1] - 0.05*beta*S[i-1]*I[i-1]/(S[i-1]+I[i-1])
        I[i] = max(0., I[i-1] + 0.05*beta*S[i-1]*I[i-1]/(S[i-1]+I[i-1]) - gamma*I[i-1])
    return S, I
S = Lambda('S', lambda SI=SI: SI[0])
I = Lambda('I', lambda SI=SI: SI[1])

# data likelihood
A = Poisson('A', mu=S, value=susceptible_data, observed=True)
B = Poisson('B', mu=I, value=infected_data, observed=True)

Часть, в которой я застрял Это связано с тем, что автор использовал метод старой (обратной) разности с фиксированным размером шага 0,05 для имитации S и I на каждом шаге, очень разумный.

Однако я не могу понять способ сделать это в pymc3, так как указанное моделирование требует запуска (обратного) разностного приближения для массива, что, как мне кажется, не то, что я могу выяснить, как это сделать в pymc3.

Например, мой Первое приближение должно было сделать: 

import numpy as np
import pymc3 as pm
import theano

t = 10

susceptible_data = array([999,997,996,994,993,992,990,989,986,984])
infected_data = array([1,2,5,6,7,8,9,11,13,15])


def si_bdiff(si_0, beta, gamma):
    s = np.zeros(t, dtype=theano.config.floatX)
    f = np.zeros(t, dtype=theano.config.floatX)

    s[0] = si_0[0]
    f[0] = si_0[1]

    for i in range(1,t):
        s[i] = s[i-1] - 0.05 * beta * s[i-1] * f[i-1] / (s[i-1] + f[i-1])
        f[i] = max(0., f[i-1] + 0.05 * beta * s[i-1] * f[i-1] / (s[i-1] + f[i-1]) - gamma * f[i-1])
    return s, f


with pm.Model() as model:
    beta = pm.Uniform('beta', lower=0, upper=40)
    gamma = pm.Uniform('gamma', lower=0, upper=1)
    si_0 = pm.Flat('SI_0', shape=2)

    s, f = pm.Deterministic('si', si_bdiff(si_0, beta, gamma))

    # data likelihood
    A = pm.Poisson('A', mu=s, observed=susceptible_data)
    B = pm.Poisson('B', mu=f, observed=infected_data)

    trace1 = pm.sample(100000, tune=50000, cores=4, init='advi')

, что, конечно, с треском провалилось. Я возился с рядом аспектов этой проблемы, но вышеупомянутая версия содержит мою ключевую проблему в том, что я не знаю, как использовать pm.Deterministic таким образом, особенно для возможности повторения (что-то? ) и вернуть то, что в основном является парой массивов.

Спасибо!

РЕДАКТИРОВАТЬ

Это сводит меня с ума, я так близко, но ...

Вот смелость моего текущего кода: 

import numpy as np
import pymc3 as pm

import theano
import theano.tensor as tt
import theano.tests.unittest_tools

from scipy.optimize import fsolve


def dSdt(Scur, Icur, beta):
    N = Icur + Scur
    eqn = -beta * Scur * Icur / N
    return eqn

def dIdt(Scur, Icur, beta, gamma):
    N = Icur + Scur
    eqn = beta * Scur * Icur / N - gamma * Icur
    return eqn

def bkwd_euler(cur, prev, beta, gamma, h):
    Scur, Icur = cur
    Sprev, Iprev = prev
    N = Scur + Icur

    eqn0 = Scur - Sprev - h * dSdt(Scur, Icur, beta)
    eqn1 = Icur - Iprev - h * dIdt(Scur, Icur, beta, gamma)

    return eqn0, eqn1


T = infected_data.shape[0]


class SimulateDE(tt.Op):
    itypes = [tt.dscalar, tt.dscalar, tt.dscalar, tt.dscalar]
    otypes = [tt.dmatrix]

    tsteps = range(T)

    def perform(self, node, inputs, output):
        s0, i0, beta, gamma, = inputs[0], inputs[1], inputs[2], inputs[3]
        for i in self.tsteps:
            if i == 0:
                si = np.array((s0, i0)).reshape(1, -1)

            else:
                prev = si[-1]
                nxt = fsolve(bkwd_euler, (100, 100), (prev, beta, gamma, 1/7))
                si = np.vstack((si, nxt))

        output[0][0] = si

    def grad(self, inputs, g):
        # this is just a placeholder, because I cannot
        # figure out how grad() should work in this case
        return [g[0][0], g[0][1], g[0][2], g[0][3]]

Проблема, с которой я сейчас застрял, заключается в том, что я просто не могу понять, что должно быть в этом grad () случай, как я думаю, я просто смог бы пройти через g, как я; gibberi sh результат, но я не могу понять, почему он не должен работать. Но, увы, это не так: 

theano.config.compute_test_value = 'ignore'
theano.tests.unittest_tools.verify_grad(SimulateDE(), [np.array(100.), np.array(1.), np.array(10.), np.array(1.)])

 <snip>

NotImplementedError: ('input nd\nApply node that caused the error: InplaceDimShuffle{x,0}(<__main__.SimulateDE object at 0x7fdcec9594d0>.0)\nToposort index: 1\nInputs types: [TensorType(float64, vector)]\nInputs shapes: [(45, 2)]\nInputs strides: [(16, 8)]\nInputs values: [\'not shown\']\nOutputs clients: [[Elemwise{mul,no_inplace}(random_projection, InplaceDimShuffle{x,0}.0)]]\n\nBacktrace when the node is created(use Theano flag traceback.limit=N to make it longer):\n  File "/usr/local/lib/python3.7/dist-packages/IPython/core/interactiveshell.py", line 2886, in _run_cell\n    return runner(coro)\n  File "/usr/local/lib/python3.7/dist-packages/IPython/core/async_helpers.py", line 68, in _pseudo_sync_runner\n    coro.send(None)\n  File "/usr/local/lib/python3.7/dist-packages/IPython/core/interactiveshell.py", line 3063, in run_cell_async\n    interactivity=interactivity, compiler=compiler, result=result)\n  File "/usr/local/lib/python3.7/dist-packages/IPython/core/interactiveshell.py", line 3254, in run_ast_nodes\n    if (await self.run_code(code, result,  async_=asy)):\n  File "/usr/local/lib/python3.7/dist-packages/IPython/core/interactiveshell.py", line 3331, in run_code\n    exec(code_obj, self.user_global_ns, self.user_ns)\n  File "<ipython-input-108-a997c24f1e9a>", line 2, in <module>\n    theano.tests.unittest_tools.verify_grad(SimulateDE(), [np.array(100.), np.array(1.), np.array(10.), np.array(1.)])\n  File "/usr/local/lib/python3.7/dist-packages/theano/tests/unittest_tools.py", line 92, in verify_grad\n    T.verify_grad(op, pt, n_tests, rng, *args, **kwargs)\n  File "/usr/local/lib/python3.7/dist-packages/theano/gradient.py", line 1759, in verify_grad\n    cost = theano.tensor.sum(t_r * o_output)\n\nHINT: Use the Theano flag \'exception_verbosity=high\' for a debugprint and storage map footprint of this apply node.', '\nThe error happened with the following inputs:', [array(100.), array(1.), array(10.), array(1.)], '\nThe value of eps is:', None, '\nThe out_type is:', None)

Так что grad(), кажется, расстраивается из-за того, что вывод bkwd_euler является массивом (45, 2), может быть? Из того, что я могу понять из документов, вывод grad() должен быть 4-мерным вектором с 4 скалярными входами, поэтому не уверен, в чем проблема в этой точке.

Любая помощь по-прежнему высоко ценится!

...