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:
Baseline comes from Poisson distribution with fixed mean $lambda_{a}$
Burst comes from Poisson distribution with varying mean $lambda_{b} (t)$
$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 #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.