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

**URL:** <https://discourse.mc-stan.org/t/confluent-hypergeometric-functions-to-estimate-population-size-and-detection-probability-from-spatially-replicated-counts/28758>\
**Category:** Modeling\
**Tags:** ecology, techniques\
**Created:** [September 6, 2022, 11:21pm UTC](https://discourse.mc-stan.org/t/confluent-hypergeometric-functions-to-estimate-population-size-and-detection-probability-from-spatially-replicated-counts/28758 "2022-09-06T23:21:53Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![nikunj410](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/nikunj410/32/13085_2.png) [@nikunj410](https://discourse.mc-stan.org/u/nikunj410)\
**Post date:** [September 6, 2022, 11:21pm UTC](https://discourse.mc-stan.org/t/confluent-hypergeometric-functions-to-estimate-population-size-and-detection-probability-from-spatially-replicated-counts/28758/1 "2022-09-06T23:21:53Z")

</div>

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)](https://onlinelibrary.wiley.com/doi/10.1111/j.0006-341X.2004.00142.x) 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 dist.} \underbrace{p(\theta | \alpha, \beta)}\_\text{Beta dist.} \underbrace{p( n\_{x,t,j}| \lambda\_{x,t} )}\_\text{Poisson dist.} d\theta.

This gives us  
p(Y= y\_{x,t,j} | \alpha, \beta, \lambda\_{x,t}) = \text{Poisson}(y\_{x,t,j} |\lambda\_{x,t}) \frac{\text{Beta}(\alpha + y\_{x,t,j}, \beta)}{\text{Beta}(\alpha, \beta)} \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](https://people.maths.ox.ac.uk/porterm/papers/hypergeometric-final.pdf) 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.

---

<div class="post-metadata">

**Author:** ![spinkney](https://avatars.discourse-cdn.com/v4/letter/s/dec6dc/32.png) [@spinkney](https://discourse.mc-stan.org/u/spinkney)\
**Post date:** [September 7, 2022, 12:30am UTC](https://discourse.mc-stan.org/t/confluent-hypergeometric-functions-to-estimate-population-size-and-detection-probability-from-spatially-replicated-counts/28758/2 "2022-09-07T00:30:43Z")

</div>

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](https://github.com/stan-dev/math/blob/develop/stan/math/prim/fun/hypergeometric_pFq.hpp). @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](https://github.com/stan-dev/math/blob/develop/stan/math/prim/fun/grad_pFq.hpp).

---

<div class="post-metadata">

**Author:** ![maxbiostat](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/maxbiostat/32/1326_2.png) [@maxbiostat](https://discourse.mc-stan.org/u/maxbiostat)\
**Post date:** [September 7, 2022, 12:33am UTC](https://discourse.mc-stan.org/t/confluent-hypergeometric-functions-to-estimate-population-size-and-detection-probability-from-spatially-replicated-counts/28758/3 "2022-09-07T00:33:05Z")

</div>

@GuidoAMoreira and I have some [Stan code](https://github.com/GuidoAMoreira/stan_summer) to do adaptive truncation, described in detail in our [preprint](https://arxiv.org/abs/2202.06121). @jsocolar, this might be the push we need to get this implemented for occupancy models, uh? Happy to collaborate!

---

<div class="post-metadata">

**Author:** ![andrjohns](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/andrjohns/32/15297_2.png) [@andrjohns](https://discourse.mc-stan.org/u/andrjohns)\
**Post date:** [September 7, 2022, 8:21am UTC](https://discourse.mc-stan.org/t/confluent-hypergeometric-functions-to-estimate-population-size-and-detection-probability-from-spatially-replicated-counts/28758/4 "2022-09-07T08:21:50Z")

</div>

I’m currently working on the confluent 1F1 and its gradients over on [this math branch](https://github.com/stan-dev/math/compare/develop...andrjohns:math:issue-2764-hyper1f1?expand=1) (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

---

<div class="post-metadata">

**Author:** ![Matthew\_Spencer1](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/matthew_spencer1/32/4616_2.png) [@Matthew\_Spencer1](https://discourse.mc-stan.org/u/Matthew_Spencer1)\
**Post date:** [September 8, 2022, 10:21pm UTC](https://discourse.mc-stan.org/t/confluent-hypergeometric-functions-to-estimate-population-size-and-detection-probability-from-spatially-replicated-counts/28758/5 "2022-09-08T22:21:13Z")

</div>

I had a similar problem a couple of months ago:  
[https://discourse.mc-stan.org/t/marginalizing-a-double-binomial-poisson-hierarchical-distribution/27801](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)](https://doi.org/10.2307/1402230), which is less general but a bit easier to follow.

---

<div class="post-metadata">

**Author:** ![nikunj410](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/nikunj410/32/13085_2.png) [@nikunj410](https://discourse.mc-stan.org/u/nikunj410)\
**Post date:** [February 7, 2023, 7:16am UTC](https://discourse.mc-stan.org/t/confluent-hypergeometric-functions-to-estimate-population-size-and-detection-probability-from-spatially-replicated-counts/28758/6 "2023-02-07T07:16:33Z")

</div>

Hi @andrjohns, was 1F1 release in the last version of stan? I could not find it in the function manual.

---

<div class="post-metadata">

**Author:** ![nikunj410](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/nikunj410/32/13085_2.png) [@nikunj410](https://discourse.mc-stan.org/u/nikunj410)\
**Post date:** [October 20, 2023, 6:20pm UTC](https://discourse.mc-stan.org/t/confluent-hypergeometric-functions-to-estimate-population-size-and-detection-probability-from-spatially-replicated-counts/28758/7 "2023-10-20T18:20:56Z")

</div>

I wrote a Stan function to calculate 1F1 on log scale based on this [paper](https://people.maths.ox.ac.uk/porterm/papers/hypergeometric-final.pdf). Posting it here in case it may be useful to someone. This may become obsolete once [pFq](https://github.com/stan-dev/math/pull/2961) are implemented in Stan.

```plaintext
  real Log_1F1(real a, real b, real z) {
    int j = 0;
    real toll = 1e-9;
    real log_Aj = 0;
    real log_Sj = log_Aj;
    
    // z > 0 ? Log_1F1(a, b, z) : z + Log_1F1(b-a, b, -z) (for numericall accuracy)
    
    if (z >= 0){
      while (log_Aj-log_Sj > log(toll)){
        j += 1;
        log_Aj += log(((a+j-1)/(b+j-1))*(z/j));
        log_Sj = log_sum_exp([log_Aj, log_Sj]);}
      return log_Sj;
    } else {
        while (log_Aj-log_Sj > log(toll)){
          j += 1;
          log_Aj += log(((b-a+j-1)/(b+j-1))*(-z/j));
          log_Sj = log_sum_exp([log_Aj, log_Sj]);}
        return z + log_Sj;
    }
  }

```

---

<div class="post-metadata">

**Author:** ![luchinoprince](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/luchinoprince/32/16355_2.png) [@luchinoprince](https://discourse.mc-stan.org/u/luchinoprince)\
**Post date:** [November 29, 2023, 4:05pm UTC](https://discourse.mc-stan.org/t/confluent-hypergeometric-functions-to-estimate-population-size-and-detection-probability-from-spatially-replicated-counts/28758/8 "2023-11-29T16:05:04Z")

</div>

So to this date there is no native function right? I could not find it ether and would like to test a model which needs it.

---

<div class="post-metadata">

**Author:** ![nikunj410](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/nikunj410/32/13085_2.png) [@nikunj410](https://discourse.mc-stan.org/u/nikunj410)\
**Post date:** [November 29, 2023, 7:29pm UTC](https://discourse.mc-stan.org/t/confluent-hypergeometric-functions-to-estimate-population-size-and-detection-probability-from-spatially-replicated-counts/28758/9 "2023-11-29T19:29:49Z")

</div>

@luchinoprince I dont think there is one but it seems to be under works. I also implemented pFq function.  
Here is the code that I wrote. I have only tested this for a\_p, b\_q and Re(z) \> 0

```no-highlight
  real log_pFq(vector a_p, vector b_q, real z) {
    int j = 0;
    real toll = log(machine_precision());
    real log_Aj = 0;
    real log_Sj = log_Aj;
    real log_z = log(z);
    
    while (log_Aj-log_Sj > toll && j< 1e6){
      j += 1;
      log_Aj += sum(log(a_p+j-1))-sum(log(b_q+j-1))+log_z-log(j);
      log_Sj = log_sum_exp([log_Aj, log_Sj]);}
    return log_Sj;
  }

```

---

<div class="post-metadata">

**Author:** ![luchinoprince](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/luchinoprince/32/16355_2.png) [@luchinoprince](https://discourse.mc-stan.org/u/luchinoprince)\
**Post date:** [November 29, 2023, 7:42pm UTC](https://discourse.mc-stan.org/t/confluent-hypergeometric-functions-to-estimate-population-size-and-detection-probability-from-spatially-replicated-counts/28758/10 "2023-11-29T19:42:32Z")

</div>

@nikunj410 thanks for the answer, I saw you implementation. I am just concerned about the speed convergence of such an implementation and its potential instability. Also not sure how the gradient would work if the pdf is an iterative function. In my case I have a tensor where the entries would be a density which uses pFq, hence I would have to calculate a lot of those at every HMC step, which I guess would be too costly. Maybe parallelising within chain could help, I will try to test it,

---

<div class="post-metadata">

**Author:** ![nikunj410](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/nikunj410/32/13085_2.png) [@nikunj410](https://discourse.mc-stan.org/u/nikunj410)\
**Post date:** [November 29, 2023, 8:13pm UTC](https://discourse.mc-stan.org/t/confluent-hypergeometric-functions-to-estimate-population-size-and-detection-probability-from-spatially-replicated-counts/28758/11 "2023-11-29T20:13:06Z")

</div>

@luchinoprince I am pretty sure that this is not a stable implementation when a,b,z \<0. For, example, when I tried implementing 1F1 for z\<0, I was getting convergence issues. For positive coefficient values, I was able to achieve similar numerical accuracy as the pFq function in Julia (because I was working on log scale).

You are right the speed is a major limiting factor. Since there are no analytical gradients, we have to rely on auto diff which is usually slower. But the work on gradients is stills in progress [here](https://github.com/stan-dev/math/issues/2509).

Maybe you want to try and implement the code I shared to see if your Stan code works with synthetic data and when the native implementation of pFq is made available in Stan, you can make appropriate changes.

---

<div class="post-metadata">

**Author:** ![luchinoprince](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/luchinoprince/32/16355_2.png) [@luchinoprince](https://discourse.mc-stan.org/u/luchinoprince)\
**Post date:** [November 30, 2023, 2:35pm UTC](https://discourse.mc-stan.org/t/confluent-hypergeometric-functions-to-estimate-population-size-and-detection-probability-from-spatially-replicated-counts/28758/12 "2023-11-30T14:35:14Z")

</div>

hey @nikunj410 , I am sure that in my case a,b would be greater than 0, while z could have ether sign(not sure). I will see if I can find other distributions which allow a better treatment, as the main issue with approximate likelihood from my perspective is that it is hard to justify/specify from what distribution we are then sampling from.

Best,  
Luca

---

<div class="post-metadata">

**Author:** ![nikunj410](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/nikunj410/32/13085_2.png) [@nikunj410](https://discourse.mc-stan.org/u/nikunj410)\
**Post date:** [March 28, 2026, 2:18am UTC](https://discourse.mc-stan.org/t/confluent-hypergeometric-functions-to-estimate-population-size-and-detection-probability-from-spatially-replicated-counts/28758/13 "2026-03-28T02:18:03Z")

</div>

If anyone is still interested, I used Claude Code to write a custom C++ code implementing a vectorized version of the Log\_1F1 function with analytical gradients. My model runs 4X faster.

In R, you can call it using

```no-highlight
file <- file.path("path_to_stan_file/stan_code.stan")
header_file <- file.path("path_to_hpp_file/hypergeometric_1f1.hpp")
[Uploading: hypergeometric_1f1.hpp…]()

mod <- cmdstan_model(
  stan_file = file,
  cpp_options = list(
    user_header = header_file,
    stan_threads = TRUE # enables within-chain parallelism
  ),
  stanc_options = list("O1","allow-undefined" = TRUE),
  force_recompile = TRUE
)

```

In stan, you can call the function using

```no-highlight
vector Log_1F1(vector a, vector b, vector z);

```

Sometimes the code causes some chains to finish unexpectedly. I change the seed, and run it again.

PS: I have not rigorously tested the code.

[hypergeometric\_1f1.hpp](https://discourse.mc-stan.org/uploads/short-url/sU6GvjyxO4i5u19HzU3lJxubaNt.hpp) (16.1 KB)