Skip to content
Permalink
Branch: master
Find file Copy path
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
261 lines (227 sloc) 7.31 KB
// Copyright ©2014 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/floats"
"gonum.org/v1/gonum/stat"
)
// Exponential represents the exponential distribution (https://en.wikipedia.org/wiki/Exponential_distribution).
type Exponential struct {
Rate float64
Src rand.Source
}
// CDF computes the value of the cumulative density function at x.
func (e Exponential) CDF(x float64) float64 {
if x < 0 {
return 0
}
return 1 - math.Exp(-e.Rate*x)
}
// ConjugateUpdate updates the parameters of the distribution from the sufficient
// statistics of a set of samples. The sufficient statistics, suffStat, have been
// observed with nSamples observations. The prior values of the distribution are those
// currently in the distribution, and have been observed with priorStrength samples.
//
// For the exponential distribution, the sufficient statistic is the inverse of
// the mean of the samples.
// The prior is having seen priorStrength[0] samples with inverse mean Exponential.Rate
// As a result of this function, Exponential.Rate is updated based on the weighted
// samples, and priorStrength is modified to include the new number of samples observed.
//
// This function panics if len(suffStat) != 1 or len(priorStrength) != 1.
func (e *Exponential) ConjugateUpdate(suffStat []float64, nSamples float64, priorStrength []float64) {
if len(suffStat) != 1 {
panic("exponential: incorrect suffStat length")
}
if len(priorStrength) != 1 {
panic("exponential: incorrect priorStrength length")
}
totalSamples := nSamples + priorStrength[0]
totalSum := nSamples / suffStat[0]
if !(priorStrength[0] == 0) {
totalSum += priorStrength[0] / e.Rate
}
e.Rate = totalSamples / totalSum
priorStrength[0] = totalSamples
}
// Entropy returns the entropy of the distribution.
func (e Exponential) Entropy() float64 {
return 1 - math.Log(e.Rate)
}
// ExKurtosis returns the excess kurtosis of the distribution.
func (Exponential) ExKurtosis() float64 {
return 6
}
// Fit sets the parameters of the probability distribution from the
// data samples x with relative weights w.
// If weights is nil, then all the weights are 1.
// If weights is not nil, then the len(weights) must equal len(samples).
func (e *Exponential) Fit(samples, weights []float64) {
suffStat := make([]float64, e.NumSuffStat())
nSamples := e.SuffStat(suffStat, samples, weights)
e.ConjugateUpdate(suffStat, nSamples, make([]float64, e.NumSuffStat()))
}
// LogProb computes the natural logarithm of the value of the probability density function at x.
func (e Exponential) LogProb(x float64) float64 {
if x < 0 {
return math.Inf(-1)
}
return math.Log(e.Rate) - e.Rate*x
}
// Mean returns the mean of the probability distribution.
func (e Exponential) Mean() float64 {
return 1 / e.Rate
}
// Median returns the median of the probability distribution.
func (e Exponential) Median() float64 {
return math.Ln2 / e.Rate
}
// Mode returns the mode of the probability distribution.
func (Exponential) Mode() float64 {
return 0
}
// NumParameters returns the number of parameters in the distribution.
func (Exponential) NumParameters() int {
return 1
}
// NumSuffStat returns the number of sufficient statistics for the distribution.
func (Exponential) NumSuffStat() int {
return 1
}
// Prob computes the value of the probability density function at x.
func (e Exponential) Prob(x float64) float64 {
return math.Exp(e.LogProb(x))
}
// Quantile returns the inverse of the cumulative probability distribution.
func (e Exponential) Quantile(p float64) float64 {
if p < 0 || p > 1 {
panic(badPercentile)
}
return -math.Log(1-p) / e.Rate
}
// Rand returns a random sample drawn from the distribution.
func (e Exponential) Rand() float64 {
var rnd float64
if e.Src == nil {
rnd = rand.ExpFloat64()
} else {
rnd = rand.New(e.Src).ExpFloat64()
}
return rnd / e.Rate
}
// Score returns the score function with respect to the parameters of the
// distribution at the input location x. The score function is the derivative
// of the log-likelihood at x with respect to the parameters
// (∂/∂θ) log(p(x;θ))
// If deriv is non-nil, len(deriv) must equal the number of parameters otherwise
// Score will panic, and the derivative is stored in-place into deriv. If deriv
// is nil a new slice will be allocated and returned.
//
// The order is [∂LogProb / ∂Rate].
//
// For more information, see https://en.wikipedia.org/wiki/Score_%28statistics%29.
//
// Special cases:
// Score(0) = [NaN]
func (e Exponential) Score(deriv []float64, x float64) []float64 {
if deriv == nil {
deriv = make([]float64, e.NumParameters())
}
if len(deriv) != e.NumParameters() {
panic(badLength)
}
if x > 0 {
deriv[0] = 1/e.Rate - x
return deriv
}
if x < 0 {
deriv[0] = 0
return deriv
}
deriv[0] = math.NaN()
return deriv
}
// ScoreInput returns the score function with respect to the input of the
// distribution at the input location specified by x. The score function is the
// derivative of the log-likelihood
// (d/dx) log(p(x)) .
// Special cases:
// ScoreInput(0) = NaN
func (e Exponential) ScoreInput(x float64) float64 {
if x > 0 {
return -e.Rate
}
if x < 0 {
return 0
}
return math.NaN()
}
// Skewness returns the skewness of the distribution.
func (Exponential) Skewness() float64 {
return 2
}
// StdDev returns the standard deviation of the probability distribution.
func (e Exponential) StdDev() float64 {
return 1 / e.Rate
}
// SuffStat computes the sufficient statistics of set of samples to update
// the distribution. The sufficient statistics are stored in place, and the
// effective number of samples are returned.
//
// The exponential distribution has one sufficient statistic, the average rate
// of the samples.
//
// If weights is nil, the weights are assumed to be 1, otherwise panics if
// len(samples) != len(weights). Panics if len(suffStat) != NumSuffStat().
func (Exponential) SuffStat(suffStat, samples, weights []float64) (nSamples float64) {
if len(weights) != 0 && len(samples) != len(weights) {
panic(badLength)
}
if len(suffStat) != (Exponential{}).NumSuffStat() {
panic(badSuffStat)
}
if len(weights) == 0 {
nSamples = float64(len(samples))
} else {
nSamples = floats.Sum(weights)
}
mean := stat.Mean(samples, weights)
suffStat[0] = 1 / mean
return nSamples
}
// Survival returns the survival function (complementary CDF) at x.
func (e Exponential) Survival(x float64) float64 {
if x < 0 {
return 1
}
return math.Exp(-e.Rate * x)
}
// setParameters modifies the parameters of the distribution.
func (e *Exponential) setParameters(p []Parameter) {
if len(p) != e.NumParameters() {
panic("exponential: incorrect number of parameters to set")
}
if p[0].Name != "Rate" {
panic("exponential: " + panicNameMismatch)
}
e.Rate = p[0].Value
}
// Variance returns the variance of the probability distribution.
func (e Exponential) Variance() float64 {
return 1 / (e.Rate * e.Rate)
}
// parameters returns the parameters of the distribution.
func (e Exponential) parameters(p []Parameter) []Parameter {
nParam := e.NumParameters()
if p == nil {
p = make([]Parameter, nParam)
} else if len(p) != nParam {
panic("exponential: improper parameter length")
}
p[0].Name = "Rate"
p[0].Value = e.Rate
return p
}
You can’t perform that action at this time.