Solved – 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
Experimental Data

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
artificial data

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      #leads per second     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) 

You will probably have more luck with PyMC3 questions on our forums:

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:

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 for more details.

Similar Posts:

Rate this post

Leave a Comment