github.com/jingcheng-WU/gonum@v0.9.1-0.20210323123734-f1a2a11a8f7b/stat/distuv/beta.go (about) 1 // Copyright ©2016 The Gonum Authors. All rights reserved. 2 // Use of this source code is governed by a BSD-style 3 // license that can be found in the LICENSE file. 4 5 package distuv 6 7 import ( 8 "math" 9 10 "golang.org/x/exp/rand" 11 12 "github.com/jingcheng-WU/gonum/mathext" 13 ) 14 15 // Beta implements the Beta distribution, a two-parameter continuous distribution 16 // with support between 0 and 1. 17 // 18 // The beta distribution has density function 19 // x^(α-1) * (1-x)^(β-1) * Γ(α+β) / (Γ(α)*Γ(β)) 20 // 21 // For more information, see https://en.wikipedia.org/wiki/Beta_distribution 22 type Beta struct { 23 // Alpha is the left shape parameter of the distribution. Alpha must be greater 24 // than 0. 25 Alpha float64 26 // Beta is the right shape parameter of the distribution. Beta must be greater 27 // than 0. 28 Beta float64 29 30 Src rand.Source 31 } 32 33 // CDF computes the value of the cumulative distribution function at x. 34 func (b Beta) CDF(x float64) float64 { 35 if x <= 0 { 36 return 0 37 } 38 if x >= 1 { 39 return 1 40 } 41 return mathext.RegIncBeta(b.Alpha, b.Beta, x) 42 } 43 44 // Entropy returns the differential entropy of the distribution. 45 func (b Beta) Entropy() float64 { 46 if b.Alpha <= 0 || b.Beta <= 0 { 47 panic("beta: negative parameters") 48 } 49 return mathext.Lbeta(b.Alpha, b.Beta) - (b.Alpha-1)*mathext.Digamma(b.Alpha) - 50 (b.Beta-1)*mathext.Digamma(b.Beta) + (b.Alpha+b.Beta-2)*mathext.Digamma(b.Alpha+b.Beta) 51 } 52 53 // ExKurtosis returns the excess kurtosis of the distribution. 54 func (b Beta) ExKurtosis() float64 { 55 num := 6 * ((b.Alpha-b.Beta)*(b.Alpha-b.Beta)*(b.Alpha+b.Beta+1) - b.Alpha*b.Beta*(b.Alpha+b.Beta+2)) 56 den := b.Alpha * b.Beta * (b.Alpha + b.Beta + 2) * (b.Alpha + b.Beta + 3) 57 return num / den 58 } 59 60 // LogProb computes the natural logarithm of the value of the probability 61 // density function at x. 62 func (b Beta) LogProb(x float64) float64 { 63 if x < 0 || x > 1 { 64 return math.Inf(-1) 65 } 66 67 if b.Alpha <= 0 || b.Beta <= 0 { 68 panic("beta: negative parameters") 69 } 70 71 lab, _ := math.Lgamma(b.Alpha + b.Beta) 72 la, _ := math.Lgamma(b.Alpha) 73 lb, _ := math.Lgamma(b.Beta) 74 var lx float64 75 if b.Alpha != 1 { 76 lx = (b.Alpha - 1) * math.Log(x) 77 } 78 var l1mx float64 79 if b.Beta != 1 { 80 l1mx = (b.Beta - 1) * math.Log(1-x) 81 } 82 return lab - la - lb + lx + l1mx 83 } 84 85 // Mean returns the mean of the probability distribution. 86 func (b Beta) Mean() float64 { 87 return b.Alpha / (b.Alpha + b.Beta) 88 } 89 90 // Mode returns the mode of the distribution. 91 // 92 // Mode returns NaN if both parametera are less than or equal to 1 as a special case, 93 // 0 if only Alpha <= 1 and 1 if only Beta <= 1. 94 func (b Beta) Mode() float64 { 95 if b.Alpha <= 1 { 96 if b.Beta <= 1 { 97 return math.NaN() 98 } 99 return 0 100 } 101 if b.Beta <= 1 { 102 return 1 103 } 104 return (b.Alpha - 1) / (b.Alpha + b.Beta - 2) 105 } 106 107 // NumParameters returns the number of parameters in the distribution. 108 func (b Beta) NumParameters() int { 109 return 2 110 } 111 112 // Prob computes the value of the probability density function at x. 113 func (b Beta) Prob(x float64) float64 { 114 return math.Exp(b.LogProb(x)) 115 } 116 117 // Quantile returns the inverse of the cumulative distribution function. 118 func (b Beta) Quantile(p float64) float64 { 119 if p < 0 || p > 1 { 120 panic(badPercentile) 121 } 122 return mathext.InvRegIncBeta(b.Alpha, b.Beta, p) 123 } 124 125 // Rand returns a random sample drawn from the distribution. 126 func (b Beta) Rand() float64 { 127 ga := Gamma{Alpha: b.Alpha, Beta: 1, Src: b.Src}.Rand() 128 gb := Gamma{Alpha: b.Beta, Beta: 1, Src: b.Src}.Rand() 129 return ga / (ga + gb) 130 } 131 132 // StdDev returns the standard deviation of the probability distribution. 133 func (b Beta) StdDev() float64 { 134 return math.Sqrt(b.Variance()) 135 } 136 137 // Survival returns the survival function (complementary CDF) at x. 138 func (b Beta) Survival(x float64) float64 { 139 switch { 140 case x <= 0: 141 return 1 142 case x >= 1: 143 return 0 144 } 145 return mathext.RegIncBeta(b.Beta, b.Alpha, 1-x) 146 } 147 148 // Variance returns the variance of the probability distribution. 149 func (b Beta) Variance() float64 { 150 return b.Alpha * b.Beta / ((b.Alpha + b.Beta) * (b.Alpha + b.Beta) * (b.Alpha + b.Beta + 1)) 151 }