# How to model time-dependent variables in pymc3

**URL:** <https://discourse.pymc.io/t/how-to-model-time-dependent-variables-in-pymc3/516>\
**Category:** Questions\
**Created:** [November 1, 2017, 10:24am UTC](https://discourse.pymc.io/t/how-to-model-time-dependent-variables-in-pymc3/516 "2017-11-01T10:24:30Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![EdgarMCR](https://avatars.discourse-cdn.com/v4/letter/e/48db29/32.png) [@EdgarMCR](https://discourse.pymc.io/u/EdgarMCR)\
**Post date:** [November 1, 2017, 10:24am UTC](https://discourse.pymc.io/t/how-to-model-time-dependent-variables-in-pymc3/516/1 "2017-11-01T10:24:30Z")

</div>

Hello,

I seem unable to find an example for this anywhere: how can a time-dependent variable be modelled within the framework of pymc3?

Currently, I want to model an inhomogeneous Poisson process. For the sake of simplicity, lets say  
\lambda (t) = a \*t + b  
where a and b are parameters I wish to determine from some observational data.

My best attempt at implementing this so far is below. I am not sure how to include the time-dependence.

```
import numpy as np
import pymc3 as pm
import theano
import theano.tensor as tt
import matplotlib.pylab as plt

# Observation
a_actual = 1.3
b_actual = 2.0
t = np.arange(10)
obs = np.random.poisson(a_actual * t + b_actual)

plt.figure()
plt.plot(t, obs)
plt.show()

a1 = tt.scalar('a')
b1 = tt.vector('b')
y = tt.vector('y')

out = a1 * y + b1
func = theano.function([a1, b1, y], [out])

niter = 3000
# Model
with pm.Model() as model:
    a = pm.Uniform(name='a', lower=0, upper=10)
    b = pm.Uniform(name='b', lower=0, upper=10)

    #Expected value
    lam = pm.Deterministic('lam', func(a, b, t))
    count = pm.Poisson(mu=lam, name='count', observed=obs)
    
    #Run Monte Carlo
    trace = pm.sample(niter)

pm.plots.plot_posterior(trace=trace["a"])
pm.plots.autocorrplot(trace=trace, varnames=["a"]);
pm.plots.plot_posterior(trace=trace["b"])
pm.plots.autocorrplot(trace=trace, varnames=["b"]);
plt.show()

```

Thank you for any help or links to examples.

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [November 1, 2017, 11:00am UTC](https://discourse.pymc.io/t/how-to-model-time-dependent-variables-in-pymc3/516/2 "2017-11-01T11:00:17Z")

</div>

In your example above, the `lam` is a linear function of `t`, so you dont need to compile a theano function. Doing `lam = pm.Deterministic('lam', a*t + b)` is sufficient.

For more time-series model, you can check out these examples:  
[http://docs.pymc.io/notebooks/GLM-rolling-regression.html](http://docs.pymc.io/notebooks/GLM-rolling-regression.html)  
[http://docs.pymc.io/notebooks/stochastic\_volatility.html](http://docs.pymc.io/notebooks/stochastic_volatility.html)  
and the random walk related examples:  
[http://docs.pymc.io/notebooks/survival\_analysis.html?highlight=random%20walk](http://docs.pymc.io/notebooks/survival_analysis.html?highlight=random%20walk)

> <https://github.com/pymc-devs/pymc3/blob/master/docs/source/notebooks/MvGaussianRandomWalk_demo.ipynb>

---

<div class="post-metadata">

**Author:** ![etyurin](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/etyurin/32/340_2.png) [@etyurin](https://discourse.pymc.io/u/etyurin)\
**Post date:** [November 5, 2017, 6:52pm UTC](https://discourse.pymc.io/t/how-to-model-time-dependent-variables-in-pymc3/516/3 "2017-11-05T18:52:27Z")

</div>

@junpenglao - and how would you implement an integrated Random Walk in pymc3? My simplistic thinking is:

```
rw1 = pm.GaussianRandomWalk('rw1',sd=1., shape=lendata)
slope = pm.Deterministic('slope', pm.math.cumsum( rw1))

```

except there’s no`pm.math.cumsum().`… My question originated after reading this vignette:

[https://cran.r-project.org/web/packages/walker/vignettes/walker.html](https://cran.r-project.org/web/packages/walker/vignettes/walker.html)

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [November 5, 2017, 7:54pm UTC](https://discourse.pymc.io/t/how-to-model-time-dependent-variables-in-pymc3/516/4 "2017-11-05T19:54:15Z")

</div>

You can try [cumsum from theano](http://deeplearning.net/software/theano/library/tensor/extra_ops.html#theano.tensor.extra_ops.cumsum):  
`from theano.tensor.extra_ops import cumsum`

---

<div class="post-metadata">

**Author:** ![EdgarMCR](https://avatars.discourse-cdn.com/v4/letter/e/48db29/32.png) [@EdgarMCR](https://discourse.pymc.io/u/EdgarMCR)\
**Post date:** [November 6, 2017, 12:12pm UTC](https://discourse.pymc.io/t/how-to-model-time-dependent-variables-in-pymc3/516/5 "2017-11-06T12:12:18Z")

</div>

Thank you for your answers!

I was wondering if you know of an easy way to do incorporate a more complex function, specifically a skewed normal distribution.

I am trying to fit the function that describes a skewed normal distribution to some data because this function appears to fit well. While there is a SkewNormal in pymc3, it is used as a probability function and I am not sure how to use pymc3 distribution as a function.

Normally, in Python I would define the function as

```
def skew(x,e=0,w=1,a=0, mag=1):
     t = (x-e) / w
    return 2 * mag * scipy.stats.norm.pdf(t) * scipy.stats.norm.cdf(a*t)

```

However, this includes a definite integral when not using scipy function.

If you know a good way to:  
a) implement this function  
b) reuse the scipy functions within pymc3  
c) reuse the pymc3 function

I would be very grateful.

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [November 6, 2017, 4:02pm UTC](https://discourse.pymc.io/t/how-to-model-time-dependent-variables-in-pymc3/516/6 "2017-11-06T16:02:09Z")

</div>

So if I get you correctly, you have observed pairs of data (x, y), and if you plot the data (x, y) it looks like a skewed normal probability density function?

In this case, I don’t think you should reuse the function in pymc3, as in pymc3 we did not provide a density function (pdf). You can wrap the scipy function in theano using `@theano.as_op` (you can do a search in discourse, there are a few examples). However, you cannot use the gradient doing so it is discouraged.

I would suggest you to rewrite the skew normal shape function in theano, the error function is available in theano, so you can try below:

```python
import theano.tensor as tt
def skewnorm_pdf(x, e=0, w=1, a=0, mag=1):
    t = (x-e) / w
    cdf = .5*(1+tt.erf(a*t/np.sqrt(2)))
    pdf = 1/np.sqrt(2*np.pi)*tt.exp(-t**2/2)
    return 2 * mag * pdf * cdf

```

Of course, you can not model the output of a function directly (as it is deterministic). A workaround is something like

```python
with pm.Model():
    e = ...
    w = ...
    a = ...
    mag = ...
    obs = pm.Normal('y', mu=skewnorm_pdf(x, e, w, a, mag), sd=.001, observed=y)

```