Scipy functions in pymc3

1

Hello all,

I am trying to reproduce some literature results, this time utilizing the special function “Confluent hypergeometric function 1F1”, OK I didn’t know it either.

The paper is: Coath, Christopher D., Robert CJ Steele, and W. Fred Lunnon. “Statistical bias in isotope ratios.” Journal of Analytical Atomic Spectrometry 28.1 (2013): 52-58.

The issue is a non-biased estimator for expected value of the ratio of two Poisson variables.

r'(X, Y)=\frac{X}{Y+1}\cdot M(1, Y+2, 0) where M is that special function above and assuming no noise.

This seems like it should be straightforward enough but I don’t understand something about Theano maybe. My attempt at a simple model is:

import pymc3 as pm
import numpy as np
import scipy
import scipy.special as sc
import theano.tensor as tt
from theano.compile.ops import as_op

@as_op(itypes=[tt.lscalar], otypes=[tt.dscalar])
def hyp1f1(value):
    return sc.hyp1f1(1, value+2, 0)


with pm.Model() as model:
    p1 = pm.Poisson('p1', 10)
    p2 = pm.Poisson('p2', 15)
    ratio = pm.Deterministic('ratio', p1/p2)
    quasi = pm.Deterministic('quasi', ((p1)/p2+1)*hyp1f1(p2))
    step = pm.Metropolis()
    trace = pm.sample(10000, step=step)

But try as I may I cannot escape from
ValueError: expected an ndarray or TypeError: For compute_test_value, one input test value does not have the requested type.

Anyone have any help/fixes for this? Might be simple.

Thanks!