I seem unable to find an example for this anywhere: how can a time-dependent variable be modelled within the framework of pymc3?

Currently, I want to model an inhomogeneous Poisson process. For the sake of simplicity, lets say
\lambda (t) = a *t + b
where a and b are parameters I wish to determine from some observational data.

My best attempt at implementing this so far is below. I am not sure how to include the time-dependence.

import numpy as np
import pymc3 as pm
import theano
import theano.tensor as tt
import matplotlib.pylab as plt
# Observation
a_actual = 1.3
b_actual = 2.0
t = np.arange(10)
obs = np.random.poisson(a_actual * t + b_actual)
plt.figure()
plt.plot(t, obs)
plt.show()
a1 = tt.scalar('a')
b1 = tt.vector('b')
y = tt.vector('y')
out = a1 * y + b1
func = theano.function([a1, b1, y], [out])
niter = 3000
# Model
with pm.Model() as model:
a = pm.Uniform(name='a', lower=0, upper=10)
b = pm.Uniform(name='b', lower=0, upper=10)
#Expected value
lam = pm.Deterministic('lam', func(a, b, t))
count = pm.Poisson(mu=lam, name='count', observed=obs)
#Run Monte Carlo
trace = pm.sample(niter)
pm.plots.plot_posterior(trace=trace["a"])
pm.plots.autocorrplot(trace=trace, varnames=["a"]);
pm.plots.plot_posterior(trace=trace["b"])
pm.plots.autocorrplot(trace=trace, varnames=["b"]);
plt.show()

In your example above, the lam is a linear function of t, so you dont need to compile a theano function. Doing lam = pm.Deterministic('lam', a*t + b) is sufficient.

I was wondering if you know of an easy way to do incorporate a more complex function, specifically a skewed normal distribution.

I am trying to fit the function that describes a skewed normal distribution to some data because this function appears to fit well. While there is a SkewNormal in pymc3, it is used as a probability function and I am not sure how to use pymc3 distribution as a function.

Normally, in Python I would define the function as

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)

However, this includes a definite integral when not using scipy function.

If you know a good way to:
a) implement this function
b) reuse the scipy functions within pymc3
c) reuse the pymc3 function

So if I get you correctly, you have observed pairs of data (x, y), and if you plot the data (x, y) it looks like a skewed normal probability density function?

In this case, I don’t think you should reuse the function in pymc3, as in pymc3 we did not provide a density function (pdf). 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, so you can try below:

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

Of course, you can not model the output of a function directly (as it is deterministic). A workaround is something like

with pm.Model():
e = ...
w = ...
a = ...
mag = ...
obs = pm.Normal('y', mu=skewnorm_pdf(x, e, w, a, mag), sd=.001, observed=y)