Use native samplers for Poisson distribution

[devmotion]

I noticed that sampler is not implemented for Poisson distributions although two native samplers PoissonCountSampler and PoissonADSampler exist. Moreover, currently R is used for sampling from Poisson distributions, which is quite slow. For now I used the same cut-off value of 6 for both samplers, as chosen in PoissonRandom.jl.

I repeated the benchmarks

using Distributions, PoissonRandom, StatsFuns
using Plots

function n_count(rng, λ, n)
  tmp = 0
  for i in 1:n
    tmp += PoissonRandom.count_rand(rng,λ)
  end
  tmp
end

function n_pois(rng,λ,n)
  tmp = 0
  for i in 1:n
    tmp += pois_rand(rng,λ)
  end
  tmp
end

function n_ad(rng, λ, n)
  tmp = 0
  for i in 1:n
    tmp += PoissonRandom.ad_rand(rng, λ)
  end
  tmp
end

function n_dist(λ,n)
  tmp = 0
  for i in 1:n
    tmp += rand(Poisson(λ))
  end
  tmp
end

function n_rfunctions(λ, n)
  tmp = 0
  for i in 1:n
    tmp += convert(Int, StatsFuns.RFunctions.poisrand(λ))
  end
  tmp
end

function n_countsampler(rng, λ::Float64, n)
  tmp = 0
  for i in 1:n
    tmp += rand(rng, Distributions.PoissonCountSampler(λ))
  end
  tmp
end

function n_adsampler(rng, λ::Float64, n)
  tmp = 0
  for i in 1:n
    tmp += rand(rng, Distributions.PoissonADSampler(λ))
  end
  tmp
end

function time_λ!(rng, times, λ::Float64, n)
  times[1] = @elapsed n_count(rng, λ, n)
  times[2] = @elapsed n_ad(rng, λ, n)
  times[3] = @elapsed n_pois(rng, λ, n)
  times[4] = @elapsed n_dist(rng, λ, n)
  times[5] = @elapsed n_rfunctions(λ, n)
  times[6] = @elapsed n_countsampler(rng, λ, n)
  times[7] = @elapsed n_adsampler(rng, λ, n)

  nothing
end

function plot_benchmark(rng)
    times = Matrix{Float64}(undef, 7, 20)

    # Compile
    time_λ!(rng, view(times, :, 1), 5, 5_000_000)

    # Run with a bunch of λ
    for λ in 1:20
        time_λ!(rng, view(times, :, λ), float(λ), 5_000_000)
    end

    plot(times',
         labels = ["count_rand" "ad_rand" "pois_rand" "Distributions" "RFunctions" "PoissonCountSampler" "PoissonADSampler"],
         lw = 3)
end

from SciML/PoissonRandom.jl#6 with this PR. I get

using Random
Random.seed!(1234)
plot_benchmark(Random.GLOBAL_RNG)
savefig("global_rng.png")

and

using RandomNumbers
plot_benchmark(Xorshifts.Xoroshiro128Plus(1234))
savefig("xoroshiro128plus.png")

Using the native samplers leads to a significant speed-up and a performance which is on par with PoissonRandom.jl.

[codecov-io]

Codecov Report

Merging #1021 into master will increase coverage by 0.25%.
The diff coverage is 100%.

@@            Coverage Diff             @@
##           master    #1021      +/-   ##
==========================================
+ Coverage   77.92%   78.17%   +0.25%     
==========================================
  Files         112      112              
  Lines        5391     5325      -66     
==========================================
- Hits         4201     4163      -38     
+ Misses       1190     1162      -28

Impacted Files	Coverage Δ
src/samplers/poisson.jl	92.18% <100%> (+1.27%)	⬆️
src/univariate/discrete/poisson.jl	66.66% <100%> (+4.56%)	⬆️
src/multivariate/mvnormalcanon.jl	78.72% <0%> (-2.13%)	⬇️
src/multivariate/mvnormal.jl	71.59% <0%> (ø)	⬆️
src/univariate/discrete/discretenonparametric.jl	98.03% <0%> (+0.01%)	⬆️
src/utils.jl	80% <0%> (+6.53%)	⬆️

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 802a65b...ef80351. Read the comment docs.

[matbesancon]

Samplers are not my area, but do we have a way to test for correctness of the sampler?

[devmotion]

It seems they are tested in

Distributions.jl/test/samplers.jl

Lines 62 to 74 in 802a65b

	## Poisson samplers
	for (S, paramlst) in [
	(PoissonCountSampler, [0.2, 0.5, 1.0, 2.0, 5.0, 10.0, 15.0, 20.0, 30.0]),
	(PoissonADSampler, [5.0, 10.0, 15.0, 20.0, 30.0])]
	local S
	println(" testing $S")
	for μ in paramlst
	test_samples(S(μ), Poisson(μ), n_tsamples)
	test_samples(S(μ), Poisson(μ), n_tsamples, rng=rng)
	end
	end

[matbesancon]

LGTM, I'll wait for another review before merging though, to get more educated opinions

[mschauer]

You can just compare with 0?

[devmotion]

I guess that could be used here as well. I think it makes only a different if one works with non-standard number types such as unitful numbers - e.g., u"1.0m" > 0 errors whereas u"1.0m" > zero(u"1.0m") works as expected.

[mschauer]

Ah, good point, let's be defensive then.

[mschauer]

So this is the sub approach we will not use anymore?

[devmotion]

Yes, instead of the quantile-based generation this PR uses the PoissonCountSampler for small rates which just counts exponentially distributed random numbers.

[mschauer]

this magic threshold appears twice, make it a constant?

[devmotion]

I guess it's even better to not define rand at all and remove the code duplication.

[matbesancon]

bump on this? Maybe naming the constant is still a nice idea

[matbesancon]

@mschauer @devmotion good to merge?

[devmotion]

Yes 👍 I guess, I adjusted the PR according to your suggestions and comments

[matbesancon]

thanks for the PR :)