I am looking for a way to simulate draws from a negative binomial distribution for a computational experiment on biological sequencing data. I am using a high performance package which only has certain distributions however, and though I know that gamma+poisson draws would give me the required simulation, the package lacks the latter. Here is a list of the distributions I have: [uniform, triangular, beta, exponential, gaussian, lognormal, vonmises, weibull, pareto, gamma]. Is there a simple approximation I can use from any of these functions that would require the least computational complexity (i.e. dont loop over creating many random draws for a single sample).
7
-
$\begingroup$ All the distribution you have are continuous. I think you need at least one discrete distribution because the negative binomial distribution is a discrete distribution. $\endgroup$– Abdoul HakiCommented Aug 2, 2021 at 20:49
-
$\begingroup$ Yes, but I don't have access to any. Although I do have access to a random choice from a list which would simulate a uniform over a discrete support. I am also not looking for an exact simulation, just similar properties and reasonably similar distribution. One idea was using exponential wait times to try to simulate Poisson, but I think the repetitions are excessive and not a great solution for the problem. $\endgroup$– GTOgodCommented Aug 2, 2021 at 21:13
-
$\begingroup$ What do you mean "but I don't have access to any" ? There are many, many free software packages, sites, and platforms you could use. If you are on a proprietary HPC platform with no external connections then you could just use C, Fortran, C++ or whatever other compilers are available. Sampling from distributions should be easy in any general purpose language. $\endgroup$– Robert LongCommented Aug 2, 2021 at 21:59
-
$\begingroup$ Another option could be to make the draws from the distribution you need in advance, and then access them from your HPC platform. $\endgroup$– Robert LongCommented Aug 2, 2021 at 22:05
-
$\begingroup$ I am using a threaded pipeline in a specialized bioinformatics language which doesn't allow me to use external code while threading. You are correct I could make the draws beforehand, but I would need millions of draws -- storing these as a list beforehand is a pretty clunky solution. $\endgroup$– GTOgodCommented Aug 2, 2021 at 22:21
|
Show 2 more comments
3 Answers
$\begingroup$
$\endgroup$
1
Here is a screenshot of a reference by NL Johnson et al. (2005) summarising the possible approaches based on a Geometric sampling:
Devroye (1985) indeed does not cover any other algorithm than the Gamma-Poisson version. Note that generating a Poisson variate is feasible from the Exponential generator, due to the connection between the Poisson process and the Exponential distribution. It is however not the recommended algorithm.
-
1$\begingroup$ The documentation for
rnbinominRsays it uses Devroye's (1986) method as a "Gamma mixture of Poissons." It looks like one really needs a way to draw Gamma variates. Ahrens & Dieter (1982) published a rejection method in CACM 25, 47-54. $\endgroup$– whuber ♦Commented Aug 3, 2021 at 18:04
$\begingroup$
$\endgroup$
One possible way is to generate a large number of value with frequencies their probabilities and then sample from a discrete uniform choice without replacement. By example, suppose your probability of success is $p=1/2$ and the number of failure until the experience has stop is $r=5$.
Then you start by computing the probability of each value using the negative binomial distribution. For the example, you will have $$P(X=0) = 0.03125, P(X=1) = 0.07813, P(X=2)=0.11719, \dots$$
Now, you create a large list of frequency corresponding to those probability. If you want to create a vector of length 100000, then this vector will have 3125 times $0$, 7813 times $1$; 11719 times $2$ and so one.
After, you just choose without replacement from that list the number $n$ of trials you want.
$\begingroup$
$\endgroup$
Why don't you use the inverse CDF method? The CDF of the negative binomial is the regularized incomplete beta function, which is just a ratio of beta functions so it should not be a problem to obtain it with your setup. You can see the detailshere and here. If you want something more fancy tailored specifically to the negative binomial distribution, the standard reference is Devroye's book Non-uniform random variate generation.
