Custom distribution

[davidbrochart]

It looks like the creation of a custom Continuous distribution is broken in PyMC3. In the example below, I get the following error:

ValueError: expected an ndarray
Apply node that caused the error: Elemwise{Composite{((i0 + (i1 * i2)) - (i3 + (i4 * i5)))}}[(0, 0)](FromFunctionOp{beta_logp}.0, TensorConstant{0.5}, Sum{acc_dtype=float64}.0, FromFunctionOp{beta_logp}.0, TensorConstant{0.5}, Sum{acc_dtype=float64}.0)
Toposort index: 10
Inputs types: [TensorType(float64, scalar), TensorType(float64, scalar), TensorType(float64, scalar), TensorType(float64, scalar), TensorType(float64, scalar), TensorType(float64, scalar)]
Inputs shapes: [(), (), (), (), (), ()]
Inputs strides: [(), (), (), (), (), ()]
Inputs values: [-0.2159736732943715, array(0.5), array(-458.9112207578818), -0.0, array(0.5), array(-400.5085556528546)]
Outputs clients: [['output']]
HINT: Re-running with most Theano optimization disabled could give you a back-trace of when this node was created. This can be done with by setting the Theano flag 'optimizer=fast_compile'. If that does not work, Theano optimizations can be disabled with 'optimizer=None'.
HINT: Use the Theano flag 'exception_verbosity=high' for a debugprint and storage map footprint of this apply node.

from pymc3.distributions import Continuous, Normal
from pymc3 import Model, sample, Metropolis
from theano import as_op
import theano.tensor as T
import numpy as np

size = 100
b = 2

X = np.linspace(0, 1, size)
Y = b * X + np.random.randn(size)

class Beta(Continuous):
    def __init__(self, mu, *args, **kwargs):
        super(Beta, self).__init__(*args, **kwargs)
        self.mu = mu
        self.mode = mu
    def logp(self, value):
        mu = self.mu
        return beta_logp(value - mu)

@as_op(itypes=[T.dscalar], otypes=[T.dscalar])
def beta_logp(value):
    return -1.5 * np.log(1 + (value) ** 2)

with Model() as model:
    # Create custom densities
    beta = Beta('slope', mu=0, testval=0)
    # Create likelihood
    like = Normal('y_est', mu=beta * X, sd=1, observed=Y)
    step = Metropolis([beta])
    trace = sample(100, step)

[bwengals]

To make your example work, there is no need to define a custom op. Do you need to make a custom theano op for your application? I commented out #@as_op and replacing np.log with T.log.

from pymc3.distributions import Continuous, Normal
from pymc3 import Model, sample, Metropolis
from theano import as_op
import theano.tensor as T
import numpy as np

size = 100 
b = 2 

X = np.linspace(0, 1, size)
Y = b * X + np.random.randn(size)

class Beta(Continuous):
    def __init__(self, mu, *args, **kwargs):
        super(Beta, self).__init__(*args, **kwargs)
        self.mu = mu
        self.mode = mu
    def logp(self, value):
        mu = self.mu
        return beta_logp(value - mu) 

#@as_op(itypes=[T.dscalar], otypes=[T.dscalar])
def beta_logp(value):
    return -1.5 * T.log(1 + (value) ** 2)

with Model() as model:
    # Create custom densities
    beta = Beta('slope', mu=0, testval=0)
    # Create likelihood
    like = Normal('y_est', mu=beta * X, sd=1, observed=Y)
    step = Metropolis([beta])
    trace = sample(100, step)

Also, consider using DensityDist. In your model, replace beta = ... with

    beta = pymc3.DensityDist('slope', lambda value: -1.5 * T.log(1 + value**2), testval=0)

Custom theano ops are required when theano does not natively support the operation you need to compute.

[davidbrochart]

Yes, that was just an example but I need to make a custom Continuous distribution. Basically I want to use the distribution of my posteriors as priors to perform an update, which I think is not possible in Theano (it involves computing a KDE or a histogram from the trace of the posterior samples).

[davidbrochart]

Returning explicitly an np.array in the beta_logp function seems to solve the issue:

@as_op(itypes=[T.dscalar], otypes=[T.dscalar])
def beta_logp(value):
    return np.array(-1.5 * np.log(1 + (value) ** 2))

[twiecki]

@davidbrochart Want to write up a short tutorial as a NB and do a PR?

[davidbrochart]

I'd love to! I'm actually writing a NB on using PyMC3 in hydrology, I'll let you know when it's ready.