Help with Compound Poisson Gamma distribution

Questions v5
1

Greetings,

I am referring to Haslett and Parnell (2008) to build an age-depth model for radiocarbon determinations. I wanted help with both statistics and how should I approach it to implement it using pymc.

Suppose a dataset with set of points (x, y) are monotonically increasing and we assume that increments between a segment on x-axis follow gamma distribution and number of increments (knots) between a segment follows Poisson distribution. We need to infer λ and β (α is fixed at 4). Process is defined as x(y).
we write, \Omega_j = x_s(y_{s,j}) - x_s(y_{s,j-1}) and \omega_j = y_{s,j} - y_{s,j-1}, so that

Ω_j = \sum_{i = 1}^{N(ω_j) + 1} G_i

with G_i \sim gamma(\alpha, \beta) and N(\omega_j) \sim Poisson(\lambda\omega_j). The posterior for \lambda and \beta is obtained from

\pi(\lambda, \beta\vert\Omega, \omega ) \propto \pi(\Omega\vert\alpha,\beta,\lambda,\omega)\pi(\lambda, \beta)

with \{\Omega, \omega\} denoting complete set of \{\Omega_j, \omega_j\} pairs. The likelihood \pi(\lambda, \beta\vert\Omega, \omega ) = \prod_{j=1}^{m-1}(\Omega_j\vert\alpha,\beta,\lambda,\omega_j) is now of the form that is given in equation below (probability density of CPG) with \Omega_j=x and \lambda y = \lambda\omega_j.

f_0(x;\lambda y, \alpha,\beta) = \exp(-\beta x)\exp(-\lambda y) \sum_{n=1}^{\infty} {\beta^{n\alpha}\over\Gamma\{n\alpha\}}x^{n\alpha-1}{(\lambda y)^n\over n!}
f(x;\lambda y,\alpha,\beta) = f_0(x;\lambda y,\alpha,\beta) + {\partial f_0(x;\lambda y,\alpha,\beta) \over \partial(\lambda y)}

Infinite sum is not available analytically, but I have found a fast routine Dunn and Smyth (2005) and for now I have implemented it using just numpy. And I approximate the derivative.

My question here is, what do we actually observe here? We definitely observe \omega_j = y_{s,j} - y_{s,j-1} but do we also observe \Omega? Because for now, we assume that x doesn’t have any uncertainity, but later x will have noise and we would then have to also infer x. Another question that I also have is, do I need to make CustomDist for CPG? Because like I said, when x will have noise, we will use CPG as a joint prior over it.

Thank you

2

Greetings,

I think I have given too much information in my last message. In short, I want to use CPG distribution as a prior on thetas (additional condition being that thetas must be monotonically increasing). I have wrote log likelihood of CPG distribution using numpy and I can also randomly sample from it.

I am just struggling with how can I convert likelihood function to PyTensor Op, since it has while loops. I am totally lost as to what to do next and any kind of guidance would be appreciated.

Thank you

3

Pytensor supports while loops using scan. You can see here and here for examples and documentation on scan generally. To make a while loop, you can use the pytensor.scan.utils.until helper, as shown here. The jax backend does not support while looping, but C and numba do.

If you don’t want to deal with any of that, you can also treat your numpy implementation as a black box Op. See here for an example. It focuses on a likelihood function, but the procedure for wrapping any Op is the same.