Solved – Modelling time-dependent rate using Bayesian statistics (pymc3)

bayesianpymctime series

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
lambdas
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)

Best Answer

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.

Related Question