# Stochastic Matrix Performance

**URL:** <https://discourse.pymc.io/t/stochastic-matrix-performance/7228>\
**Category:** Questions\
**Created:** [April 12, 2021, 9:00pm UTC](https://discourse.pymc.io/t/stochastic-matrix-performance/7228 "2021-04-12T21:00:15Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![justinrporter](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/justinrporter/32/3924_2.png) [@justinrporter](https://discourse.pymc.io/u/justinrporter)\
**Post date:** [April 12, 2021, 9:00pm UTC](https://discourse.pymc.io/t/stochastic-matrix-performance/7228/1 "2021-04-12T21:00:15Z")

</div>

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](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/7/77e914f371c549a28dda4fb0632b38810d972673.png)

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!

---

<div class="post-metadata">

**Author:** ![justinrporter](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/justinrporter/32/3924_2.png) [@justinrporter](https://discourse.pymc.io/u/justinrporter)\
**Post date:** [April 13, 2021, 2:42am UTC](https://discourse.pymc.io/t/stochastic-matrix-performance/7228/2 "2021-04-13T02:42:42Z")

</div>

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

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

```python

    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.

---

<div class="post-metadata">

**Author:** ![twiecki](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/twiecki/32/6930_2.png) [@twiecki](https://discourse.pymc.io/u/twiecki)\
**Post date:** [April 14, 2021, 10:26am UTC](https://discourse.pymc.io/t/stochastic-matrix-performance/7228/3 "2021-04-14T10:26:29Z")

</div>

This is an ODE, no?

In that case you can use SunODE with PyMC3: [GitHub - aseyboldt/sunode: Solve ODEs fast, with support for PyMC3](https://github.com/aseyboldt/sunode)

---

<div class="post-metadata">

**Author:** ![justinrporter](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/justinrporter/32/3924_2.png) [@justinrporter](https://discourse.pymc.io/u/justinrporter)\
**Post date:** [April 14, 2021, 3:22pm UTC](https://discourse.pymc.io/t/stochastic-matrix-performance/7228/4 "2021-04-14T15:22:55Z")

</div>

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!