1
$\begingroup$

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)
$\endgroup$
  • 2
    $\begingroup$ Did you look at the "disasters" example in PyMC3 example? I would start by trying to model a mixture: constant vs Poisson components. $\endgroup$ – Vladislavs Dovgalecs Oct 31 '17 at 17:03
  • $\begingroup$ @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? $\endgroup$ – Edgar H Nov 1 '17 at 8:59
  • $\begingroup$ Should I ask this question on StackOverFlow as it is more about the programming part than the statistics part? $\endgroup$ – Edgar H Nov 1 '17 at 9:31
2
$\begingroup$

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.

$\endgroup$
1
$\begingroup$

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)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.