The R Project SVN R

Rev

Rev 74265 | Rev 79365 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
42333 ripley 1
% File src/library/stats/man/GammaDist.Rd
68948 ripley 2
% Part of the R package, https://www.R-project.org
74265 hornik 3
% Copyright 1995-2018 R Core Team
42333 ripley 4
% Distributed under GPL 2 or later
5
 
27442 ripley 6
\name{GammaDist}
56186 murdoch 7
\alias{GammaDist}
27442 ripley 8
\alias{dgamma}
9
\alias{pgamma}
10
\alias{qgamma}
11
\alias{rgamma}
28897 maechler 12
\concept{incomplete gamma function}
27442 ripley 13
\title{The Gamma Distribution}
14
\description{
15
  Density, distribution function, quantile function and random
16
  generation for the Gamma distribution with parameters \code{shape} and
17
  \code{scale}.
18
}
19
\usage{
20
dgamma(x, shape, rate = 1, scale = 1/rate, log = FALSE)
30912 ripley 21
pgamma(q, shape, rate = 1, scale = 1/rate, lower.tail = TRUE,
22
       log.p = FALSE)
23
qgamma(p, shape, rate = 1, scale = 1/rate, lower.tail = TRUE,
24
       log.p = FALSE)
27442 ripley 25
rgamma(n, shape, rate = 1, scale = 1/rate)
26
}
27
\arguments{
28
  \item{x, q}{vector of quantiles.}
29
  \item{p}{vector of probabilities.}
30
  \item{n}{number of observations. If \code{length(n) > 1}, the length
31
    is taken to be the number required.}
32
  \item{rate}{an alternative way to specify the scale.}
43895 maechler 33
  \item{shape, scale}{shape and scale parameters.  Must be positive,
34
    \code{scale} strictly.}
36220 ripley 35
  \item{log, log.p}{logical; if \code{TRUE}, probabilities/densities \eqn{p}
36
    are returned as \eqn{log(p)}.}
27442 ripley 37
  \item{lower.tail}{logical; if TRUE (default), probabilities are
50805 ripley 38
    \eqn{P[X \le x]}, otherwise, \eqn{P[X > x]}.}
27442 ripley 39
}
40
\value{
41
  \code{dgamma} gives the density,
34848 ripley 42
  \code{pgamma} gives the distribution function,
27442 ripley 43
  \code{qgamma} gives the quantile function, and
44
  \code{rgamma} generates random deviates.
37274 ripley 45
 
46
  Invalid arguments will result in return value \code{NaN}, with a warning.
61975 murdoch 47
 
48
  The length of the result is determined by \code{n} for
49
  \code{rgamma}, and is the maximum of the lengths of the
65782 ripley 50
  numerical arguments for the other functions.  
61975 murdoch 51
 
65782 ripley 52
  The numerical arguments other than \code{n} are recycled to the
61975 murdoch 53
  length of the result.  Only the first elements of the logical
65782 ripley 54
  arguments are used.
27442 ripley 55
}
56
\details{
57
  If \code{scale} is omitted, it assumes the default value of \code{1}.
58
 
59
  The Gamma distribution with parameters \code{shape} \eqn{=\alpha}{= a}
60
  and \code{scale} \eqn{=\sigma}{= s} has density
61
  \deqn{
62
    f(x)= \frac{1}{{\sigma}^{\alpha}\Gamma(\alpha)} {x}^{\alpha-1} e^{-x/\sigma}%
63
  }{f(x)= 1/(s^a Gamma(a)) x^(a-1) e^-(x/s)}
50805 ripley 64
  for \eqn{x \ge 0}, \eqn{\alpha > 0}{a > 0} and \eqn{\sigma > 0}{s > 0}.
34851 ripley 65
  (Here \eqn{\Gamma(\alpha)}{Gamma(a)} is the function implemented by \R's
59979 ripley 66
  \code{\link{gamma}()} and defined in its help.  Note that \eqn{a = 0}
43899 maechler 67
  corresponds to the trivial distribution with all mass at point 0.)
34848 ripley 68
 
27442 ripley 69
  The mean and variance are
70
  \eqn{E(X) = \alpha\sigma}{E(X) = a*s} and
71
  \eqn{Var(X) = \alpha\sigma^2}{Var(X) = a*s^2}.
28231 maechler 72
 
34848 ripley 73
  The cumulative hazard \eqn{H(t) = - \log(1 - F(t))}{H(t) = - log(1 - F(t))}
61173 ripley 74
  is
75
\preformatted{-pgamma(t, ..., lower = FALSE, log = TRUE)
76
}
39015 ripley 77
 
46325 ripley 78
  Note that for smallish values of \code{shape} (and moderate
79
  \code{scale}) a large parts of the mass of the Gamma distribution is
80
  on values of \eqn{x} so near zero that they will be represented as
61177 ripley 81
  zero in computer arithmetic.  So \code{rgamma} may well return values
46325 ripley 82
  which will be represented as zero.  (This will also happen for very
83
  large values of \code{scale} since the actual generation is done for
61160 ripley 84
  \code{scale = 1}.)
50414 ripley 85
}
50399 maechler 86
% Have caught all currently known problems; hence no longer say:
87
%   Similarly, \code{qgamma} has a very hard job for
88
%   small \code{scale}, and warns of potential unreliability for
89
%   \code{scale < 1e-10}.
27442 ripley 90
\note{
74363 maechler 91
  The S (Becker \emph{et al}, 1988) parametrization was via \code{shape}
71883 ripley 92
  and \code{rate}: S had no \code{scale} parameter. It is an error
93
  to supply and \code{scale} and \code{rate}.
28231 maechler 94
 
28789 ripley 95
  \code{pgamma} is closely related to the incomplete gamma function.  As
42706 ripley 96
  defined by Abramowitz and Stegun 6.5.1 (and by \sQuote{Numerical
97
    Recipes}) this is
61177 ripley 98
  \deqn{P(a,x) = \frac{1}{\Gamma(a)} \int_0^x t^{a-1} e^{-t} dt}{P(a,x) = 1/Gamma(a) integral_0^x t^(a-1) exp(-t) dt}
28789 ripley 99
  \eqn{P(a, x)} is \code{pgamma(x, a)}.  Other authors (for example
100
  Karl Pearson in his 1922 tables) omit the normalizing factor,
61042 maechler 101
  defining the incomplete gamma function \eqn{\gamma(a,x)} as
102
  \eqn{\gamma(a,x) = \int_0^x t^{a-1} e^{-t} dt,}{gamma(a,x) =
103
    integral_0^x t^(a-1) exp(-t) dt,} i.e., \code{pgamma(x, a) * gamma(a)}.
104
  Yet other use the \sQuote{upper} incomplete gamma function,
61177 ripley 105
  \deqn{\Gamma(a,x) = \int_x^\infty t^{a-1} e^{-t} dt,}{Gamma(a,x) = integral_x^Inf t^(a-1) exp(-t) dt,}
61042 maechler 106
  which can be computed by
61160 ripley 107
  \code{pgamma(x, a, lower = FALSE) * gamma(a)}.
61042 maechler 108
 
109
  Note however that \code{pgamma(x, a, ..)} currently requires \eqn{a > 0},
110
  whereas the incomplete gamma function is also defined for negative
111
  \eqn{a}.  In that case, you can use \code{gamma_inc(a,x)} (for
62120 hornik 112
  \eqn{\Gamma(a,x)}) from package \CRANpkg{gsl}.
61042 maechler 113
 
114
  See also
68939 ripley 115
  \url{https://en.wikipedia.org/wiki/Incomplete_gamma_function}, or
61042 maechler 116
  \url{http://dlmf.nist.gov/8.2#i}.
37222 ripley 117
}
118
\source{
37226 ripley 119
  \code{dgamma} is computed via the Poisson density, using code contributed
37222 ripley 120
  by Catherine Loader (see \code{\link{dbinom}}).
43895 maechler 121
 
42987 ripley 122
  \code{pgamma} uses an unpublished (and not otherwise documented)
123
  algorithm \sQuote{mainly by Morten Welinder}.
34848 ripley 124
 
37222 ripley 125
  \code{qgamma} is based on a C translation of
126
 
127
  Best, D. J. and D. E. Roberts (1975).
128
  Algorithm AS91. Percentage points of the chi-squared distribution.
129
  \emph{Applied Statistics},  \bold{24}, 385--388.
130
 
131
  plus a final Newton step to improve the approximation.
43895 maechler 132
 
37222 ripley 133
  \code{rgamma} for \code{shape >= 1} uses
134
 
135
  Ahrens, J. H. and Dieter, U. (1982).
136
  Generating gamma variates by a modified rejection technique.
137
  \emph{Communications of the ACM}, \bold{25}, 47--54,
138
 
139
  and for \code{0 < shape < 1} uses
140
 
141
  Ahrens, J. H. and Dieter, U. (1974).
142
  Computer methods for sampling from gamma, beta, Poisson and binomial
143
  distributions. \emph{Computing}, \bold{12}, 223--246.
27442 ripley 144
}
145
\references{
74265 hornik 146
  Becker, R. A., Chambers, J. M. and Wilks, A. R. (1988).
27442 ripley 147
  \emph{The New S Language}.
47262 ripley 148
  Wadsworth & Brooks/Cole.
28231 maechler 149
 
74265 hornik 150
  Shea, B. L. (1988).
151
  Algorithm AS 239: Chi-squared and incomplete Gamma integral,
152
  \emph{Applied Statistics (JRSS C)}, \bold{37}, 466--473.
153
  \doi{10.2307/2347328}.
28789 ripley 154
 
155
  Abramowitz, M. and Stegun, I. A. (1972)
156
  \emph{Handbook of Mathematical Functions.} New York: Dover.
157
  Chapter 6: Gamma and Related Functions.
61042 maechler 158
 
159
  NIST Digital Library of Mathematical Functions.
160
  \url{http://dlmf.nist.gov/}, section 8.2.
27442 ripley 161
}
162
\seealso{
52770 ripley 163
  \code{\link{gamma}} for the gamma function.
164
 
165
  \link{Distributions} for other standard distributions, including
166
  \code{\link{dbeta}} for the Beta distribution and \code{\link{dchisq}}
167
  for the chi-squared distribution which is a special case of the Gamma
168
  distribution.
27442 ripley 169
}
170
\examples{
61160 ripley 171
-log(dgamma(1:4, shape = 1))
27442 ripley 172
p <- (1:9)/10
61168 ripley 173
pgamma(qgamma(p, shape = 2), shape = 2)
61160 ripley 174
1 - 1/exp(qgamma(p, shape = 1))
39015 ripley 175
 
55251 ripley 176
\donttest{# even for shape = 0.001 about half the mass is on numbers
39015 ripley 177
# that cannot be represented accurately (and most of those as zero)
178
pgamma(.Machine$double.xmin, 0.001)
52417 maechler 179
pgamma(5e-324, 0.001)  # on most machines 5e-324 is the smallest
39015 ripley 180
                       # representable non-zero number
181
table(rgamma(1e4, 0.001) == 0)/1e4
55251 ripley 182
}}
27442 ripley 183
\keyword{distribution}