stat/distuv/laplace.go

// 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"
	"math/rand"
	"sort"

	"gonum.org/v1/gonum/floats"
	"gonum.org/v1/gonum/stat"
)

// Laplace represents the Laplace distribution (https://en.wikipedia.org/wiki/Laplace_distribution).
type Laplace struct {
	Mu     float64 // Mean of the Laplace distribution
	Scale  float64 // Scale of the Laplace distribution
	Source *rand.Rand
}

// CDF computes the value of the cumulative density function at x.
func (l Laplace) CDF(x float64) float64 {
	if x < l.Mu {
		return 0.5 * math.Exp((x-l.Mu)/l.Scale)
	}
	return 1 - 0.5*math.Exp(-(x-l.Mu)/l.Scale)
}

// Entropy returns the entropy of the distribution.
func (l Laplace) Entropy() float64 {
	return 1 + math.Log(2*l.Scale)
}

// ExKurtosis returns the excess kurtosis of the distribution.
func (l Laplace) ExKurtosis() float64 {
	return 3
}

// 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).
//
// Note: Laplace distribution has no FitPrior because it has no sufficient
// statistics.
func (l *Laplace) Fit(samples, weights []float64) {
	if len(samples) != len(weights) {
		panic(badLength)
	}

	if len(samples) == 0 {
		panic(badNoSamples)
	}
	if len(samples) == 1 {
		l.Mu = samples[0]
		l.Scale = 0
		return
	}

	var (
		sortedSamples []float64
		sortedWeights []float64
	)
	if sort.Float64sAreSorted(samples) {
		sortedSamples = samples
		sortedWeights = weights
	} else {
		// Need to copy variables so the input variables aren't effected by the sorting
		sortedSamples = make([]float64, len(samples))
		copy(sortedSamples, samples)
		sortedWeights := make([]float64, len(samples))
		copy(sortedWeights, weights)

		stat.SortWeighted(sortedSamples, sortedWeights)
	}

	// The (weighted) median of the samples is the maximum likelihood estimate
	// of the mean parameter
	// TODO: Rethink quantile type when stat has more options
	l.Mu = stat.Quantile(0.5, stat.Empirical, sortedSamples, sortedWeights)

	sumWeights := floats.Sum(weights)

	// The scale parameter is the average absolute distance
	// between the sample and the mean
	absError := stat.MomentAbout(1, samples, l.Mu, weights)

	l.Scale = absError / sumWeights
}

// LogProb computes the natural logarithm of the value of the probability density
// function at x.
func (l Laplace) LogProb(x float64) float64 {
	return -math.Ln2 - math.Log(l.Scale) - math.Abs(x-l.Mu)/l.Scale
}

// MarshalParameters implements the ParameterMarshaler interface
func (l Laplace) MarshalParameters(p []Parameter) {
	if len(p) != l.NumParameters() {
		panic(badLength)
	}
	p[0].Name = "Mu"
	p[0].Value = l.Mu
	p[1].Name = "Scale"
	p[1].Value = l.Scale
}

// Mean returns the mean of the probability distribution.
func (l Laplace) Mean() float64 {
	return l.Mu
}

// Median returns the median of the LaPlace distribution.
func (l Laplace) Median() float64 {
	return l.Mu
}

// Mode returns the mode of the LaPlace distribution.
func (l Laplace) Mode() float64 {
	return l.Mu
}

// NumParameters returns the number of parameters in the distribution.
func (l Laplace) NumParameters() int {
	return 2
}

// Quantile returns the inverse of the cumulative probability distribution.
func (l Laplace) Quantile(p float64) float64 {
	if p < 0 || p > 1 {
		panic(badPercentile)
	}
	if p < 0.5 {
		return l.Mu + l.Scale*math.Log(1+2*(p-0.5))
	}
	return l.Mu - l.Scale*math.Log(1-2*(p-0.5))
}

// Prob computes the value of the probability density function at x.
func (l Laplace) Prob(x float64) float64 {
	return math.Exp(l.LogProb(x))
}

// Rand returns a random sample drawn from the distribution.
func (l Laplace) Rand() float64 {
	var rnd float64
	if l.Source == nil {
		rnd = rand.Float64()
	} else {
		rnd = l.Source.Float64()
	}
	u := rnd - 0.5
	if u < 0 {
		return l.Mu + l.Scale*math.Log(1+2*u)
	}
	return l.Mu - l.Scale*math.Log(1-2*u)
}

// 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 / ∂Mu, ∂LogProb / ∂Scale].
//
// For more information, see https://en.wikipedia.org/wiki/Score_%28statistics%29.
//
// Special cases:
//  Score(0) = [0, -0.5/l.Scale]
func (l Laplace) Score(deriv []float64, x float64) []float64 {
	if deriv == nil {
		deriv = make([]float64, l.NumParameters())
	}
	if len(deriv) != l.NumParameters() {
		panic(badLength)
	}
	diff := x - l.Mu
	if diff > 0 {
		deriv[0] = 1 / l.Scale
	} else if diff < 0 {
		deriv[0] = -1 / l.Scale
	} else if diff == 0 {
		deriv[0] = 0
	} else {
		// must be NaN
		deriv[0] = math.NaN()
	}

	deriv[1] = math.Abs(diff)/(l.Scale*l.Scale) - 0.5/(l.Scale)
	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(l.Mu) = 0
func (l Laplace) ScoreInput(x float64) float64 {
	diff := x - l.Mu
	if diff == 0 {
		return 0
	}
	if diff > 0 {
		return -1 / l.Scale
	}
	return 1 / l.Scale
}

// Skewness returns the skewness of the distribution.
func (Laplace) Skewness() float64 {
	return 0
}

// StdDev returns the standard deviation of the distribution.
func (l Laplace) StdDev() float64 {
	return math.Sqrt2 * l.Scale
}

// Survival returns the survival function (complementary CDF) at x.
func (l Laplace) Survival(x float64) float64 {
	if x < l.Mu {
		return 1 - 0.5*math.Exp((x-l.Mu)/l.Scale)
	}
	return 0.5 * math.Exp(-(x-l.Mu)/l.Scale)
}

// UnmarshalParameters implements the ParameterMarshaler interface
func (l *Laplace) UnmarshalParameters(p []Parameter) {
	if len(p) != l.NumParameters() {
		panic(badLength)
	}
	if p[0].Name != "Mu" {
		panic("laplace: " + panicNameMismatch)
	}
	if p[1].Name != "Scale" {
		panic("laplace: " + panicNameMismatch)
	}
	l.Mu = p[0].Value
	l.Scale = p[1].Value
}

// Variance returns the variance of the probability distribution.
func (l Laplace) Variance() float64 {
	return 2 * l.Scale * l.Scale
}