{"id":8063,"date":"2014-11-15T08:27:52","date_gmt":"2014-11-15T07:27:52","guid":{"rendered":"http:\/\/djalil.chafai.net\/blog\/?p=8063"},"modified":"2014-11-17T11:14:32","modified_gmt":"2014-11-17T10:14:32","slug":"about-random-generators-of-geometric-distribution","status":"publish","type":"post","link":"https:\/\/djalil.chafai.net\/blog\/2014\/11\/15\/about-random-generators-of-geometric-distribution\/","title":{"rendered":"About random generators of geometric distribution"},"content":{"rendered":"

\"Dice\"<\/a><\/p>\n

Recently, a friend of mine, Arnaud Guyader<\/a>, discovered by accident during a teaching session that the geometric distribution generator of GNU R<\/a>, implemented in the rgeom<\/em> function, is sub-optimal. The command help rgeom<\/em> says rgeom uses the derivation as an exponential mixture of Poissons, see Devroye, L. (1986) Non-Uniform Random Variate Generation. Springer-Verlag, New York. Page 480.<\/em> Moreover, page 480 of this book says nothing about the geometric distribution.<\/p>\n

\/* Source code of rgeom.c in latest version of GNU R (2014-11-15).
\n*\u00a0 Mathlib : A C Library of Special Functions
\n*\u00a0 Copyright (C) 1998 Ross Ihaka and the R Core Team.
\n*\u00a0 Copyright (C) 2000 The R Core Team
\n*
\n*\u00a0 This program is free software; you can redistribute it and\/or modify
\n*\u00a0 it under the terms of the GNU General Public License as published by
\n*\u00a0 the Free Software Foundation; either version 2 of the License, or
\n*\u00a0 (at your option) any later version.
\n*
\n*\u00a0 This program is distributed in the hope that it will be useful,
\n*\u00a0 but WITHOUT ANY WARRANTY; without even the implied warranty of
\n*\u00a0 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.\u00a0 See the
\n*\u00a0 GNU General Public License for more details.
\n*
\n*\u00a0 You should have received a copy of the GNU General Public License
\n*\u00a0 along with this program; if not, a copy is available at
\n*\u00a0 http:\/\/www.r-project.org\/Licenses\/
\n*
\n*\u00a0 SYNOPSIS
\n*
\n*\u00a0\u00a0\u00a0 #include <Rmath.h>
\n*\u00a0\u00a0\u00a0 double rgeom(double p);
\n*
\n*\u00a0 DESCRIPTION
\n*
\n*\u00a0\u00a0\u00a0 Random variates from the geometric distribution.
\n*
\n*\u00a0 NOTES
\n*
\n*\u00a0\u00a0\u00a0 We generate lambda as exponential with scale parameter
\n*\u00a0\u00a0\u00a0 p \/ (1 - p).\u00a0 Return a Poisson deviate with mean lambda.
\n*
\n*\u00a0 REFERENCE
\n*
\n*\u00a0\u00a0\u00a0 Devroye, L. (1986).
\n*\u00a0\u00a0\u00a0 Non-Uniform Random Variate Generation.
\n*\u00a0\u00a0\u00a0 New York: Springer-Verlag.
\n*\u00a0\u00a0\u00a0 Page 480.
\n*\/<\/p>\n

#include \"nmath.h\"<\/p>\n

double rgeom(double p)
\n{
\nif (!R_FINITE(p) || p <= 0 || p > 1) ML_ERR_return_NAN;<\/p>\n

return rpois(exp_rand() * ((1 - p) \/ p));
\n}<\/p><\/blockquote>\n

The fact that the geometric distribution is an exponential mixture of Poisson is a special case of the more general fact that the negative binomial distribution, i.e. the convolution power of the geometric distribution also known as the Pascal distribution, is a gamma mixture of Poisson. Namely if \\( {X\\sim\\mathrm{Gamma}(n,\\ell)} \\) with \\( {\\ell=p\/(1-p)} \\), and if \\( {Y|X\\sim\\mathrm{Poisson}(X)} \\) then \\( {Y\\sim\\mathrm{NegativeBinomial}(n,p)} \\), since for every \\( {k\\geq n} \\), \\( {\\mathbb{P}(Y=k)} \\) writes<\/p>\n

\\[ \\int_0^\\infty\\!e^{-\\lambda} \\frac{\\lambda^k}{k!} \\frac{\\ell^n\\lambda^{n-1}}{\\Gamma(n)}e^{-\\lambda\\ell}\\,d\\lambda =\\frac{\\ell^n}{k!\\Gamma(n)} \\int_0^\\infty\\!\\lambda^{n+k-1}e^{-\\lambda(\\ell+1)}\\,d\\lambda =\\frac{\\Gamma(n+k)}{k!\\Gamma(n)}(1-p)^kp^n. \\]<\/p>\n

For \\( {n=1} \\) we recover the geometric distribution. This mixture representation is useful for the simulation of negative binomial distributions with high integer parameter, but is not the best choice for the simulation of the special case of the geometric distribution.<\/p>\n

The simulation of the geometric distribution is better done by using the so-called inversion method (generalized inverse of the cdf), which boils down to take the integer part of the standard exponential distribution generator: ceil(log(rand())\/log(1-p))<\/em> since if \\( {X\\sim\\mathrm{Exponential}(\\lambda)} \\) with \\( {\\lambda=-\\log(1-p)} \\) then<\/p>\n

\\[ \\mathbb{P}(k\\leq X\\leq k+1) =e^{-k\\lambda}-e^{-(k+1)\\lambda} =e^{-\\lambda k}(1-e^{-\\lambda}) =(1-p)^kp. \\]<\/p>\n

Patch. <\/b>Here is a possible patch<\/a> to the rgeom.c of GNU R.<\/p>\n

\n

31,32c31
\n<\u00a0 *\u00a0\u00a0\u00a0 We generate lambda as exponential with scale parameter
\n<\u00a0 *\u00a0\u00a0\u00a0 p \/ (1 - p).\u00a0 Return a Poisson deviate with mean lambda.
\n---
\n>\u00a0 *\u00a0\u00a0\u00a0 returns 1+ integer part of exponential with scale parameter -1 \/ log(1 - p).
\n39c38
\n<\u00a0 *\u00a0\u00a0\u00a0 Page 480.
\n---
\n>\u00a0 *\u00a0\u00a0\u00a0 Section X.2 page 498.
\n48c47
\n<\u00a0\u00a0\u00a0\u00a0 return rpois(exp_rand() * ((1 - p) \/ p));
\n---
\n>\u00a0\u00a0\u00a0\u00a0 return ceil(log(rand()) \/ log(1 - p));<\/p>\n<\/blockquote>\n

Further reading.<\/b> The simulation of the geometric distribution is discussed in Section X.2 page 498 of the book of Devroye: the inversion method is considered, but no words on the exponential mixture of Poisson. The simulation of the negative binomial distribution is discussed in Section X.4.7 page 543, and the Gamma mixture of Poisson is considered.<\/p>\n

Other software.<\/b> The random generator for the geometric distribution proposed presently by the GNU Scientific Library<\/a>, by the GNU Octave<\/a>, and by the Scilab<\/a> software are based on the inversion method. The Distributions.jl<\/a> package for the Julia<\/a> software also uses the inversion method. The good news with open source software is that even if the help does not mention the algorithm that is used, you can always check the source code.<\/p>\n","protected":false},"excerpt":{"rendered":"

Recently, a friend of mine, Arnaud Guyader, discovered by accident during a teaching session that the geometric distribution generator of GNU R, implemented in the…<\/p>\n

Continue readingAbout random generators of geometric distribution<\/span><\/a><\/div>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"iawp_total_views":112},"categories":[1],"tags":[],"_links":{"self":[{"href":"https:\/\/djalil.chafai.net\/blog\/wp-json\/wp\/v2\/posts\/8063"}],"collection":[{"href":"https:\/\/djalil.chafai.net\/blog\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/djalil.chafai.net\/blog\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/djalil.chafai.net\/blog\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/djalil.chafai.net\/blog\/wp-json\/wp\/v2\/comments?post=8063"}],"version-history":[{"count":34,"href":"https:\/\/djalil.chafai.net\/blog\/wp-json\/wp\/v2\/posts\/8063\/revisions"}],"predecessor-version":[{"id":8133,"href":"https:\/\/djalil.chafai.net\/blog\/wp-json\/wp\/v2\/posts\/8063\/revisions\/8133"}],"wp:attachment":[{"href":"https:\/\/djalil.chafai.net\/blog\/wp-json\/wp\/v2\/media?parent=8063"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/djalil.chafai.net\/blog\/wp-json\/wp\/v2\/categories?post=8063"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/djalil.chafai.net\/blog\/wp-json\/wp\/v2\/tags?post=8063"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}