Estimate Parameters of Gamma Distribution

Description

Estimate the shape and scale parameters (or the mean and coefficient of variation) of a Gamma distribution.

Usage

1
2
3
4
5
6
7
  egamma(x, method = "mle", ci = FALSE, 
    ci.type = "two-sided", ci.method = "normal.approx", 
    normal.approx.transform = "kulkarni.powar", conf.level = 0.95)

  egammaAlt(x, method = "mle", ci = FALSE, 
    ci.type = "two-sided", ci.method = "normal.approx", 
    normal.approx.transform = "kulkarni.powar", conf.level = 0.95)

Arguments

x

numeric vector of non-negative observations. Missing (NA), undefined (NaN), and infinite (Inf, -Inf) values are allowed but will be removed.

method

character string specifying the method of estimation. The possible values are:
"mle" (maximum likelihood; the default),
"bcmle" (bias-corrected mle),
"mme" (method of moments), and
"mmue" (method of moments based on the unbiased estimator of variance).
See the DETAILS section for more information.

ci

logical scalar indicating whether to compute a confidence interval for the mean. The default value is ci=FALSE.

ci.type

character string indicating what kind of confidence interval to compute. The possible values are "two-sided" (the default), "lower", and "upper". This argument is ignored if ci=FALSE.

ci.method

character string indicating which method to use to construct the confidence interval. Possible values are "normal.approx" (the default), "profile.likelihood", and "chisq.approx". This argument is ignored if ci=FALSE.

normal.approx.transform

character string indicating which power transformation to use when
ci.method="normal.approx". Possible values are
"kulkarni.powar" (the default), "cube.root", and "fourth.root". See the DETAILS section for more informaiton. This argument is ignored if ci=FALSE or ci.method="chisq.approx".

conf.level

a scalar between 0 and 1 indicating the confidence level of the confidence interval. The default value is conf.level=0.95. This argument is ignored if ci=FALSE.

Details

If x contains any missing (NA), undefined (NaN) or infinite (Inf, -Inf) values, they will be removed prior to performing the estimation.

Let \underline{x} = x_1, x_2, …, x_n denote a random sample of n observations from a gamma distribution with parameters shape=α and scale=β. The relationship between these parameters and the mean (mean=μ) and coefficient of variation (cv=τ) of this distribution is given by:

α = τ^{-2} \;\;\;\;\;\; (1)

β = μ/α \;\;\;\;\;\; (2)

μ = α \; β \;\;\;\;\;\; (3)

τ = α^{-1/2} \;\;\;\;\;\; (4)

The function egamma returns estimates of the shape and scale parameters. The function egammaAlt returns estimates of the mean (μ) and coefficient of variation (cv) based on the estimates of the shape and scale parameters.

Estimation

Maximum Likelihood Estimation (method="mle")
The maximum likelihood estimators (mle's) of the shape and scale parameters α and β are solutions of the simultaneous equations:

\hat{α}_{mle} = \frac{1}{n}∑_{i=1}^n log(x_i) - log(\bar{x}) = ψ(\hat{α}_{mle}) - log(\hat{α}_{mle}) \ \;\;\;\;\;\; (5)

\hat{β}_{mle} = \bar{x} / \hat{α} \;\;\;\;\;\; (6)

where ψ denotes the digamma function, and \bar{x} denotes the sample mean:

\bar{x} = \frac{1}{n}∑_{i=1}^n x_i \;\;\;\;\;\; (7)

(Forbes et al., 2011, chapter 22; Johnson et al., 1994, chapter 17).

Bias-Corrected Maximum Likelihood Estimation (method="bcmle")
The “bias-corrected” maximum likelihood estimator of the shape parameter is based on the suggestion of Anderson and Ray (1975; see also Johnon et al., 1994, p.366 and Singh et al., 2010b, p.48), who noted that the bias of the maximum likelihood estimator of the shape parameter can be considerable when the sample size is small. This estimator is given by:

\hat{α}_{bcmle} = \frac{n-3}{n}\hat{α}_{mle} + \frac{2}{3n} \;\;\;\;\;\; (8)

The estimate of the scale paramter is not modified (i.e., the mle of β is returned).

Method of Moments Estimation (method="mme")
The method of moments estimators (mme's) of the shape and scale parameters α and β are:

\hat{α}_{mme} = (\bar{x}/s_m)^2 \;\;\;\;\;\; (9)

\hat{β}_{mme} = s_m^2 / \bar{x} \;\;\;\;\;\; (10)

where s_m^2 denotes the method of moments estimator of variance:

s_m^2 = \frac{1}{n} ∑_{i=1}^n (x_i - \bar{x})^2 \;\;\;\;\;\; (11)


Method of Moments Estimation Based on the Unbiased Estimator of Variance (method="mmue")
The method of moments estimators based on the unbiased estimator of variance are exactly the same as the method of moments estimators, except that the method of moments estimator of variance is replaced with the unbiased estimator of variance:

\hat{α}_{mmue} = (\bar{x}/s)^2 \;\;\;\;\;\; (12)

\hat{β}_{mmue} = s^2 / \bar{x} \;\;\;\;\;\; (13)

where s^2 denotes the unbiased estimator of variance:

s^2 = \frac{1}{n-1} ∑_{i=1}^n (x_i - \bar{x})^2 \;\;\;\;\;\; (14)


Confidence Intervals
This section discusses how confidence intervals for the mean μ are computed.

Normal Approximation (ci.method="normal.approx")
The normal approximation method is based on the method of Kulkarni and Powar (2010), who use a power transformation of the the original data to approximate a sample from a normal distribuiton, compute the confidence interval for the mean on the transformed scale using the usual formula for a confidence interval for the mean of a normal distribuiton, and then tranform the limits back to the original space using equations based on the expected value of a gamma random variable raised to a power.

The particular power used for the normal approximation is defined by the argument
normal.approx.transform. The value normal.approx.transform="cube.root" uses the cube root transformation suggested by Wilson and Hilferty (1931), and the value
"fourth.root" uses the fourth root transformation suggested by Hawkins and Wixley (1986). The default value "kulkarni.powar" uses the “Optimum Power Normal Approximation Method” of Kulkarni and Powar (2010), who show this method performs the best in terms of maintining coverage and minimizing confidence interval width compared to eight other methods. The “optimum” power p is determined by:

p = -0.0705 - 0.178\hat{α} + 0.475 √{\hat{α}} if \hat{α} ≤ 1.5
p = 0.246 if \hat{α} > 1.5 (15)

where \hat{α} denotes the estimate of the shape parameter. Kulkarni and Powar (2010) derived this equation by determining what power transformation yields a skew closest to 0 and a kurtosis closest to 3 for a gamma random variable with a given shape parameter. Although Kulkarni and Powar (2010) use the maximum likelihood estimate of shape to determine the power to use to induce approximate normality, for the functions egamma and egammaAlt the power is based on whatever estimate of shape is used (e.g., method="mle", method="bcmle", etc.).

Likelihood Profile (ci.method="profile.likelihood")
This method was proposed by Cox (1970, p.88), and Venzon and Moolgavkar (1988) introduced an efficient method of computation. This method is also discussed by Stryhn and Christensen (2003) and Royston (2007). The idea behind this method is to invert the likelihood-ratio test to obtain a confidence interval for the mean μ while treating the coefficient of variation τ as a nuisance parameter.

The likelihood function is given by:

L(θ, τ | \underline{x}) = ∏_{i=1}^n \frac{x_i^{α-1} e^{-x_i/β}}{β^α Γ(α)} \;\;\;\;\;\; (16)

where α and β are defined in Equations (1) and (2) above, and Γ(t) denotes the Gamma function evaluated at t.

Following Stryhn and Christensen (2003), denote the maximum likelihood estimates of the mean and coefficient of variation by (θ^*, τ^*). The likelihood ratio test statistic (G^2) of the hypothesis H_0: θ = θ_0 (where θ_0 is a fixed value) equals the drop in 2 log(L) between the “full” model and the reduced model with θ fixed at θ_0, i.e.,

G^2 = 2 \{log[L(θ^*, τ^*)] - log[L(θ_0, τ_0^*)]\} \;\;\;\;\;\; (17)

where τ_0^* is the maximum likelihood estimate of τ for the reduced model (i.e., when θ = θ_0). Under the null hypothesis, the test statistic G^2 follows a chi-squared distribution with 1 degree of freedom.

Alternatively, we may express the test statistic in terms of the profile likelihood function L_1 for the mean θ, which is obtained from the usual likelihood function by maximizing over the parameter τ, i.e.,

L_1(θ) = max_{τ} L(θ, τ) \;\;\;\;\;\; (18)

Then we have

G^2 = 2 \{log[L_1(θ^*)] - log[L_1(θ_0)]\} \;\;\;\;\;\; (19)

A two-sided (1-α)100\% confidence interval for the mean θ consists of all values of θ_0 for which the test is not significant at level alpha:

θ_0: G^2 ≤ χ^2_{1, {1-α}} \;\;\;\;\;\; (20)

where χ^2_{ν, p} denotes the p'th quantile of the chi-squared distribution with ν degrees of freedom.

Chi-Square Approximation (ci.method="chisq.approx")
This method is based on the relationship between the sample mean of the gamma distribution and the chi-squared distribution (Grice and Bain, 1980). Because this method is exact only when the shape parameter α is known, the method used here is called the “chi-square approximation” method because the estimate of the shape parameter is used. This method is not the method proposed by Grice and Bain (1980) in which the confidence interval is adjusted based on adjusting for the fact that the shape parameter is estimated. The chi-square approximation method used by egamma and egammaAlt is equivalent to the “Approximate gamma” method of Singh et al. (2010b, equation (2-34), p.51).

Value

a list of class "estimate" containing the estimated parameters and other information. See
estimate.object for details.

Warning

When ci=TRUE and ci.method="normal.approx", it is possible for the lower confidence limit based on the transformed data to be less than 0. In this case, the lower confidence limit on the original scale is set to 0 and a warning is issued stating that the normal approximation is not accurate in this case.

Note

The gamma distribution takes values on the positive real line. Special cases of the gamma are the exponential distribution and the chi-square distributions. Applications of the gamma include life testing, statistical ecology, queuing theory, inventory control, and precipitation processes. A gamma distribution starts to resemble a normal distribution as the shape parameter a tends to infinity.

Some EPA guidance documents (e.g., Singh et al., 2002; Singh et al., 2010a,b) strongly recommend against using a lognormal model for environmental data and recommend trying a gamma distribuiton instead.

Author(s)

Steven P. Millard (EnvStats@ProbStatInfo.com)

References

Anderson, C.W., and W.D. Ray. (1975). Improved Maximum Likelihood Estimators for the Gamma Distribution. Communications in Statistics, 4, 437–448.

Forbes, C., M. Evans, N. Hastings, and B. Peacock. (2011). Statistical Distributions, Fourth Edition. John Wiley and Sons, Hoboken, NJ.

Grice, J.V., and L.J. Bain. (1980). Inferences Concerning the Mean of the Gamma Distribution. Journal of the American Statistician, 75, 929-933.

Hawkins, D. M., and R.A.J. Wixley. (1986). A Note on the Transformation of Chi-Squared Variables to Normality. The American Statistician, 40, 296–298.

Johnson, N.L., S. Kotz, and N. Balakrishnan. (1994). Continuous Univariate Distributions, Volume 1. Second Edition. John Wiley and Sons, New York, Chapter 17.

Kulkarni, H.V., and S.K. Powar. (2010). A New Method for Interval Estimation of the Mean of the Gamma Distribution. Lifetime Data Analysis, 16, 431–447.

Singh, A., A.K. Singh, and R.J. Iaci. (2002). Estimation of the Exposure Point Concentration Term Using a Gamma Distribution. EPA/600/R-02/084. October 2002. Technology Support Center for Monitoring and Site Characterization, Office of Research and Development, Office of Solid Waste and Emergency Response, U.S. Environmental Protection Agency, Washington, D.C.

Singh, A., R. Maichle, and N. Armbya. (2010a). ProUCL Version 4.1.00 User Guide (Draft). EPA/600/R-07/041, May 2010. Office of Research and Development, U.S. Environmental Protection Agency, Washington, D.C.

Singh, A., N. Armbya, and A. Singh. (2010b). ProUCL Version 4.1.00 Technical Guide (Draft). EPA/600/R-07/041, May 2010. Office of Research and Development, U.S. Environmental Protection Agency, Washington, D.C.

Wilson, E.B., and M.M. Hilferty. (1931). The Distribution of Chi-Squares. Proceedings of the National Academy of Sciences, 17, 684–688.

See Also

GammaDist, estimate.object, eqgamma, predIntGamma, tolIntGamma.

Examples

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
  # Generate 20 observations from a gamma distribution with parameters 
  # shape=3 and scale=2, then estimate the parameters. 
  # (Note: the call to set.seed simply allows you to reproduce this 
  # example.)

  set.seed(250) 
  dat <- rgamma(20, shape = 3, scale = 2) 
  egamma(dat, ci = TRUE)

  #Results of Distribution Parameter Estimation
  #--------------------------------------------
  #
  #Assumed Distribution:            Gamma
  #
  #Estimated Parameter(s):          shape = 2.203862
  #                                 scale = 2.174928
  #
  #Estimation Method:               mle
  #
  #Data:                            dat
  #
  #Sample Size:                     20
  #
  #Confidence Interval for:         mean
  #
  #Confidence Interval Method:      Optimum Power Normal Approximation
  #                                 of Kulkarni & Powar (2010)
  #                                 using mle of 'shape'
  #
  #Normal Transform Power:          0.246
  #
  #Confidence Interval Type:        two-sided
  #
  #Confidence Level:                95%
  #
  #Confidence Interval:             LCL = 3.361652
  #                                 UCL = 6.746794

  # Clean up
  rm(dat)

  #----------------------------------------------------------------

  # Using the reference area TcCB data in EPA.94b.tccb.df, assume a 
  # gamma distribution, estimate the parameters based on the 
  # bias-corrected mle of shape, and compute a one-sided upper 90% 
  # confidence interval for the mean.

  attach(EPA.94b.tccb.df)

  #----------
  # First test to see whether the data appear to follow a gamma 
  # distribution.

  gofTest(TcCB[Area == "Reference"], dist = "gamma", 
    est.arg.list = list(method = "bcmle"))

  #Results of Goodness-of-Fit Test
  #-------------------------------
  #
  #Test Method:                     Shapiro-Wilk GOF Based on 
  #                                 Chen & Balakrisnan (1995)
  #
  #Hypothesized Distribution:       Gamma
  #
  #Estimated Parameter(s):          shape = 4.5695247
  #                                 scale = 0.1309788
  #
  #Estimation Method:               bcmle
  #
  #Data:                            TcCB[Area == "Reference"]
  #
  #Sample Size:                     47
  #
  #Test Statistic:                  W = 0.9703827
  #
  #Test Statistic Parameter:        n = 47
  #
  #P-value:                         0.2739512
  #
  #Alternative Hypothesis:          True cdf does not equal the
  #                                 Gamma Distribution.

  #----------
  # Now estimate the paramters and compute the upper confidence 
  # limit.

  egamma(TcCB[Area == "Reference"], method = "bcmle", ci = TRUE, 
    ci.type = "upper", conf.level = 0.9) 

  #Results of Distribution Parameter Estimation
  #--------------------------------------------
  #
  #Assumed Distribution:            Gamma
  #
  #Estimated Parameter(s):          shape = 4.5695247
  #                                 scale = 0.1309788
  #
  #Estimation Method:               bcmle
  #
  #Data:                            TcCB[Area == "Reference"]
  #
  #Sample Size:                     47
  #
  #Confidence Interval for:         mean
  #
  #Confidence Interval Method:      Optimum Power Normal Approximation
  #                                 of Kulkarni & Powar (2010)
  #                                 using bcmle of 'shape'
  #
  #Normal Transform Power:          0.246
  #
  #Confidence Interval Type:        upper
  #
  #Confidence Level:                90%
  #
  #Confidence Interval:             LCL = 0.0000000
  #                                 UCL = 0.6561838

  #----------------------------------------------------------------

  # Repeat the above example but use the alternative 
  # parameterization.

  egammaAlt(TcCB[Area == "Reference"], method = "bcmle", ci = TRUE, 
    ci.type = "upper", conf.level = 0.9) 

  #Results of Distribution Parameter Estimation
  #--------------------------------------------
  #
  #Assumed Distribution:            Gamma
  #
  #Estimated Parameter(s):          mean = 0.5985106
  #                                 cv   = 0.4678046
  #
  #Estimation Method:               bcmle of 'shape'
  #
  #Data:                            TcCB[Area == "Reference"]
  #
  #Sample Size:                     47
  #
  #Confidence Interval for:         mean
  #
  #Confidence Interval Method:      Optimum Power Normal Approximation
  #                                 of Kulkarni & Powar (2010)
  #                                 using bcmle of 'shape'
  #
  #Normal Transform Power:          0.246
  #
  #Confidence Interval Type:        upper
  #
  #Confidence Level:                90%
  #
  #Confidence Interval:             LCL = 0.0000000
  #                                 UCL = 0.6561838

  #----------

  # Clean up
  #---------
  detach("EPA.94b.tccb.df") 

Want to suggest features or report bugs for rdrr.io? Use the GitHub issue tracker.