stat/distuv: add logistic distribution

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 (

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
+}

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)
+	}
+}