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
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
|
// 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"
"golang.org/x/exp/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
Src rand.Source
}
// 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.GammaIncReg(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.GammaIncRegInv(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 {
const (
// The 0.2 threshold is from https://www4.stat.ncsu.edu/~rmartin/Codes/rgamss.R
// described in detail in https://arxiv.org/abs/1302.1884.
smallAlphaThresh = 0.2
)
if g.Beta <= 0 {
panic("gamma: beta <= 0")
}
unifrnd := rand.Float64
exprnd := rand.ExpFloat64
normrnd := rand.NormFloat64
if g.Src != nil {
rnd := rand.New(g.Src)
unifrnd = rnd.Float64
exprnd = rnd.ExpFloat64
normrnd = rnd.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 < smallAlphaThresh:
// 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
lr := -math.Log1p(1 / lambda / math.E)
for {
e := exprnd()
var z float64
if e >= -lr {
z = e + lr
} else {
z = -exprnd() / lambda
}
eza := math.Exp(-z / a)
lh := -z - eza
var lEta float64
if z >= 0 {
lEta = -z
} else {
lEta = -1 + lambda*z
}
if lh-lEta > -exprnd() {
return eza / b
}
}
case a >= smallAlphaThresh:
// 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
m := 1.0
if a < 1 {
d += 1.0
m = math.Pow(unifrnd(), 1/a)
}
c := 1 / (3 * math.Sqrt(d))
for {
x := normrnd()
v := 1 + x*c
if v <= 0.0 {
continue
}
v = v * v * v
u := unifrnd()
if u < 1.0-0.0331*(x*x)*(x*x) {
return m * d * v / b
}
if math.Log(u) < 0.5*x*x+d*(1-v+math.Log(v)) {
return m * 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.GammaIncRegComp(g.Alpha, g.Beta*x)
}
// StdDev returns the standard deviation of the probability distribution.
func (g Gamma) StdDev() float64 {
return math.Sqrt(g.Alpha) / g.Beta
}
// Variance returns the variance of the probability distribution.
func (g Gamma) Variance() float64 {
return g.Alpha / g.Beta / g.Beta
}
|
