# Fitting a distribution with custom functions

**URL:** <https://discourse.pymc.io/t/fitting-a-distribution-with-custom-functions/210>\
**Category:** Questions\
**Tags:** theano\
**Created:** [August 1, 2017, 4:41pm UTC](https://discourse.pymc.io/t/fitting-a-distribution-with-custom-functions/210 "2017-08-01T16:41:41Z")\
**Posts on this page:** 1\
**Showing post:** 3

<div class="post-metadata">

**Author:** ![sammosummo](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/sammosummo/32/96_2.png) [@sammosummo](https://discourse.pymc.io/u/sammosummo)\
**Post date:** [August 2, 2017, 1:24pm UTC](https://discourse.pymc.io/t/fitting-a-distribution-with-custom-functions/210/3 "2017-08-02T13:24:47Z")

</div>

Thanks! So the following seems to get past Theano:

```auto
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:

```auto
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.

---

_[View the full topic](https://discourse.pymc.io/t/fitting-a-distribution-with-custom-functions/210)._