Generating negative binomial counts using a Poisson mixture of Gamma distributions

Questions
,
1

Hello!

I’m trying to simulate counts of (pseudo-) species with different characteristics (rare, medium, abundant). From this point of view, it is quite natural to use Gamma-Poissson mixtures, with convenient parameters for the Gamma distributions.
It works very well for generating a SINGLE sample of each king of species, but my code completely fails in simulating several species :

Sampled median values of \[Gamma] : {0.000231012,0.000198812,0.000166503,0.000251062,0.000187207,0.000241421,0.000192554,0.000186855,0.000168828,0.000209722,0.000183093,4.94303,0.490392,0.490566,0.00018395,4.96136,0.000237141,0.488084,0.493871,0.00020586,4.95461,0.493009,4.92685,0.495393,0.000231126,4.93443,0.493136,0.496181,0.492077,0.491182}

Sampled median counts of species : {0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,5,0,0,0,0,5,0,5,0,0,5,0,0,0,0}

Theoretical median counts of species : {2,2,2,2,2,2,2,2,2,2,2,49967,4967,4967,2,49967,2,4967,4967,2,49967,4967,49967,4967,2,49967,4967,4967,4967,4967}

Here, theoretical median counts are of course right, sample median values of Gamma are correct, but sample median counts obtained from random values of Gamma are completely wrong! There is certainly a mistake somewhere in my code, but I’m not able to find it… Could somebody help?

PourMath.nb (192 KB)

2

I think there are potentially 3 things to change:

  1. P should be a number between 0 and 1.
  2. Fisher2Math should be Fisher2Math[k_, p_] := {k, (1 - p)/p}
  3. Grille should not include 0 and 1 and should be Grille = Range[0.05, 0.95, 0.05];

Then all of the means and medians match up. In short, to create a negative binomial distribution with parameters p and k using a mixture of Poisson and Gamma distributions, consider the following:

p = 1/100;
k = 0.2;
SeedRandom[12345];
r0 = RandomVariate[NegativeBinomialDistribution[k, p], 10^6];
r1 = RandomVariate[PoissonDistribution[#]] & /@ RandomVariate[GammaDistribution[k, (1 - p)/p], 10^6];
({#[r0], #[r1]} // N) & /@ {Mean, Median, Variance}
(* {{19.7711, 19.87}, {2., 2.}, {1977.72, 2000.76}} *)

I’m also not following why you’d want to generate random samples indirectly when you could do so directly. You certainly can do so. Maybe the gamma distribution parameters are later to be considered functions of environmental variables?