Stochastic Matrix Performance

Questions
1

HI all,

I am making a model for multicompartment pharmacokinetics. I am interested in more complex models, but I am having trouble formulating even simple ones in a way that pymc3 is happy with. The basic model is like this:

model

At every timepoint, k01 x c0 mass flows from compartment 0 to compartment 1, k10 x c1 mass flows from compartment 1 to compartment 0, and u x c0 mass is excreted (destroyed).

I know that there are closed form solutions to this kind of model, but I am interested in more complex behavior, such as injections of mass at arbitrary times, where there are not closed-form solutions for c(t).

Therefore, my plan was to use a stochastic matrix to build a vector of concentrations over time, since c*K[t] = c[t+1]

Here’s the code:

tsteps = 10
tobs = np.arange(tsteps)
yobs = np.exp(-tobs/10) + np.random.normal(scale=0.1)    

with pm.Model() as model:

    theta = pm.Uniform('theta', lower=0, upper=1, shape=(3,))
    k01, k10, u = theta[0], theta[1], theta[2]

    K = tt.stack([tt.stack([1 - k01 - u,     k01]),
                  tt.stack([        k10, 1 - k10])])

    c = [[1, 0]] + [None]*(len(tobs)-1)

    for t in tobs[1:]:
        c[t] = tt.dot(c[t-1], K)

    c = tt.stack([c[t] for t in tobs])

    obs = pm.Normal("obs", mu=c[:, 0], sigma=0.1, observed=yobs)

    idata = pm.sample(1000)

This works well for tsteps < 20, but hangs before even starting to sample for tsteps = 30.

I strongly suspect that I am doing something (or several somethings) wrong here (for instance, the creation of lists that are then casted to theano tensors seems bad, for instance), but I’m not sure what to do differently.

Any advice would be appreciated!

2

Update: major progress has been made using theano.scan, for something like:

c, updates = theano.scan(
    fn=lambda c_prev, K: tt.dot(c_prev, K),
    outputs_info=c0,
    non_sequences=[K],
    n_steps=10
)

This results in a dramatic speedup over the above code.

I’d be curious if anyone knows how to generalize this to 3d tensors, so that I can fit multiple traces of c(t) simultaneously. I tried using tensordot but couldn’t seem to get it to come up with the right dimensions…

Right now, I just have a loop:


    cs = []
    for i, K in enumerate(Ks):
            c, updates = theano.scan(
                fn=lambda c_prev, K: tt.dot(c_prev, K),
                outputs_info=c0[i],
                non_sequences=[K],
                n_steps=Yobs.shape[1]
            )
        cs.append(c)

Again, advice is appreciated.

3

This is an ODE, no?

In that case you can use SunODE with PyMC3: GitHub - aseyboldt/sunode: Solve ODEs fast, with support for PyMC3

4

Yes, thank you for the suggestion! It definitely is, and I did try to do things with pyMC’s built-in ODE support, and it worked splendidly for the case above.

The issue is that I need to be able to model arbitrary u(t)—corresponding in this case to changes to the infusion rate of the drug—and I didn’t see how to do this. Do you think it would be possible?

Thanks!