Fitting a distribution with custom functions

Questions
3

Thanks! So the following seems to get past Theano:

import pymc3 as pm
import numpy as np
import theano.tensor as tt
import theano
from scipy.special import hyp2f1, gammaln
from scipy.stats import multivariate_normal, pearsonr
import matplotlib.pyplot as plt


@theano.as_op(itypes=[tt.dscalar, tt.dvector], otypes=[tt.dvector])
def _hyp2f1(c, z):
    return hyp2f1(0.5, 0.5, c, z)


class SampleCorrelationDistribution(pm.Continuous):
    """Distribution of sample correlations.

    """

    def __init__(self, rho, n, *args, **kwargs):

        self._rho = tt.as_tensor_variable(rho)
        self._n = tt.as_tensor_variable(n)

        super(SampleCorrelationDistribution, self).__init__(*args, **kwargs)

    def logp(self, value):

        if value.ndim != 1:
            raise ValueError('Only one-dimensional values are supported.')
        rho = self._rho
        n = self._n

        lognum = tt.log(n - 2) + tt.gammaln(n - 1) + ((n - 1) / 2) * tt.log(
            1 - rho ** 2) + ((n - 4) / 2) * tt.log(1 - value ** 2)

        logden = tt.log(tt.sqrt(2 * np.pi)) + tt.gammaln(n - 0.5) + (
            n - 1.5) * tt.log(1 - rho * value)

        hyp = _hyp2f1((2 * n - 1) / 2, (rho * value) / 2)

        loghyp = tt.log(hyp)

        return tt.sum(lognum - logden * loghyp)


with pm.Model():

    true_rho = 0
    true_n = 100
    data = []

    for it in range(1000):

        xy = multivariate_normal.rvs([0, 0], [[1, true_rho], [true_rho, 1]],
                                     size=true_n)
        x, y = xy[:, 0], xy[:, 1]
        r, p = pearsonr(x, y)
        data.append(r)

    rho = pm.Uniform(name='rho')
    y = SampleCorrelationDistribution(name='y', rho=rho, n=true_n, observed=data)
    trace = pm.sample()
    pm.traceplot(trace)
    plt.savefig('traceplot.png')

But it fails in PyMC3:

AttributeError: 'FromFunctionOp' object has no attribute 'grad'

Perhaps there should be some catch in PyMC3 to use non-gradient methods instead of raising an error?

For this particular problem I know what the gradient should be, it it’s not super clear to me how to define that in Theano.