# Help with Compound Poisson Gamma distribution

**URL:** <https://discourse.pymc.io/t/help-with-compound-poisson-gamma-distribution/16621>\
**Category:** v5\
**Tags:** modeling\
**Created:** [March 5, 2025, 4:19am UTC](https://discourse.pymc.io/t/help-with-compound-poisson-gamma-distribution/16621 "2025-03-05T04:19:59Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![Shivraj1906](https://avatars.discourse-cdn.com/v4/letter/s/e95f7d/32.png) [@Shivraj1906](https://discourse.pymc.io/u/Shivraj1906)\
**Post date:** [March 5, 2025, 4:19am UTC](https://discourse.pymc.io/t/help-with-compound-poisson-gamma-distribution/16621/1 "2025-03-05T04:19:59Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Shivraj1906](https://avatars.discourse-cdn.com/v4/letter/s/e95f7d/32.png) [@Shivraj1906](https://discourse.pymc.io/u/Shivraj1906)\
**Post date:** [March 19, 2025, 5:47am UTC](https://discourse.pymc.io/t/help-with-compound-poisson-gamma-distribution/16621/2 "2025-03-19T05:47:27Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![jessegrabowski](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jessegrabowski/32/5010_2.png) [@jessegrabowski](https://discourse.pymc.io/u/jessegrabowski)\
**Post date:** [March 19, 2025, 6:26am UTC](https://discourse.pymc.io/t/help-with-compound-poisson-gamma-distribution/16621/3 "2025-03-19T06:26:56Z")

</div>

Pytensor supports while loops using `scan`. You can see [here](https://pytensor.readthedocs.io/en/latest/library/scan.html) and [here](https://pytensor.readthedocs.io/en/latest/gallery/scan/scan_tutorial.html) for examples and documentation on scan generally. To make a while loop, you can use the `pytensor.scan.utils.until` helper, as shown [here](https://pytensor.readthedocs.io/en/latest/library/scan.html#conditional-ending-of-scan). 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](https://www.pymc.io/projects/examples/en/latest/howto/blackbox_external_likelihood_numpy.html) for an example. It focuses on a likelihood function, but the procedure for wrapping any Op is the same.