From 855d38ebe85e404549f47a422bbe36a5bed31da7 Mon Sep 17 00:00:00 2001 From: gudvinr Date: Thu, 4 Nov 2021 12:08:21 +0300 Subject: [PATCH] stat/distuv: add logistic distribution --- stat/distuv/constants.go | 4 + stat/distuv/logistic.go | 97 ++++++++++++++++++++ stat/distuv/logistic_test.go | 172 +++++++++++++++++++++++++++++++++++ 3 files changed, 273 insertions(+) create mode 100644 stat/distuv/logistic.go create mode 100644 stat/distuv/logistic_test.go diff --git a/stat/distuv/constants.go b/stat/distuv/constants.go index 374d1993..3ebe6350 100644 --- a/stat/distuv/constants.go +++ b/stat/distuv/constants.go @@ -17,6 +17,10 @@ const ( // Euler–Mascheroni constant. eulerGamma = 0.5772156649015328606065120900824024310421593359399235988057672348848677267776646709369470632917467495146314472498070824809605 + + // sqrt3 is the value of sqrt(3) + // https://www.wolframalpha.com/input/?i=sqrt%283%29 + sqrt3 = 1.7320508075688772935274463415058723669428052538103806280558069794519330169088000370811461867572485756756261414154067030299699450 ) const ( diff --git a/stat/distuv/logistic.go b/stat/distuv/logistic.go new file mode 100644 index 00000000..50aa6f2a --- /dev/null +++ b/stat/distuv/logistic.go @@ -0,0 +1,97 @@ +// Copyright ©2021 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" +) + +// Logistic implements the Logistic distribution, a two-parameter distribution with support on the real axis. +// Its cumulative distribution function is the logistic function. +// +// General form of probability density fuction for Logistic distribution is +// E(x) / (s * (1 + E(x))^2) +// where E(x) = exp(-(x-μ)/s) +// +// For more information, see https://en.wikipedia.org/wiki/Logistic_distribution. +type Logistic struct { + Mu float64 // Mean value + S float64 // Scale parameter proportional to standard deviation +} + +// CDF computes the value of the cumulative density function at x. +func (l Logistic) CDF(x float64) float64 { + return 1 / (1 + math.Exp(-(x-l.Mu)/l.S)) +} + +// ExKurtosis returns the excess kurtosis of the distribution. +func (l Logistic) ExKurtosis() float64 { + return 6.0 / 5.0 +} + +// LogProb computes the natural logarithm of the value of the probability +// density function at x. +func (l Logistic) LogProb(x float64) float64 { + return x - 2*math.Log(math.Exp(x)+1) +} + +// Mean returns the mean of the probability distribution. +func (l Logistic) Mean() float64 { + return l.Mu +} + +// Mode returns the mode of the distribution. +// +// It is same as Mean for Logistic distribution. +func (l Logistic) Mode() float64 { + return l.Mu +} + +// Median returns the median of the distribution. +// +// It is same as Mean for Logistic distribution. +func (l Logistic) Median() float64 { + return l.Mu +} + +// NumParameters returns the number of parameters in the distribution. +// +// Always returns 2. +func (l Logistic) NumParameters() int { + return 2 +} + +// Prob computes the value of the probability density function at x. +func (l Logistic) Prob(x float64) float64 { + E := math.Exp(-(x - l.Mu) / l.S) + return E / (l.S * math.Pow(1+E, 2)) +} + +// Quantile returns the inverse of the cumulative distribution function. +func (l Logistic) Quantile(p float64) float64 { + return l.Mu + l.S*math.Log(p/(1-p)) +} + +// Skewness returns the skewness of the distribution. +// +// Always 0 for Logistic distribution. +func (l Logistic) Skewness() float64 { + return 0 +} + +// StdDev returns the standard deviation of the probability distribution. +func (l Logistic) StdDev() float64 { + return l.S * math.Pi / sqrt3 +} + +// Survival returns the survival function (complementary CDF) at x. +func (l Logistic) Survival(x float64) float64 { + return 1 - l.CDF(x) +} + +// Variance returns the variance of the probability distribution. +func (l Logistic) Variance() float64 { + return l.S * l.S * math.Pi * math.Pi / 3 +} diff --git a/stat/distuv/logistic_test.go b/stat/distuv/logistic_test.go new file mode 100644 index 00000000..30607149 --- /dev/null +++ b/stat/distuv/logistic_test.go @@ -0,0 +1,172 @@ +// Copyright ©2021 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" + "testing" + + "gonum.org/v1/gonum/floats/scalar" +) + +func TestLogisticParameters(t *testing.T) { + t.Parallel() + + var want float64 + + l := Logistic{Mu: 1, S: 0} + + want = 2 + if result := l.NumParameters(); result != int(want) { + t.Errorf("Wrong number of parameters: %d != %.0f", result, want) + } + + want = 6.0 / 5.0 + if result := l.ExKurtosis(); result != want { + t.Errorf("Wrong excess kurtosis: %f != %f", result, want) + } + + want = 0.0 + if result := l.Skewness(); result != want { + t.Errorf("Wrong skewness: %f != %f", result, want) + } + + want = l.Mu + if result := l.Mean(); result != want { + t.Errorf("Wrong mean value: %f != %f", result, want) + } + + want = l.Mu + if result := l.Median(); result != want { + t.Errorf("Wrong median value: %f != %f", result, want) + } + + want = l.Mu + if result := l.Mode(); result != want { + t.Errorf("Wrong mode value: %f != %f", result, want) + } +} + +func TestLogisticStdDev(t *testing.T) { + t.Parallel() + + l := Logistic{Mu: 0, S: sqrt3 / math.Pi} + + want := 1.0 + if result := l.StdDev(); !scalar.EqualWithinAbs(result, want, 1e-10) { + t.Errorf("Wrong StdDev with Mu=%f, S=%f: %f != %f", l.Mu, l.S, result, want) + } + + want = 1.0 + if result := l.Variance(); !scalar.EqualWithinAbs(result, want, 1e-10) { + t.Errorf("Wrong Variance with Mu=%f, S=%f: %f != %f", l.Mu, l.S, result, want) + } +} + +func TestLogisticCDF(t *testing.T) { + t.Parallel() + + // Values for "want" are taken from WolframAlpha: CDF[LogisticDistribution[mu,s], input] to 10 digits. + for _, v := range []struct { + mu, s, input, want float64 + }{ + {0.0, 0.0, 1.0, 1.0}, + {0.0, 1.0, 0.0, 0.5}, + {-0.5, 1.0, 0.0, 0.6224593312}, + {69.0, 420.0, 42.0, 0.4839341039}, + } { + l := Logistic{Mu: v.mu, S: v.s} + if result := l.CDF(v.input); !scalar.EqualWithinAbs(result, v.want, 1e-10) { + t.Errorf("Wrong CDF(%f) with Mu=%f, S=%f: %f != %f", v.input, l.Mu, l.S, result, v.want) + } + } + + // Edge case of zero in denominator. + l := Logistic{Mu: 0, S: 0} + + input := 0.0 + if result := l.CDF(input); !math.IsNaN(result) { + t.Errorf("Wrong CDF(%f) with Mu=%f, S=%f: %f is not NaN", input, l.Mu, l.S, result) + } +} + +// TestLogisticSurvival doesn't need excessive testing since it's just 1-CDF. +func TestLogisticSurvival(t *testing.T) { + t.Parallel() + + l := Logistic{Mu: 0, S: 1} + + input, want := 0.0, 0.5 + if result := l.Survival(input); result != want { + t.Errorf("Wrong Survival(%f) with Mu=%f, S=%f: %f != %f", input, l.Mu, l.S, result, want) + } +} + +func TestLogisticProb(t *testing.T) { + t.Parallel() + + // Values for "want" are taken from WolframAlpha: PDF[LogisticDistribution[mu,s], input] to 10 digits. + for _, v := range []struct { + mu, s, input, want float64 + }{ + {0.0, 1.0, 0.0, 0.25}, + {-0.5, 1.0, 0.0, 0.2350037122}, + {69.0, 420.0, 42.0, 0.0005946235404}, + } { + l := Logistic{Mu: v.mu, S: v.s} + if result := l.Prob(v.input); !scalar.EqualWithinAbs(result, v.want, 1e-10) { + t.Errorf("Wrong Prob(%f) with Mu=%f, S=%f: %.09f != %.09f", v.input, l.Mu, l.S, result, v.want) + } + } + + // Edge case of zero in denominator. + l := Logistic{Mu: 0, S: 0} + + input := 0.0 + if result := l.Prob(input); !math.IsNaN(result) { + t.Errorf("Wrong Prob(%f) with Mu=%f, S=%f: %f is not NaN", input, l.Mu, l.S, result) + } + + input = 1.0 + if result := l.Prob(input); !math.IsNaN(result) { + t.Errorf("Wrong Prob(%f) with Mu=%f, S=%f: %f is not NaN", input, l.Mu, l.S, result) + } +} + +func TestLogisticLogProb(t *testing.T) { + t.Parallel() + + l := Logistic{Mu: 0, S: 1} + + input, want := 0.0, -math.Log(4) + if result := l.LogProb(input); result != want { + t.Errorf("Wrong LogProb(%f) with Mu=%f, S=%f: %f != %f", input, l.Mu, l.S, result, want) + } +} + +func TestQuantile(t *testing.T) { + t.Parallel() + + for _, v := range []struct { + mu, s, input, want float64 + }{ + {0.0, 1.0, 0.5, 0.0}, + {0.0, 1.0, 0.0, math.Inf(-1)}, + {0.0, 1.0, 1.0, math.Inf(+1)}, + } { + l := Logistic{Mu: v.mu, S: v.s} + if result := l.Quantile(v.input); result != v.want { + t.Errorf("Wrong Quantile(%f) with Mu=%f, S=%f: %f != %f", v.input, l.Mu, l.S, result, v.want) + } + } + + // Edge case with NaN. + l := Logistic{Mu: 0, S: 0} + + input := 0.0 + if result := l.Quantile(input); !math.IsNaN(result) { + t.Errorf("Wrong Quantile(%f) with Mu=%f, S=%f: %f is not NaN", input, l.Mu, l.S, result) + } +}