diff --git a/stat/distuv/chi.go b/stat/distuv/chi.go new file mode 100644 index 00000000..82099d58 --- /dev/null +++ b/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) +} diff --git a/stat/distuv/chi_test.go b/stat/distuv/chi_test.go new file mode 100644 index 00000000..51b7a066 --- /dev/null +++ b/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) + } +}