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
|
// 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"
)
// Beta implements the Beta distribution, a two-parameter continuous distribution
// with support between 0 and 1.
//
// The beta distribution has density function
//
// x^(α-1) * (1-x)^(β-1) * Γ(α+β) / (Γ(α)*Γ(β))
//
// For more information, see https://en.wikipedia.org/wiki/Beta_distribution
type Beta struct {
// Alpha is the left shape parameter of the distribution. Alpha must be greater
// than 0.
Alpha float64
// Beta is the right shape parameter of the distribution. Beta must be greater
// than 0.
Beta float64
Src rand.Source
}
// CDF computes the value of the cumulative distribution function at x.
func (b Beta) CDF(x float64) float64 {
if x <= 0 {
return 0
}
if x >= 1 {
return 1
}
return mathext.RegIncBeta(b.Alpha, b.Beta, x)
}
// Entropy returns the differential entropy of the distribution.
func (b Beta) Entropy() float64 {
if b.Alpha <= 0 || b.Beta <= 0 {
panic("beta: negative parameters")
}
return mathext.Lbeta(b.Alpha, b.Beta) - (b.Alpha-1)*mathext.Digamma(b.Alpha) -
(b.Beta-1)*mathext.Digamma(b.Beta) + (b.Alpha+b.Beta-2)*mathext.Digamma(b.Alpha+b.Beta)
}
// ExKurtosis returns the excess kurtosis of the distribution.
func (b Beta) ExKurtosis() float64 {
num := 6 * ((b.Alpha-b.Beta)*(b.Alpha-b.Beta)*(b.Alpha+b.Beta+1) - b.Alpha*b.Beta*(b.Alpha+b.Beta+2))
den := b.Alpha * b.Beta * (b.Alpha + b.Beta + 2) * (b.Alpha + b.Beta + 3)
return num / den
}
// LogProb computes the natural logarithm of the value of the probability
// density function at x.
func (b Beta) LogProb(x float64) float64 {
if x < 0 || x > 1 {
return math.Inf(-1)
}
if b.Alpha <= 0 || b.Beta <= 0 {
panic("beta: negative parameters")
}
lab, _ := math.Lgamma(b.Alpha + b.Beta)
la, _ := math.Lgamma(b.Alpha)
lb, _ := math.Lgamma(b.Beta)
var lx float64
if b.Alpha != 1 {
lx = (b.Alpha - 1) * math.Log(x)
}
var l1mx float64
if b.Beta != 1 {
l1mx = (b.Beta - 1) * math.Log(1-x)
}
return lab - la - lb + lx + l1mx
}
// Mean returns the mean of the probability distribution.
func (b Beta) Mean() float64 {
return b.Alpha / (b.Alpha + b.Beta)
}
// Mode returns the mode of the distribution.
//
// Mode returns NaN if both parameters are less than or equal to 1 as a special case,
// 0 if only Alpha <= 1 and 1 if only Beta <= 1.
func (b Beta) Mode() float64 {
if b.Alpha <= 1 {
if b.Beta <= 1 {
return math.NaN()
}
return 0
}
if b.Beta <= 1 {
return 1
}
return (b.Alpha - 1) / (b.Alpha + b.Beta - 2)
}
// NumParameters returns the number of parameters in the distribution.
func (b Beta) NumParameters() int {
return 2
}
// Prob computes the value of the probability density function at x.
func (b Beta) Prob(x float64) float64 {
return math.Exp(b.LogProb(x))
}
// Quantile returns the inverse of the cumulative distribution function.
func (b Beta) Quantile(p float64) float64 {
if p < 0 || p > 1 {
panic(badPercentile)
}
return mathext.InvRegIncBeta(b.Alpha, b.Beta, p)
}
// Rand returns a random sample drawn from the distribution.
func (b Beta) Rand() float64 {
ga := Gamma{Alpha: b.Alpha, Beta: 1, Src: b.Src}.Rand()
gb := Gamma{Alpha: b.Beta, Beta: 1, Src: b.Src}.Rand()
return ga / (ga + gb)
}
// StdDev returns the standard deviation of the probability distribution.
func (b Beta) StdDev() float64 {
return math.Sqrt(b.Variance())
}
// Survival returns the survival function (complementary CDF) at x.
func (b Beta) Survival(x float64) float64 {
switch {
case x <= 0:
return 1
case x >= 1:
return 0
}
return mathext.RegIncBeta(b.Beta, b.Alpha, 1-x)
}
// Variance returns the variance of the probability distribution.
func (b Beta) Variance() float64 {
return b.Alpha * b.Beta / ((b.Alpha + b.Beta) * (b.Alpha + b.Beta) * (b.Alpha + b.Beta + 1))
}
|
