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  }