stat/distuv: add chi distribution

stat/distuv/chi.go

@@ -0,0 +1,124 @@
+// 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"
+
+	"golang.org/x/exp/rand"
+
+	"gonum.org/v1/gonum/mathext"
+)
+
+// Chi implements the χ distribution, a one parameter distribution
+// with support on the positive numbers.
+//
+// The density function is given by
+//  1/(2^{k/2-1} * Γ(k/2)) * x^{k - 1} * e^{-x^2/2}
+//
+// For more information, see https://en.wikipedia.org/wiki/Chi_distribution.
+type Chi struct {
+	// K is the shape parameter, corresponding to the degrees of freedom. Must
+	// be greater than 0.
+	K float64
+
+	Src rand.Source
+}
+
+// CDF computes the value of the cumulative density function at x.
+func (c Chi) CDF(x float64) float64 {
+	return mathext.GammaIncReg(c.K/2, (x*x)/2)
+}
+
+// Entropy returns the differential entropy of the distribution.
+func (c Chi) Entropy() float64 {
+	lg, _ := math.Lgamma(c.K / 2)
+	return lg + 0.5*(c.K-math.Ln2-(c.K-1)*mathext.Digamma(c.K/2))
+}
+
+// ExKurtosis returns the excess kurtosis of the distribution.
+func (c Chi) ExKurtosis() float64 {
+	v := c.Variance()
+	s := math.Sqrt(v)
+	return 2 / v * (1 - c.Mean()*s*c.Skewness() - v)
+}
+
+// LogProb computes the natural logarithm of the value of the probability
+// density function at x.
+func (c Chi) LogProb(x float64) float64 {
+	if x < 0 {
+		return math.Inf(-1)
+	}
+	lg, _ := math.Lgamma(c.K / 2)
+	return (c.K-1)*math.Log(x) - (x*x)/2 - (c.K/2-1)*math.Ln2 - lg
+}
+
+// Mean returns the mean of the probability distribution.
+func (c Chi) Mean() float64 {
+	lg1, _ := math.Lgamma((c.K + 1) / 2)
+	lg, _ := math.Lgamma(c.K / 2)
+	return math.Sqrt2 * math.Exp(lg1-lg)
+}
+
+// Median returns the median of the distribution.
+func (c Chi) Median() float64 {
+	return c.Quantile(0.5)
+}
+
+// Mode returns the mode of the distribution.
+//
+// Mode returns NaN if K is less than one.
+func (c Chi) Mode() float64 {
+	return math.Sqrt(c.K - 1)
+}
+
+// NumParameters returns the number of parameters in the distribution.
+func (c Chi) NumParameters() int {
+	return 1
+}
+
+// Prob computes the value of the probability density function at x.
+func (c Chi) Prob(x float64) float64 {
+	return math.Exp(c.LogProb(x))
+}
+
+// Rand returns a random sample drawn from the distribution.
+func (c Chi) Rand() float64 {
+	return math.Sqrt(Gamma{c.K / 2, 0.5, c.Src}.Rand())
+}
+
+// Quantile returns the inverse of the cumulative distribution function.
+func (c Chi) Quantile(p float64) float64 {
+	if p < 0 || 1 < p {
+		panic(badPercentile)
+	}
+	return math.Sqrt(2 * mathext.GammaIncRegInv(0.5*c.K, p))
+}
+
+// Skewness returns the skewness of the distribution.
+func (c Chi) Skewness() float64 {
+	v := c.Variance()
+	s := math.Sqrt(v)
+	return c.Mean() / (s * v) * (1 - 2*v)
+}
+
+// StdDev returns the standard deviation of the probability distribution.
+func (c Chi) StdDev() float64 {
+	return math.Sqrt(c.Variance())
+}
+
+// Survival returns the survival function (complementary CDF) at x.
+func (c Chi) Survival(x float64) float64 {
+	if x < 0 {
+		return 1
+	}
+	return mathext.GammaIncRegComp(0.5*c.K, 0.5*(x*x))
+}
+
+// Variance returns the variance of the probability distribution.
+func (c Chi) Variance() float64 {
+	m := c.Mean()
+	return math.Max(0, c.K-m*m)
+}

stat/distuv/chi_test.go

@@ -0,0 +1,99 @@
+// 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"
+	"sort"
+	"testing"
+
+	"golang.org/x/exp/rand"
+
+	"gonum.org/v1/gonum/floats/scalar"
+)
+
+func TestChiProb(t *testing.T) {
+	t.Parallel()
+	for _, test := range []struct {
+		x, k, want float64
+	}{
+		{10, 3, 1.538919725341288e-20},
+		{2.3, 3, 0.2997000593061405},
+		{0.8, 0.2, 0.1702707693447167},
+	} {
+		pdf := Chi{test.k, nil}.Prob(test.x)
+		if !scalar.EqualWithinAbsOrRel(pdf, test.want, 1e-10, 1e-10) {
+			t.Errorf("Pdf mismatch, x = %v, K = %v. Got %v, want %v", test.x, test.k, pdf, test.want)
+		}
+	}
+}
+
+func TestChiCDF(t *testing.T) {
+	t.Parallel()
+	for _, test := range []struct {
+		x, k, want float64
+	}{
+		// Values calculated with scipy.stats.chi.cdf
+		{0, 1, 0},
+		{0.01, 5, 5.319040436531812e-12},
+		{0.05, 3, 3.3220267268523235e-05},
+		{0.5, 2, 0.1175030974154046},
+		{0.95, 3, 0.17517554009157732},
+		{0.99, 5, 0.035845177452671864},
+		{1, 1, 0.6826894921370859},
+		{1.5, 4, 0.3101135068635068},
+		{10, 10, 1},
+		{25, 15, 1},
+	} {
+		cdf := Chi{test.k, nil}.CDF(test.x)
+		if !scalar.EqualWithinAbsOrRel(cdf, test.want, 1e-10, 1e-10) {
+			t.Errorf("CDF mismatch, x = %v, K = %v. Got %v, want %v", test.x, test.k, cdf, test.want)
+		}
+	}
+}
+
+func TestChi(t *testing.T) {
+	t.Parallel()
+	src := rand.New(rand.NewSource(1))
+	for i, b := range []Chi{
+		{3, src},
+		{1.5, src},
+		{0.9, src},
+	} {
+		testChi(t, b, i)
+	}
+}
+
+func testChi(t *testing.T, c Chi, i int) {
+	const (
+		tol  = 1e-2
+		n    = 1e6
+		bins = 50
+	)
+	x := make([]float64, n)
+	generateSamples(x, c)
+	sort.Float64s(x)
+
+	testRandLogProbContinuous(t, i, 0, x, c, tol, bins)
+	checkMean(t, i, x, c, tol)
+	checkMedian(t, i, x, c, tol)
+	checkVarAndStd(t, i, x, c, tol)
+	checkEntropy(t, i, x, c, tol)
+	checkExKurtosis(t, i, x, c, 7e-2)
+	checkProbContinuous(t, i, x, 0, math.Inf(1), c, 1e-5)
+	checkQuantileCDFSurvival(t, i, x, c, 1e-2)
+
+	expectedMode := math.Sqrt(c.K - 1)
+	if !scalar.Same(c.Mode(), expectedMode) {
+		t.Errorf("Mode is not equal to sqrt(k - 1). Got %v, want %v", c.Mode(), expectedMode)
+	}
+	if c.NumParameters() != 1 {
+		t.Errorf("NumParameters is not 1. Got %v", c.NumParameters())
+	}
+	survival := c.Survival(-0.00001)
+	if survival != 1 {
+		t.Errorf("Survival is not 1 for negative argument. Got %v", survival)
+	}
+}