# Modelling time-dependent rate using Bayesian statistics (pymc3)

How to model time-dependent variables explicitly? (or alternatively, a better approach to modelling)

I measure events over time and there are two sources: a) constant rate baseline and b) a time-dependent burst as seen below

I want to quantify the difference between these two sources of events and fit the likely parameters from the experimental data.

My modelling assumptions are:

1. Baseline comes from Poisson distribution with fixed mean $\lambda_{a}$

2. Burst comes from Poisson distribution with varying mean $\lambda_{b} (t)$

3. $\lambda_{b} (t)$ is described by a skewed normal distribution with parameters, mean $\mu$, standard deviation $\sigma$, skewness $\alpha$, magnitude $mag$

Generating artificial data as a test input. $\lambda_a$ and $\lambda_b$ look like this and leads to artificial datasets like this

What I struggle with is including the explicit time dependence. In pymc3

    import numpy as np
import scipy.stats
import pymc3 as pm

#Generate a training set
lambda_a_true = 10

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)

time = np.linspace(0.0, 300.0, 1000)
a_t =10 ; mu_t = 35; sigma_t = 25; mag_t = 35
lambda_b_true = skew(time, mu_t, sigma_t, a_t, mag_t)

ts=np.arange(301); count=np.zeros(301, dtype = np.uint16)
for ii in ts:
bl = np.random.poisson(lam=lambda_a_true)
tv = np.random.poisson(lam=skew(ii, mu_t, sigma_t, a_t, mag_t))
count[ii]=bl+tv

niter = 2000
with pm.Model() as model:
#Baseline lambda - one of our unknown
lambda_bl = pm.Uniform('lambda_bl', 0., 20)

#Parameters for skewed Gaussians - also unknowns
alpha = pm.Uniform('lambda_bl', 0., 20)
mu = pm.Uniform('lambda_bl', 0., 300)
sigma = pm.Uniform('lambda_bl', 0., 300)
mag = pm.Uniform('lambda_bl', 0., 100)

#How to include time dependence here?
lambda_tv = skew(t, mu=mu, sd=sigma, tau=None, alpha=alpha, mag=mag)

• Did you look at the "disasters" example in PyMC3 example? I would start by trying to model a mixture: constant vs Poisson components. Oct 31, 2017 at 17:03
• @VladislavsDovgalecs The problem is that it only considers two time periods. If I wanted to use that approach I would have to insert 300 switch statements. Surely there is a better way? Nov 1, 2017 at 8:59
• Should I ask this question on StackOverFlow as it is more about the programming part than the statistics part? Nov 1, 2017 at 9:31

You will probably have more luck with PyMC3 questions on our forums: http://discourse.pymc.io

It sounds like you have a mixture model here where you want to infer which samples come from which distribution. There quite a few examples which you can look at: http://docs.pymc.io/examples.html#mixture-models

While your model should look a bit different afterwards, it's also important to know that using python code in the model creation does not do what you think it does: lambda_tv = skew(t, mu=mu, sd=sigma, tau=None, alpha=alpha, mag=mag). See http://docs.pymc.io/theano.html for more details.

As @twiecki suggested, the devs at pymc3 provided me with an answer.

The key to the solution is that theano provides an error function, allowing the skewed distribution to be coded.

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, [...]

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


and hence

with pm.Model() as model:
#Baseline lambda - one of our unknown
lambda_bl = pm.Uniform('lambda_bl', 0., 20)

#Parameters for skewed Gaussians - also unknowns
mu = pm.Uniform('mu', 0., 300, shape=N_SPOTS)
sigma = pm.Uniform('sigma', 0., 50, shape=N_SPOTS)
mag = pm.Uniform('mag', 0., 100, shape=N_SPOTS)
alpha = pm.Uniform('alpha', 0., 50, shape=N_SPOTS)

lam = lambda_bl + tt.sum( skew_theano(ts,e=mu,w=sigma,a=alpha, mag=mag) )
cs = pm.Poisson(mu=lam, name='cs', observed=count)