Hello

I’m having trouble sampling from a model. The likelihood of the data comes from a Poisson single hit process which simplifies to: \text{Bernoulli}(p=f(x, v)) with f(x, v) = 1-2^{-10^{x+v}}.

Here’s some test data, and what this function looks like with v=4.

```
import numpy as np
import matplotlib.pyplot as plt
import pymc as pm
def f(v, x):
return 1.0 - 2.0 ** (-10.0 ** (x + v))
status = np.array([1, 1, 1, 1, 1, 0, 0, 0])
x = np.array([0, -1, -2, -3, -4, -5, -6, -7])
x_plot = np.linspace(-8, 0)
plt.plot(x_plot, f(v=4, x=x_plot), label="f(v=4, x=x)")
plt.scatter(x, status, label="Data")
plt.ylabel("P(positive)")
plt.xlabel("x")
plt.legend()
```

Running this model, even with a prior bang on the expected value, gives 100% divergences:

```
with pm.Model():
v = pm.Normal("v", 4.0, 1.0)
pm.Bernoulli("lik", p=f(v=v, x=x), observed=status)
pm.sample()
```

Presumably the issue comes from the fact that `v`

is an exponent of an exponent. I can’t think of how to reparameterize this function. Does anyone have any ideas? Or is there anything else I can try?

Thanks in advance

David