Replace bessel interpolant to scipy function

I think that when we set up the phase-marginalized likelihood the required scipy function i0e was not vectorized and so we used an interpolant to efficiently evaluate.

The new function is faster and doesn't require an interpolant.

I tried to make it so that everything won't break for users who have subclassed this likelihood, but I think that needs a little touching up.

In [1]: interpolant = interp1d(np.logspace(-5, 10, int(1e6)), np.log(i0e(np.logspace(-5, 10, int(1e6)))) + np.logspace(-5, 10, int(1e6)))

In [2]: %timeit np.log(interpolant(np.exp(np.random.uniform(-5, 10))))
20.2 µs ± 115 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

In [3]: %timeit ln_i0(np.exp(np.random.uniform(-5, 10)))
3.83 µs ± 18.9 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)