# Continuous versions of Poisson, Binomial, NegativeBinomial

**URL:** <https://discourse.julialang.org/t/continuous-versions-of-poisson-binomial-negativebinomial/115781>\
**Category:** Statistics\
**Tags:** distributions\
**Created:** [June 17, 2024, 11:32pm UTC](https://discourse.julialang.org/t/continuous-versions-of-poisson-binomial-negativebinomial/115781 "2024-06-17T23:32:56Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![sdwfrost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sdwfrost/32/2831_2.png) [@sdwfrost](https://discourse.julialang.org/u/sdwfrost)\
**Post date:** [June 17, 2024, 11:32pm UTC](https://discourse.julialang.org/t/continuous-versions-of-poisson-binomial-negativebinomial/115781/1 "2024-06-17T23:32:56Z")

</div>

Dear All,

For one of my projects, I need to have continuous versions of the Poisson, Binomial and NegativeBinomial distributions, along the lines of [Padellini and Rue](https://repository.kaust.edu.sa/server/api/core/bitstreams/4bd7c7bb-aee4-4bad-ad48-08562484880c/content). Before I start the tedious task of porting these over to a Distributions.jl-compatible API, has anyone else done this to save me the trouble?

---

<div class="post-metadata">

**Author:** ![sethaxen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sethaxen/32/35604_2.png) [@sethaxen](https://discourse.julialang.org/u/sethaxen)\
**Post date:** [June 26, 2024, 8:41pm UTC](https://discourse.julialang.org/t/continuous-versions-of-poisson-binomial-negativebinomial/115781/2 "2024-06-26T20:41:46Z")

</div>

As far as I understand it, these continuous analogs are formed by finding a continuous distribution whose CDF evaluates to the same as the CDF of the discrete version at all points in the discrete support. In that case, here are some basic implementations to get you started:

```julia
using StatsFuns, Distributions, ADTypes, FiniteDiff
import DifferentiationInterface as DI

struct ContinuousPoisson{T<:Real,B<:ADTypes.AbstractADType} <: ContinuousUnivariateDistribution
    λ::T
    ad_type::B
end

ContinuousPoisson(λ) = ContinuousPoisson(λ, ADTypes.AutoFiniteDiff())

_discretize(d::ContinuousPoisson) = Poisson(d.λ)

Base.minimum(::ContinuousPoisson) = 0
Base.maximum(::ContinuousPoisson) = Inf

Distributions.cdf(d::ContinuousPoisson, x::Real) = StatsFuns.gammaccdf(x + 1, 1, d.λ)

function Distributions.logcdf(d::ContinuousPoisson, x::Real)
    return StatsFuns.gammalogccdf(x + 1, 1, d.λ)
end

struct ContinuousBinomial{T<:Real,B<:ADTypes.AbstractADType} <: ContinuousUnivariateDistribution
    n::Int
    p::T
    ad_type::B
end

ContinuousBinomial(n, p) = ContinuousBinomial(n, p, ADTypes.AutoFiniteDiff())

_discretize(d::ContinuousBinomial) = Binomial(d.n, d.p)

Base.minimum(::ContinuousBinomial) = 0
Base.maximum(::ContinuousBinomial) = d.n

function Distributions.cdf(d::ContinuousBinomial, x::Real)
    return StatsFuns.betaccdf(min(x + 1, d.n), max(d.n - x, 0), d.p)
end

function Distributions.logcdf(d::ContinuousBinomial, x::Real)
    return StatsFuns.betalogccdf(min(x + 1, d.n), max(d.n - x, 0), d.p)
end

struct ContinuousNegativeBinomial{T<:Real,B<:ADTypes.AbstractADType} <: ContinuousUnivariateDistribution
    r::Int
    p::T
    ad_type::B
end

ContinuousNegativeBinomial(r, p) = ContinuousNegativeBinomial(r, p, ADTypes.AutoFiniteDiff())

_discretize(d::ContinuousNegativeBinomial) = NegativeBinomial(d.r, d.p)

Base.minimum(::ContinuousNegativeBinomial) = 0
Base.maximum(::ContinuousNegativeBinomial) = Inf

function Distributions.cdf(d::ContinuousNegativeBinomial, x::Real)
    return StatsFuns.betacdf(d.r, x + 1, d.p)
end

function Distributions.logcdf(d::ContinuousNegativeBinomial, x::Real)
    return StatsFuns.betalogcdf(d.r, x + 1, d.p)
end

# shared code
for T in [:ContinuousPoisson, :ContinuousBinomial, :ContinuousNegativeBinomial]
    @eval begin
        function Distributions.quantile(d::$T, p::Real)
            upper = Distributions.quantile(_discretize(d), p)
            lower = upper - 1
            # NOTE: quantile_bisect is an internal function
            return Distributions.quantile_bisect(d, p, lower, upper)
        end

        function Distributions.pdf(d::$T, x::Real)
            return DI.derivative(Base.Fix1(Distributions.cdf, d), d.ad_type, x)
        end
        
        function Distributions.logpdf(d::$T, x::Real)
            return log(DI.derivative(Base.Fix1(Distributions.logcdf, d), d.ad_type, x)) + logcdf(d, x)
        end        
    end
end

```

And some quick plots to check that things look good:

```julia
using StatsPlots
dists = [ContinuousPoisson(5), ContinuousBinomial(10, 0.3), ContinuousNegativeBinomial(5, 0.3)]
plts = []
for dc in dists
    d = _discretize(dc)
    p1 = plot(d; func=cdf, label="discrete", title="CDF")
    plot!(p1, dc; func=cdf, label="continuous")
    p2 = plot(d; label="discrete", title="PDF")
    plot!(p2, dc; label="continuous")
    append!(plts, [p1, p2])
end
plot(plts...; layout=(3, 2), dpi=300)

```

 ![](https://global.discourse-cdn.com/julialang/original/3X/8/1/814e103df41ea333a1286cd7a05ac02b23db17d3.png)

Some major limitations are that the incomplete beta and gamma functions that are used in the CDFs are implemented in SpecialFunctions and restricted to `AbstractFloat` types, so they can only be differentiated by source-to-source ADs like Tapir or Enzyme. I’m not sure if either of these work on these functions at the moment, which is why these implementations use finite differences.

---

<div class="post-metadata">

**Author:** ![sdwfrost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sdwfrost/32/2831_2.png) [@sdwfrost](https://discourse.julialang.org/u/sdwfrost)\
**Post date:** [June 27, 2024, 11:57pm UTC](https://discourse.julialang.org/t/continuous-versions-of-poisson-binomial-negativebinomial/115781/3 "2024-06-27T23:57:21Z")

</div>

Thanks @sethaxen! I had to define `mode` for the distributions to get this working. It isn’t quite matching the PDFs I expect - the support for x\<0 should be 0. Here’s a comparison with the R library `cbinom` (which is based on [this paper](http://ac.inf.elte.hu/Vol_039_2013/137_39.pdf)  
[continuous\_distributions.jl](https://discourse.julialang.org/uploads/short-url/jguGAwJpUi9b9Ci3Pm2lIlK3Uy.jl) (2.8 KB)  
:

```Julia
include("./continuous_distributions.jl")
using RCall
@rlibrary cbinom

dc = ContinuousBinomial(10, 0.3)
x = collect(0:0.1:15)
y = pdf(dc,x)
yr = R"cbinom::dcbinom($x, 10, 0.3, log=FALSE)"
yr = convert(Array,yr)

plot(x,y,color="blue",label="Julia ContinuousBinomial")
plot!(x,yr,color="red",label="R cbinom")

```

![cbinom](https://global.discourse-cdn.com/julialang/original/3X/a/8/a826d4e0b2914150367f600bdf909b98191bd8be.png)

I’ll go through the fine print - many thanks for the pointers to the AD aspects of defining distributions, this was something that I had no idea how to get started!

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [June 28, 2024, 6:22am UTC](https://discourse.julialang.org/t/continuous-versions-of-poisson-binomial-negativebinomial/115781/4 "2024-06-28T06:22:54Z")

</div>

> [@sdwfrost](#):
>
> It isn’t quite matching the PDFs I expect - the support for x\<0 should be 0.

At first glance it looks like an off-by-one error?

---

<div class="post-metadata">

**Author:** ![sdwfrost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sdwfrost/32/2831_2.png) [@sdwfrost](https://discourse.julialang.org/u/sdwfrost)\
**Post date:** [June 28, 2024, 5:26pm UTC](https://discourse.julialang.org/t/continuous-versions-of-poisson-binomial-negativebinomial/115781/5 "2024-06-28T17:26:03Z")

</div>

Yes, that’s right; I was just getting my head around the parameterizations. The R package uses the [Ilienko (2013)](https://search.r-project.org/CRAN/refmans/cbinom/html/cbinom.html) approach, while the solution @sethaxen gave used the approach of [Padellini and Rue](https://repository.kaust.edu.sa/server/api/core/bitstreams/4bd7c7bb-aee4-4bad-ad48-08562484880c/content) to interpolate the CDFs. The latter approach gives a lower domain of -1 rather than 0. Rather than re-write the distributions, I think I’ll just have a note to myself to remember that generating random numbers from these can lead to negative values. I had to fix the mode calculations too; the corrected version is attached.  
[continuous\_distributions.jl](https://discourse.julialang.org/uploads/short-url/b3krX4j1nuofi61cGHRyfgZbIr0.jl) (2.8 KB)

---

<div class="post-metadata">

**Author:** ![sethaxen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sethaxen/32/35604_2.png) [@sethaxen](https://discourse.julialang.org/u/sethaxen)\
**Post date:** [July 4, 2024, 8:44am UTC](https://discourse.julialang.org/t/continuous-versions-of-poisson-binomial-negativebinomial/115781/6 "2024-07-04T08:44:08Z")

</div>

Yeah so the issue is that if we’re matching the CDFs at all points in the support of the discrete distribution, then we need to have nonzero probability mass at x \le 0. We can do that by adding negative numbers to the support or by adding a discrete “atom” at x=0 (e.g. how `Distributions.censored` works). But since that adds a discrete jump in the “density”, plots will look a little strange, and you’d need to modify `logpdf` to only use AD when x\>0. Shifting everything by 1 makes the CDFs no longer match, but in some cases makes the PDFs might match better.