# Stochastic Matrix Performance

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:

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.

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):
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)
``````

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?