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
e = ...
w = ...
a = ...
mag = ...
obs = pm.Normal('y', mu=skewnorm_pdf(x, e, w, a, mag), sd=.001, observed=y)