How to model time-dependent variables in pymc3

Questions
1

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.

2

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/stochastic_volatility.html
and the random walk related examples:
http://docs.pymc.io/notebooks/survival_analysis.html?highlight=random%20walk

1 Like
3

@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 nopm.math.cumsum().… My question originated after reading this vignette:

https://cran.r-project.org/web/packages/walker/vignettes/walker.html

4

You can try cumsum from theano:
from theano.tensor.extra_ops import cumsum

5

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.

6

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:

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

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