Scipy.stats distributions cdfs

Questions
1

Does pymc3 support scipy.stats distributions cdfs and ppfs?
From what I could find, yes. But I am not being able to implement myself.
Specially, I want to perform a transformation on a MvNormal variable, but pymc3 seems not to accept the use of scipy distributions.

Am I doing something wrong?
The code below is just an example and is incomplete because I will do other things with the transformed variable

import numpy as np
import pymc3 as pm
import theano.tensor as tt
from theano.compile.ops import as_op
from scipy.stats import norm, beta

S = 10
alphas = np.random.poisson(4,S)
betas = np.random.poisson(10,S)
hyper_mu = np.zeros(S)

@as_op(itypes=[tt.dmatrix,tt.dvector,tt.dvector],otypes=[tt.dmatrix])
def fn (P,alphas,betas):
    N = P.shape[0]
    out = np.zeros(shape=(N,N))
    for i in range(0, N):
        out[i] = beta.ppf(norm.cdf(P[i]), alphas[i], betas[i])
    return out

with pm.Model() as model:
    packed_L = pm.LKJCholeskyCov('packed_L',
                                 n=S, eta=1,
                                 sd_dist=pm.HalfCauchy.dist(2.5))
    chol = pm.expand_packed_triangular(S, packed_L)
    M = pm.MvNormal('M', mu=hyper_mu, chol=chol, shape=(S, S))
    Q = fn(M, alphas, betas)

Error message:

Traceback (most recent call last):
  File "/usr/lib/python3/dist-packages/IPython/core/interactiveshell.py", line 2882, in run_code
    exec(code_obj, self.user_global_ns, self.user_ns)
  File "<ipython-input-2-a61aec453dcb>", line 26, in <module>
    Q = fn(M, alphas, betas)
  File "/usr/local/lib/python3.6/dist-packages/theano/gof/op.py", line 615, in __call__
    node = self.make_node(*inputs, **kwargs)
  File "/usr/local/lib/python3.6/dist-packages/theano/gof/op.py", line 983, in make_node
    if not all(inp.type == it for inp, it in zip(inputs, self.itypes)):
  File "/usr/local/lib/python3.6/dist-packages/theano/gof/op.py", line 983, in <genexpr>
    if not all(inp.type == it for inp, it in zip(inputs, self.itypes)):
AttributeError: 'numpy.ndarray' object has no attribute 'type'
2

Use as_op is likely very inefficient, there is an effort to introduce cdf to the distribution class but progress is a bit slow: https://github.com/pymc-devs/pymc3/pull/2688

Having said that, it is possible to make it work with as_op as well, you can try:

@as_op(itypes=[tt.dmatrix, tt.dvector, tt.dvector], otypes=[tt.dmatrix])
def fn(P, alphas, betas):
    return beta.ppf(norm.cdf(P), alphas, betas)

import theano

with pm.Model() as model:
    packed_L = pm.LKJCholeskyCov('packed_L',
                                 n=S, eta=1,
                                 sd_dist=pm.HalfCauchy.dist(2.5))
    chol = pm.expand_packed_triangular(S, packed_L)
    M = pm.MvNormal('M', mu=hyper_mu, chol=chol, shape=(S, S))
    Q = fn(M, theano.shared(alphas.astype('float64')), theano.shared(betas.astype('float64')))
1 Like