Repository: gonum/gonum
Pull request #1660 (2021-06-10 23:55:29)

stat/distuv: add chi distribution
Please take a look.

<!--
Checklist:

- API changes have been discussed
- code is goformated correctly (goimports)
- packages with generated code have had code generation run
- tests pass locally
- linked to relevant issues

Please make sure your commit message summary line and pull request
title match the Go convention; a one-line summary of the change,
prefixed by the primary affected package that should complete the
sentence, "This change modifies Gonum to _____."
-->


kortschak: @romanwerpachowski This fails for reasons that I don't understand. The `Prob`, `Quantile` and `Rand` functions seem good, but very clearly `ExKurtosis`, `Mean` and `Variance` lie.
codecov-commenter: # [Codecov](https://codecov.io/gh/gonum/gonum/pull/1660?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum) Report
> Merging [#1660](https://codecov.io/gh/gonum/gonum/pull/1660?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum) (4ae061e) into [master](https://codecov.io/gh/gonum/gonum/commit/97f387b38f8e1baed300c7cc3d4b2dd9938820da?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum) (97f387b) will **decrease** coverage by `0.14%`.
> The diff coverage is `92.59%`.

> :exclamation: Current head 4ae061e differs from pull request most recent head a3b9dd0. Consider uploading reports for the commit a3b9dd0 to get more accurate results
[![Impacted file tree graph](https://codecov.io/gh/gonum/gonum/pull/1660/graphs/tree.svg?width=650&height=150&src=pr&token=B7qHnqSPNL&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum)](https://codecov.io/gh/gonum/gonum/pull/1660?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum)

```diff
@@            Coverage Diff             @@
##           master    #1660      +/-   ##
==========================================
- Coverage   73.22%   73.07%   -0.15%     
==========================================
  Files         504      489      -15     
  Lines       59811    59305     -506     
==========================================
- Hits        43794    43339     -455     
+ Misses      13537    13512      -25     
+ Partials     2480     2454      -26     
```


| [Impacted Files](https://codecov.io/gh/gonum/gonum/pull/1660?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum) | Coverage Δ | |
|---|---|---|
| [stat/distuv/chi.go](https://codecov.io/gh/gonum/gonum/pull/1660/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-c3RhdC9kaXN0dXYvY2hpLmdv) | `92.59% <92.59%> (ø)` | |
| [internal/asm/c64/conj.go](https://codecov.io/gh/gonum/gonum/pull/1660/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-aW50ZXJuYWwvYXNtL2M2NC9jb25qLmdv) | `0.00% <0.00%> (-100.00%)` | :arrow_down: |
| [internal/asm/c128/stubs\_noasm.go](https://codecov.io/gh/gonum/gonum/pull/1660/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-aW50ZXJuYWwvYXNtL2MxMjgvc3R1YnNfbm9hc20uZ28=) | | |
| [internal/asm/f64/l2norm\_noasm.go](https://codecov.io/gh/gonum/gonum/pull/1660/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-aW50ZXJuYWwvYXNtL2Y2NC9sMm5vcm1fbm9hc20uZ28=) | | |
| [internal/asm/f64/stubs\_noasm.go](https://codecov.io/gh/gonum/gonum/pull/1660/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-aW50ZXJuYWwvYXNtL2Y2NC9zdHVic19ub2FzbS5nbw==) | | |
| [internal/math32/sqrt.go](https://codecov.io/gh/gonum/gonum/pull/1660/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-aW50ZXJuYWwvbWF0aDMyL3NxcnQuZ28=) | | |
| [mat/offset\_appengine.go](https://codecov.io/gh/gonum/gonum/pull/1660/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-bWF0L29mZnNldF9hcHBlbmdpbmUuZ28=) | | |
| [internal/asm/f32/stubs\_noasm.go](https://codecov.io/gh/gonum/gonum/pull/1660/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-aW50ZXJuYWwvYXNtL2YzMi9zdHVic19ub2FzbS5nbw==) | | |
| [mat/index\_bound\_checks.go](https://codecov.io/gh/gonum/gonum/pull/1660/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-bWF0L2luZGV4X2JvdW5kX2NoZWNrcy5nbw==) | | |
| [internal/asm/f64/axpy.go](https://codecov.io/gh/gonum/gonum/pull/1660/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-aW50ZXJuYWwvYXNtL2Y2NC9heHB5Lmdv) | | |
| ... and [10 more](https://codecov.io/gh/gonum/gonum/pull/1660/diff?src=pr&el=tree-more&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum) | |

------

[Continue to review full report at Codecov](https://codecov.io/gh/gonum/gonum/pull/1660?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/gonum/gonum/pull/1660?src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum). Last update [97f387b...a3b9dd0](https://codecov.io/gh/gonum/gonum/pull/1660?src=pr&el=lastupdated&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum).

kortschak: PTAL

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


Test patch:
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)
+	}
+}