Message Boards Message Boards

Generating negative binomial counts using a Poisson mixture of Gamma distributions

Claude Mante, Retired
Posted 1 month ago

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?

Attachments:
PourMath.nb
POSTED BY: Claude Mante
1 Reply
Sort By:
Jim Baldwin
Jim Baldwin
Posted 17 days ago

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?
POSTED BY: Jim Baldwin
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Tag limit exceeded
Note: Only the first five people you tag will receive an email notification; the other tagged names will appear as links to their profiles.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract