Confluent hypergeometric functions to estimate population size and detection probability from spatially replicated counts

I am modeling spatial population dynamics which consist of two layers. The inner layer is the population abundance of the species which is modeled using an a spatial population dynamics model. The outer layer is the observation process.

Here is the brief description of the model

Assuming count data y_{x,t,j}\sim \text{Binomial}(n_{x,t,j}, \theta), were x is the patch location, t is time, j is the replicate sampling ID, and \theta is the detection probability. We assume \theta \sim \text{Beta}(\alpha , \beta) and n_{x,t,j} \sim \text{Poisson} (\lambda _{x,t}), where \lambda _{x,t} comes from the population dynamics model e.g. \lambda _{x,t+1} = \int^{\infty}_{- \infty} k(x-y) f(\lambda _{y,t}) dy.

Royle (2004) suggests that if we have replicate counts, we can estimate p, n_{x,t,j}, and \lambda _{x,t}. However, this problem cannot be coded directly in Stan because n_{x,t,j} is a discrete parameter and, therefore, has to be marginalized out. Mathematically, we can marginalize \theta and n_{x,t,j} in the following manner

p(Y= y_{x,t,j} | \alpha, \beta, \lambda_{x,t}) = \sum_{n_{x,t,j} =0}^{\infty} \int_{0}^{1} \underbrace{p(y_{x,t,j} | \theta, n_{x,t,j} )}_\text{Binomial} \underbrace{p(\theta | \alpha, \beta)}_\text{Beta} \underbrace{p( n_{x,t,j}| \lambda_{x,t} )}_\text{Poisson} d\theta.

This stack exchange article solves the above integral to show
p(Y= y_{x,t,j} | \alpha, \beta, \lambda_{x,t}) = \text{Poisson}(y_{x,t,j} |\lambda_{x,t}) \frac{(\alpha)_{y_{x,t,j}}}{(\alpha +\beta)_{y_{x,t,j}}} \ce{_{1}F_{1}} (\beta, \alpha + \beta + y_{x,t,j}, \lambda_{x,t}),

where \ce{_{1}F_{1}} is the confluent hypergeometric function of the first kind. We have successfully marginalized the \theta and n_{x,t,j}, which can be recovered in the generated quantities block. Additionally, there might be computational advantages in terms of both accuracy and speed (not tested yet).

I am stuck trying to evaluate \ce{_{1}F_{1}}. Currently Stan does not support confluent hypergeometric functions. There are several infinite summation and integral methods to estimate \ce{_{1}F_{1}} but they may not be stable and scalable. Is the Stan development team working on implementing a scalable and reliable method/function to calculate confluent hypergeometric functions anytime soon?

I believe that the technique I presented in this post could be very useful in patch occupancy modeling.

1 Like

Hypergeometric functions are in Stan math and should be exposed in the Stan language with the next release. See math/hypergeometric_pFq.hpp at develop · stan-dev/math · GitHub. @andrjohns implemented this and will know more.

The interesting parts are in the derivative s that he implemented at math/grad_pFq.hpp at develop · stan-dev/math · GitHub.

2 Likes

@GuidoAMoreira and I have some Stan code to do adaptive truncation, described in detail in our preprint. @jsocolar, this might be the push we need to get this implemented for occupancy models, uh? Happy to collaborate!

1 Like

I’m currently working on the confluent 1F1 and its gradients over on this math branch (since it’s the foundation of a bunch of numerical stability improvements for the log-CDFs).

So this should be available for you in the next release

3 Likes

I had a similar problem a couple of months ago:
https://discourse.mc-stan.org/t/marginalizing-a-double-binomial-poisson-hierarchical-distribution/27801
I had just two replicate counts, although it sounds as though you may have more. I didn’t marginalize out the detection probabilities, which may affect whether my approach is relevant to you. I initially looked at ways of evaluating the confluent hypergeometric function, including the adaptive truncation suggested by @maxbiostat, but in the end found that there was a relatively easy exact solution if you work with the probability generating functions. There is Stan code at the end of the post in the link above (I may have made a few minor changes since then, but that version was working). I recommend first having a look at section 36.8 in Johnson, Kotz and Balakrishnan, Discrete Multivariate Distributions, which is quite general, and then Hamdan and Tsokos (1971), which is less general but a bit easier to follow.

1 Like