How to model time-dependent variables in pymc3

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)
2 Likes