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:


*

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

 A: 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.
A: 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)

