| // Copyright ©2016 The gonum Authors. All rights reserved. |
| // Use of this source code is governed by a BSD-style |
| // license that can be found in the LICENSE file. |
| |
| package distuv |
| |
| import ( |
| "math" |
| "math/rand" |
| |
| "gonum.org/v1/gonum/mathext" |
| ) |
| |
| // Gamma implements the Gamma distribution, a two-parameter continuous distribution |
| // with support over the positive real numbers. |
| // |
| // The gamma distribution has density function |
| // β^α / Γ(α) x^(α-1)e^(-βx) |
| // |
| // For more information, see https://en.wikipedia.org/wiki/Gamma_distribution |
| type Gamma struct { |
| // Alpha is the shape parameter of the distribution. Alpha must be greater |
| // than 0. If Alpha == 1, this is equivalent to an exponential distribution. |
| Alpha float64 |
| // Beta is the rate parameter of the distribution. Beta must be greater than 0. |
| // If Beta == 2, this is equivalent to a Chi-Squared distribution. |
| Beta float64 |
| |
| Source *rand.Rand |
| } |
| |
| // CDF computes the value of the cumulative distribution function at x. |
| func (g Gamma) CDF(x float64) float64 { |
| if x < 0 { |
| return 0 |
| } |
| return mathext.GammaInc(g.Alpha, g.Beta*x) |
| } |
| |
| // ExKurtosis returns the excess kurtosis of the distribution. |
| func (g Gamma) ExKurtosis() float64 { |
| return 6 / g.Alpha |
| } |
| |
| // LogProb computes the natural logarithm of the value of the probability |
| // density function at x. |
| func (g Gamma) LogProb(x float64) float64 { |
| if x <= 0 { |
| return math.Inf(-1) |
| } |
| a := g.Alpha |
| b := g.Beta |
| lg, _ := math.Lgamma(a) |
| return a*math.Log(b) - lg + (a-1)*math.Log(x) - b*x |
| } |
| |
| // Mean returns the mean of the probability distribution. |
| func (g Gamma) Mean() float64 { |
| return g.Alpha / g.Beta |
| } |
| |
| // Mode returns the mode of the normal distribution. |
| // |
| // The mode is NaN in the special case where the Alpha (shape) parameter |
| // is less than 1. |
| func (g Gamma) Mode() float64 { |
| if g.Alpha < 1 { |
| return math.NaN() |
| } |
| return (g.Alpha - 1) / g.Beta |
| } |
| |
| // NumParameters returns the number of parameters in the distribution. |
| func (Gamma) NumParameters() int { |
| return 2 |
| } |
| |
| // Prob computes the value of the probability density function at x. |
| func (g Gamma) Prob(x float64) float64 { |
| return math.Exp(g.LogProb(x)) |
| } |
| |
| // Quantile returns the inverse of the cumulative distribution function. |
| func (g Gamma) Quantile(p float64) float64 { |
| if p < 0 || p > 1 { |
| panic(badPercentile) |
| } |
| return mathext.GammaIncInv(g.Alpha, p) / g.Beta |
| } |
| |
| // Rand returns a random sample drawn from the distribution. |
| // |
| // Rand panics if either alpha or beta is <= 0. |
| func (g Gamma) Rand() float64 { |
| if g.Beta <= 0 { |
| panic("gamma: beta <= 0") |
| } |
| |
| unifrnd := rand.Float64 |
| exprnd := rand.ExpFloat64 |
| normrnd := rand.NormFloat64 |
| if g.Source != nil { |
| unifrnd = g.Source.Float64 |
| exprnd = g.Source.ExpFloat64 |
| normrnd = g.Source.NormFloat64 |
| } |
| |
| a := g.Alpha |
| b := g.Beta |
| switch { |
| case a <= 0: |
| panic("gamma: alpha < 0") |
| case a == 1: |
| // Generate from exponential |
| return exprnd() / b |
| case a < 0.3: |
| // Generate using |
| // Liu, Chuanhai, Martin, Ryan and Syring, Nick. "Simulating from a |
| // gamma distribution with small shape parameter" |
| // https://arxiv.org/abs/1302.1884 |
| // use this reference: http://link.springer.com/article/10.1007/s00180-016-0692-0 |
| |
| // Algorithm adjusted to work in log space as much as possible. |
| lambda := 1/a - 1 |
| lw := math.Log(a) - 1 - math.Log(1-a) |
| lr := -math.Log(1 + math.Exp(lw)) |
| lc, _ := math.Lgamma(a + 1) |
| for { |
| e := exprnd() |
| var z float64 |
| if e >= -lr { |
| z = e + lr |
| } else { |
| z = -exprnd() / lambda |
| } |
| lh := lc - z - math.Exp(-z/a) |
| var lEta float64 |
| if z >= 0 { |
| lEta = lc - z |
| } else { |
| lEta = lc + lw + math.Log(lambda) + lambda*z |
| } |
| if lh-lEta > -exprnd() { |
| return math.Exp(-z/a) / b |
| } |
| } |
| case a >= 0.3 && a < 1: |
| // Generate using: |
| // Kundu, Debasis, and Rameshwar D. Gupta. "A convenient way of generating |
| // gamma random variables using generalized exponential distribution." |
| // Computational Statistics & Data Analysis 51.6 (2007): 2796-2802. |
| |
| // TODO(btracey): Change to using Algorithm 3 if we can find the bug in |
| // the implementation below. |
| |
| // Algorithm 2. |
| alpha := g.Alpha |
| a := math.Pow(1-expNegOneHalf, alpha) / (math.Pow(1-expNegOneHalf, alpha) + alpha*math.Exp(-1)/math.Pow(2, alpha)) |
| b := math.Pow(1-expNegOneHalf, alpha) + alpha/math.E/math.Pow(2, alpha) |
| var x float64 |
| for { |
| u := unifrnd() |
| if u <= a { |
| x = -2 * math.Log(1-math.Pow(u*b, 1/alpha)) |
| } else { |
| x = -math.Log(math.Pow(2, alpha) / alpha * b * (1 - u)) |
| } |
| v := unifrnd() |
| if x <= 1 { |
| if v <= math.Pow(x, alpha-1)*math.Exp(-x/2)/(math.Pow(2, alpha-1)*math.Pow(1-math.Exp(-x/2), alpha-1)) { |
| break |
| } |
| } else { |
| if v <= math.Pow(x, alpha-1) { |
| break |
| } |
| } |
| } |
| return x / g.Beta |
| |
| /* |
| // Algorithm 3. |
| d := 1.0334 - 0.0766*math.Exp(2.2942*alpha) |
| a := math.Pow(2, alpha) * math.Pow(1-math.Exp(-d/2), alpha) |
| b := alpha * math.Pow(d, alpha-1) * math.Exp(-d) |
| c := a + b |
| var x float64 |
| for { |
| u := unifrnd() |
| if u <= a/(a+b) { |
| x = -2 * math.Log(1-math.Pow(c*u, 1/a)/2) |
| } else { |
| x = -math.Log(c * (1 - u) / (alpha * math.Pow(d, alpha-1))) |
| } |
| v := unifrnd() |
| if x <= d { |
| if v <= (math.Pow(x, alpha-1)*math.Exp(-x/2))/(math.Pow(2, alpha-1)*math.Pow(1-math.Exp(-x/2), alpha-1)) { |
| break |
| } |
| } else { |
| if v <= math.Pow(d/x, 1-alpha) { |
| break |
| } |
| } |
| } |
| return x / g.Beta |
| */ |
| case a > 1: |
| // Generate using: |
| // Marsaglia, George, and Wai Wan Tsang. "A simple method for generating |
| // gamma variables." ACM Transactions on Mathematical Software (TOMS) |
| // 26.3 (2000): 363-372. |
| d := a - 1.0/3 |
| c := 1 / (3 * math.Sqrt(d)) |
| for { |
| u := -exprnd() |
| x := normrnd() |
| v := 1 + x*c |
| v = v * v * v |
| if u < 0.5*x*x+d*(1-v+math.Log(v)) { |
| return d * v / b |
| } |
| } |
| } |
| panic("unreachable") |
| } |
| |
| // Survival returns the survival function (complementary CDF) at x. |
| func (g Gamma) Survival(x float64) float64 { |
| if x < 0 { |
| return 1 |
| } |
| return mathext.GammaIncComp(g.Alpha, g.Beta*x) |
| } |
| |
| // StdDev returns the standard deviation of the probability distribution. |
| func (g Gamma) StdDev() float64 { |
| return math.Sqrt(g.Variance()) |
| } |
| |
| // Variance returns the variance of the probability distribution. |
| func (g Gamma) Variance() float64 { |
| return g.Alpha / g.Beta / g.Beta |
| } |