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.