Repository: gonum/gonum
Pull request #1653 (2021-06-05 17:29:18)

lapack/gonum: prevent division of untyped ints always yielding zero in Dlasy2
Please take a look at the line I changed. The following operation:
```go
scale = 1 / 8 / maxbtmp
```
always yields zero. I added decimal points to `1` and `8` to fix this undesired behavior.



codecov-commenter: # [Codecov](https://codecov.io/gh/gonum/gonum/pull/1653?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum) Report
> Merging [#1653](https://codecov.io/gh/gonum/gonum/pull/1653?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum) (ff49396) into [master](https://codecov.io/gh/gonum/gonum/commit/2b838f4d6e91ac5a35b0961f5386468df1b44014?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum) (2b838f4) will **decrease** coverage by `0.15%`.
> The diff coverage is `0.00%`.

[![Impacted file tree graph](https://codecov.io/gh/gonum/gonum/pull/1653/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/1653?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum)

```diff
@@            Coverage Diff             @@
##           master    #1653      +/-   ##
==========================================
- Coverage   73.22%   73.06%   -0.16%     
==========================================
  Files         504      488      -16     
  Lines       59811    59274     -537     
==========================================
- Hits        43794    43310     -484     
+ Misses      13537    13511      -26     
+ Partials     2480     2453      -27     
```


| [Impacted Files](https://codecov.io/gh/gonum/gonum/pull/1653?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum) | Coverage Δ | |
|---|---|---|
| [lapack/gonum/dlasy2.go](https://codecov.io/gh/gonum/gonum/pull/1653/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-bGFwYWNrL2dvbnVtL2RsYXN5Mi5nbw==) | `88.00% <0.00%> (ø)` | |
| [internal/asm/c64/conj.go](https://codecov.io/gh/gonum/gonum/pull/1653/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/f64/ge\_noasm.go](https://codecov.io/gh/gonum/gonum/pull/1653/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-aW50ZXJuYWwvYXNtL2Y2NC9nZV9ub2FzbS5nbw==) | | |
| [mat/index\_bound\_checks.go](https://codecov.io/gh/gonum/gonum/pull/1653/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-bWF0L2luZGV4X2JvdW5kX2NoZWNrcy5nbw==) | | |
| [mat/offset\_appengine.go](https://codecov.io/gh/gonum/gonum/pull/1653/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-bWF0L29mZnNldF9hcHBlbmdpbmUuZ28=) | | |
| [internal/math32/sqrt.go](https://codecov.io/gh/gonum/gonum/pull/1653/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-aW50ZXJuYWwvbWF0aDMyL3NxcnQuZ28=) | | |
| [internal/asm/c128/stubs\_noasm.go](https://codecov.io/gh/gonum/gonum/pull/1653/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/dot.go](https://codecov.io/gh/gonum/gonum/pull/1653/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-aW50ZXJuYWwvYXNtL2Y2NC9kb3QuZ28=) | | |
| [internal/asm/c64/stubs\_noasm.go](https://codecov.io/gh/gonum/gonum/pull/1653/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-aW50ZXJuYWwvYXNtL2M2NC9zdHVic19ub2FzbS5nbw==) | | |
| [internal/asm/f64/axpy.go](https://codecov.io/gh/gonum/gonum/pull/1653/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum#diff-aW50ZXJuYWwvYXNtL2Y2NC9heHB5Lmdv) | | |
| ... and [8 more](https://codecov.io/gh/gonum/gonum/pull/1653/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/1653?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/1653?src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gonum). Last update [2b838f4...ff49396](https://codecov.io/gh/gonum/gonum/pull/1653?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).


Patch:
diff --git a/AUTHORS b/AUTHORS
index 198f685e..31a9d386 100644
--- a/AUTHORS
+++ b/AUTHORS
@@ -81,6 +81,7 @@ Nick Potts <nick@the-potts.com>
 Nils Wogatzky <odog@netcologne.de>
 Olivier Wulveryck <olivier.wulveryck@gmail.com>
 Or Rikon <rikonor@gmail.com>
+Patricio Whittingslow <graded.sp@gmail.com>
 Patrick DeVivo <patrick@tickgit.com>
 Pontus Melke <pontusmelke@gmail.com>
 Renee French
diff --git a/CONTRIBUTORS b/CONTRIBUTORS
index 259fa7cb..2777ee30 100644
--- a/CONTRIBUTORS
+++ b/CONTRIBUTORS
@@ -88,6 +88,7 @@ Nick Potts <nick@the-potts.com>
 Nils Wogatzky <odog@netcologne.de>
 Olivier Wulveryck <olivier.wulveryck@gmail.com>
 Or Rikon <rikonor@gmail.com>
+Patricio Whittingslow <graded.sp@gmail.com>
 Patrick DeVivo <patrick@tickgit.com>
 Pontus Melke <pontusmelke@gmail.com>
 Renee French
diff --git a/interp/cubic.go b/interp/cubic.go
index 8275ebb1..4afb5e59 100644
--- a/interp/cubic.go
+++ b/interp/cubic.go
@@ -104,8 +104,7 @@ func (pc *PiecewiseCubic) FitWithDerivatives(xs, ys, dydxs []float64) {
 		pc.coeffs.Set(i, 2, (3*dy-(2*dydxs[i]+dydxs[i+1])*dx)/dx/dx)
 		pc.coeffs.Set(i, 3, (-2*dy+(dydxs[i]+dydxs[i+1])*dx)/dx/dx/dx)
 	}
-	pc.xs = make([]float64, n)
-	copy(pc.xs, xs)
+	pc.xs = append(pc.xs[:0], xs...)
 	pc.lastY = ys[m]
 	pc.lastDyDx = dydxs[m]
 }
@@ -297,3 +296,119 @@ func fritschButlandEdgeDerivative(xs, ys, slopes []float64, leftEdge bool) float
 	}
 	return g
 }
+
+// fitWithSecondDerivatives fits a piecewise cubic predictor to (X, Y, d^2Y/dX^2) value
+// triples provided as three slices.
+// It panics if any of these is true:
+// - len(xs) < 2,
+// - elements of xs are not strictly increasing,
+// - len(xs) != len(ys),
+// - len(xs) != len(d2ydx2s).
+// Note that this method does not guarantee on its own the continuity of first derivatives.
+func (pc *PiecewiseCubic) fitWithSecondDerivatives(xs, ys, d2ydx2s []float64) {
+	n := len(xs)
+	switch {
+	case len(ys) != n, len(d2ydx2s) != n:
+		panic(differentLengths)
+	case n < 2:
+		panic(tooFewPoints)
+	}
+	m := n - 1
+	pc.coeffs.Reset()
+	pc.coeffs.ReuseAs(m, 4)
+	for i := 0; i < m; i++ {
+		dx := xs[i+1] - xs[i]
+		if dx <= 0 {
+			panic(xsNotStrictlyIncreasing)
+		}
+		dy := ys[i+1] - ys[i]
+		dm := d2ydx2s[i+1] - d2ydx2s[i]
+		pc.coeffs.Set(i, 0, ys[i])                             // a_0
+		pc.coeffs.Set(i, 1, (dy-(d2ydx2s[i]+dm/3)*dx*dx/2)/dx) // a_1
+		pc.coeffs.Set(i, 2, d2ydx2s[i]/2)                      // a_2
+		pc.coeffs.Set(i, 3, dm/6/dx)                           // a_3
+	}
+	pc.xs = append(pc.xs[:0], xs...)
+	pc.lastY = ys[m]
+	lastDx := xs[m] - xs[m-1]
+	pc.lastDyDx = pc.coeffs.At(m-1, 1) + 2*pc.coeffs.At(m-1, 2)*lastDx + 3*pc.coeffs.At(m-1, 3)*lastDx*lastDx
+}
+
+// makeCubicSplineSecondDerivativeEquations generates the basic system of linear equations
+// which have to be satisfied by the second derivatives to make the first derivatives of a
+// cubic spline continuous. It panics if elements of xs are not strictly increasing, or
+// len(xs) != len(ys).
+// makeCubicSplineSecondDerivativeEquations returns a tri-diagonal matrix A and a vector b
+// defining a system of linear equations A*m = b for second derivatives vector m.
+func makeCubicSplineSecondDerivativeEquations(xs, ys []float64) (*mat.Tridiag, mat.MutableVector) {
+	n := len(xs)
+	if len(ys) != n {
+		panic(differentLengths)
+	}
+	b := make([]float64, n)
+	m := n - 1
+	// Diagonal of A:
+	d := make([]float64, n)
+	// Sub-diagonal of A:
+	dl := make([]float64, m)
+	// Super-diagonal of A:
+	du := make([]float64, m)
+	if n > 2 {
+		for i := 0; i < m; i++ {
+			dx := xs[i+1] - xs[i]
+			if dx <= 0 {
+				panic(xsNotStrictlyIncreasing)
+			}
+			slope := (ys[i+1] - ys[i]) / dx
+			if i > 0 {
+				b[i] += slope
+				d[i] += dx / 3
+				du[i] = dx / 6
+			}
+			if i < m-1 {
+				b[i+1] -= slope
+				d[i+1] += dx / 3
+				dl[i] = dx / 6
+			}
+		}
+	}
+	return mat.NewTridiag(n, dl, d, du), mat.NewVecDense(n, b)
+}
+
+// NaturalCubic is a piecewise cubic 1-dimensional interpolator with
+// continuous value, first and second derivatives, which can be fitted to (X, Y)
+// value pairs without providing derivatives. See e.g. https://www.math.drexel.edu/~tolya/cubicspline.pdf
+// for details.
+type NaturalCubic struct {
+	cubic PiecewiseCubic
+}
+
+// Predict returns the interpolation value at x.
+func (nc *NaturalCubic) Predict(x float64) float64 {
+	return nc.cubic.Predict(x)
+}
+
+// PredictDerivative returns the predicted derivative at x.
+func (nc *NaturalCubic) PredictDerivative(x float64) float64 {
+	return nc.cubic.PredictDerivative(x)
+}
+
+// Fit fits a predictor to (X, Y) value pairs provided as two slices.
+// It panics if len(xs) < 2, elements of xs are not strictly increasing
+// or len(xs) != len(ys). It returns an error if solving the required system
+// of linear equations fails.
+func (nc *NaturalCubic) Fit(xs, ys []float64) error {
+	a, b := makeCubicSplineSecondDerivativeEquations(xs, ys)
+	// Add boundary conditions y''(left) = y''(right) = 0:
+	n := len(xs)
+	b.SetVec(0, 0)
+	b.SetVec(n-1, 0)
+	a.SetBand(0, 0, 1)
+	a.SetBand(n-1, n-1, 1)
+	x := mat.NewVecDense(n, nil)
+	err := a.SolveVecTo(x, false, b)
+	if err == nil {
+		nc.cubic.fitWithSecondDerivatives(xs, ys, x.RawVector().Data)
+	}
+	return err
+}
diff --git a/lapack/gonum/dgetc2.go b/lapack/gonum/dgetc2.go
new file mode 100644
index 00000000..2ba2a91e
--- /dev/null
+++ b/lapack/gonum/dgetc2.go
@@ -0,0 +1,117 @@
+// 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 gonum
+
+import (
+	"math"
+
+	"gonum.org/v1/gonum/blas/blas64"
+)
+
+// Dgetc2 computes an LU factorization with complete pivoting of the
+// n×n matrix A. The factorization has the form
+//  A = P * L * U * Q,
+// where P and Q are permutation matrices, L is lower triangular with
+// unit diagonal elements and U is upper triangular.
+//
+// a is modified to the information to construct L and U.
+// The lower triangle of a contains the matrix L (not including diagonal).
+// The upper triangle contains the matrix U. The matrices P and Q can
+// be constructed from ipiv and jpiv, respectively. k is non-negative if U(k, k)
+// is likely to produce overflow when we try to solve for x in Ax = b.
+// U is perturbed in this case to avoid the overflow.
+//
+// Dgetc2 is an internal routine. It is exported for testing purposes.
+func (impl Implementation) Dgetc2(n int, a []float64, lda int, ipiv, jpiv []int) (k int) {
+	switch {
+	case n < 0:
+		panic(nLT0)
+	case lda < max(1, n):
+		panic(badLdA)
+	}
+
+	// Negative k indicates U was not perturbed.
+	k = -1
+	// Quick return if possible.
+	if n == 0 {
+		return k
+	}
+
+	switch {
+	case len(a) < (n-1)*lda+n:
+		panic(shortA)
+	case len(ipiv) != n:
+		panic(badLenIpiv)
+	case len(jpiv) != n:
+		panic(badLenJpvt)
+	}
+
+	const (
+		eps    = dlamchP
+		smlnum = dlamchS / eps
+	)
+	if n == 1 {
+		ipiv[0], jpiv[0] = 0, 0
+		if math.Abs(a[0]) < smlnum {
+			a[0] = smlnum
+			k = 0
+		}
+		return k
+	}
+
+	// Factorize A using complete pivoting.
+	// Set pivots less than lc to lc.
+	var lc float64
+	var ipv, jpv int
+	bi := blas64.Implementation()
+	for i := 0; i < n-1; i++ {
+		xmax := 0.0
+		for ip := i; ip < n; ip++ {
+			for jp := i; jp < n; jp++ {
+				if math.Abs(a[ip*lda+jp]) >= xmax {
+					xmax = math.Abs(a[ip*lda+jp])
+					ipv = ip
+					jpv = jp
+				}
+			}
+		}
+		if i == 0 {
+			lc = math.Max(eps*xmax, smlnum)
+		}
+
+		// Swap rows.
+		if ipv != i {
+			bi.Dswap(n, a[ipv*lda:], 1, a[i*lda:], 1)
+		}
+		ipiv[i] = ipv
+
+		// Swap columns.
+		if jpv != i {
+			bi.Dswap(n, a[jpv:], lda, a[i:], lda)
+		}
+		jpiv[i] = jpv
+
+		// Check for singularity.
+		if math.Abs(a[i*lda+i]) < lc {
+			k = i
+			a[i*lda+i] = lc
+		}
+
+		for j := i + 1; j < n; j++ {
+			a[j*lda+i] /= a[i*lda+i]
+		}
+		bi.Dger(n-i-1, n-i-1, -1, a[(i+1)*lda+i:], lda, a[i*lda+i+1:], 1, a[(i+1)*lda+i+1:], lda)
+	}
+
+	if math.Abs(a[(n-1)*lda+n-1]) < lc {
+		k = n - 1
+		a[(n-1)*lda+(n-1)] = lc
+	}
+
+	// Set last pivots to last index.
+	ipiv[n-1] = n - 1
+	jpiv[n-1] = n - 1
+	return k
+}
diff --git a/lapack/gonum/dlaswp.go b/lapack/gonum/dlaswp.go
index b207d121..392477b4 100644
--- a/lapack/gonum/dlaswp.go
+++ b/lapack/gonum/dlaswp.go
@@ -42,11 +42,17 @@ func (impl Implementation) Dlaswp(n int, a []float64, lda int, k1, k2 int, ipiv
 	bi := blas64.Implementation()
 	if incX == 1 {
 		for k := k1; k <= k2; k++ {
+			if k == ipiv[k] {
+				continue
+			}
 			bi.Dswap(n, a[k*lda:], 1, a[ipiv[k]*lda:], 1)
 		}
 		return
 	}
 	for k := k2; k >= k1; k-- {
+		if k == ipiv[k] {
+			continue
+		}
 		bi.Dswap(n, a[k*lda:], 1, a[ipiv[k]*lda:], 1)
 	}
 }
diff --git a/lapack/gonum/dlasy2.go b/lapack/gonum/dlasy2.go
index d670f40d..d95402ca 100644
--- a/lapack/gonum/dlasy2.go
+++ b/lapack/gonum/dlasy2.go
@@ -261,7 +261,7 @@ func (impl Implementation) Dlasy2(tranl, tranr bool, isgn, n1, n2 int, tl []floa
 
 		maxbtmp := math.Max(math.Abs(btmp[0]), math.Abs(btmp[1]))
 		maxbtmp = math.Max(maxbtmp, math.Max(math.Abs(btmp[2]), math.Abs(btmp[3])))
-		scale = 1 / 8 / maxbtmp
+		scale = (1.0 / 8.0) / maxbtmp
 		btmp[0] *= scale
 		btmp[1] *= scale
 		btmp[2] *= scale
diff --git a/mat/band.go b/mat/band.go
index 25aca349..d6bfdfd9 100644
--- a/mat/band.go
+++ b/mat/band.go
@@ -51,6 +51,9 @@ type RawBander interface {
 // A MutableBanded can set elements of a band matrix.
 type MutableBanded interface {
 	Banded
+
+	// SetBand sets the element at row i, column j to the value v.
+	// It panics if the location is outside the appropriate region of the matrix.
 	SetBand(i, j int, v float64)
 }
 
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/interp/cubic_test.go b/interp/cubic_test.go
index 8f05aeee..b59df3e8 100644
--- a/interp/cubic_test.go
+++ b/interp/cubic_test.go
@@ -671,3 +671,134 @@ func TestFritschButlandErrors(t *testing.T) {
 		}
 	}
 }
+
+func TestPiecewiseCubicFitWithSecondDerivatives(t *testing.T) {
+	t.Parallel()
+	const tol = 1e-14
+	xs := []float64{-2, 0, 3}
+	ys := []float64{2.5, 1, 2.5}
+	d2ydx2s := []float64{1, 2, 3}
+	var pc PiecewiseCubic
+	pc.fitWithSecondDerivatives(xs, ys, d2ydx2s)
+	m := len(xs) - 1
+	if pc.lastY != ys[m] {
+		t.Errorf("Mismatch in lastY: got %v, want %g", pc.lastY, ys[m])
+	}
+	if !floats.Equal(pc.xs, xs) {
+		t.Errorf("Mismatch in xs: got %v, want %v", pc.xs, xs)
+	}
+	for i := 0; i < len(xs); i++ {
+		yHat := pc.Predict(xs[i])
+		if math.Abs(yHat-ys[i]) > tol {
+			t.Errorf("Mismatch in predicted Y[%d]: got %v, want %g", i, yHat, ys[i])
+		}
+		var d2ydx2Hat float64
+		if i < m {
+			d2ydx2Hat = 2 * pc.coeffs.At(i, 2)
+		} else {
+			d2ydx2Hat = 2*pc.coeffs.At(m-1, 2) + 6*pc.coeffs.At(m-1, 3)*(xs[m]-xs[m-1])
+		}
+		if math.Abs(d2ydx2Hat-d2ydx2s[i]) > tol {
+			t.Errorf("Mismatch in predicted d2Y/dX2[%d]: got %v, want %g", i, d2ydx2Hat, d2ydx2s[i])
+		}
+	}
+	// Test pc.lastDyDx without copying verbatim the calculation from the tested method:
+	lastDyDx := pc.PredictDerivative(xs[m] - tol/1000)
+	if math.Abs(lastDyDx-pc.lastDyDx) > tol {
+		t.Errorf("Mismatch in lastDxDy: got %v, want %g", pc.lastDyDx, lastDyDx)
+	}
+}
+
+func TestPiecewiseCubicFitWithSecondDerivativesErrors(t *testing.T) {
+	t.Parallel()
+	for _, test := range []struct {
+		xs, ys, d2ydx2s []float64
+	}{
+		{
+			xs:      []float64{0, 1, 2},
+			ys:      []float64{10, 20},
+			d2ydx2s: []float64{0, 0, 0},
+		},
+		{
+			xs:      []float64{0, 1, 1},
+			ys:      []float64{10, 20, 30},
+			d2ydx2s: []float64{0, 0, 0, 0},
+		},
+		{
+			xs:      []float64{0},
+			ys:      []float64{0},
+			d2ydx2s: []float64{0},
+		},
+		{
+			xs:      []float64{0, 1, 1},
+			ys:      []float64{10, 20, 10},
+			d2ydx2s: []float64{0, 0, 0},
+		},
+	} {
+		var pc PiecewiseCubic
+		if !panics(func() { pc.fitWithSecondDerivatives(test.xs, test.ys, test.d2ydx2s) }) {
+			t.Errorf("expected panic for xs: %v, ys: %v and d2ydx2s: %v", test.xs, test.ys, test.d2ydx2s)
+		}
+	}
+}
+
+func TestMakeCubicSplineSecondDerivativeEquations(t *testing.T) {
+	t.Parallel()
+	const tol = 1e-15
+	xs := []float64{-1, 0, 2}
+	ys := []float64{2, 0, 2}
+	n := len(xs)
+	A, b := makeCubicSplineSecondDerivativeEquations(xs, ys)
+	if b.Len() != n {
+		t.Errorf("Mismatch in b size: got %v, want %d", b.Len(), n)
+	}
+	r, c := A.Dims()
+	if r != n || c != n {
+		t.Errorf("Mismatch in A size: got %d x %d, want %d x %d", r, c, n, n)
+	}
+	expectedB := mat.NewVecDense(3, []float64{0, 3, 0})
+	var diffB mat.VecDense
+	diffB.SubVec(b, expectedB)
+	if diffB.Norm(math.Inf(1)) > tol {
+		t.Errorf("Mismatch in b values: got %v, want %v", b, expectedB)
+	}
+	expectedA := mat.NewDense(3, 3, []float64{0, 0, 0, 1 / 6., 1, 2 / 6., 0, 0, 0})
+	var diffA mat.Dense
+	diffA.Sub(A, expectedA)
+	if diffA.Norm(math.Inf(1)) > tol {
+		t.Errorf("Mismatch in A values: got %v, want %v", A, expectedA)
+	}
+}
+
+func TestNaturalCubicFit(t *testing.T) {
+	t.Parallel()
+	xs := []float64{-1, 0, 2, 3.5}
+	ys := []float64{2, 0, 2, 1.5}
+	var nc NaturalCubic
+	err := nc.Fit(xs, ys)
+	if err != nil {
+		t.Errorf("Error when fitting NaturalCubic: %v", err)
+	}
+	testXs := []float64{-1, -0.99, -0.5, 0, 0.5, 1, 1.5, 2, 2.5, 3, 3.49, 3.5}
+	// From scipy.interpolate.CubicSpline:
+	want := []float64{
+		2.0,
+		1.9737725526315788,
+		0.7664473684210527,
+		0.0,
+		0.027960526315789477,
+		0.6184210526315789,
+		1.3996710526315788,
+		2.0,
+		2.1403508771929824,
+		1.9122807017543857,
+		1.508859403508772,
+		1.5}
+	for i := 0; i < len(testXs); i++ {
+		got := nc.Predict(testXs[i])
+		if math.Abs(got-want[i]) > 1e-14 {
+			t.Errorf("Mismatch in predicted Y value for x = %g: got %v, want %g", testXs[i], got, want[i])
+		}
+	}
+
+}
diff --git a/lapack/gonum/lapack_test.go b/lapack/gonum/lapack_test.go
index 4d76525b..2dded735 100644
--- a/lapack/gonum/lapack_test.go
+++ b/lapack/gonum/lapack_test.go
@@ -118,6 +118,11 @@ func TestDgesvd(t *testing.T) {
 	testlapack.DgesvdTest(t, impl, tol)
 }
 
+func TestDgetc2(t *testing.T) {
+	t.Parallel()
+	testlapack.Dgetc2Test(t, impl)
+}
+
 func TestDgetri(t *testing.T) {
 	t.Parallel()
 	testlapack.DgetriTest(t, impl)
diff --git a/lapack/testlapack/dgetc2.go b/lapack/testlapack/dgetc2.go
new file mode 100644
index 00000000..461c28df
--- /dev/null
+++ b/lapack/testlapack/dgetc2.go
@@ -0,0 +1,98 @@
+// 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 testlapack
+
+import (
+	"fmt"
+	"math"
+	"testing"
+
+	"golang.org/x/exp/rand"
+
+	"gonum.org/v1/gonum/blas"
+	"gonum.org/v1/gonum/blas/blas64"
+)
+
+type Dgetc2er interface {
+	Dgetc2(n int, a []float64, lda int, ipiv, jpiv []int) (k int)
+}
+
+func Dgetc2Test(t *testing.T, impl Dgetc2er) {
+	const tol = 1e-12
+	rnd := rand.New(rand.NewSource(1))
+	for _, n := range []int{0, 1, 2, 3, 4, 5, 10, 20} {
+		for _, lda := range []int{n, n + 5} {
+			dgetc2Test(t, impl, rnd, n, lda, tol)
+		}
+	}
+}
+
+func dgetc2Test(t *testing.T, impl Dgetc2er, rnd *rand.Rand, n, lda int, tol float64) {
+	name := fmt.Sprintf("n=%v,lda=%v", n, lda)
+	if lda == 0 {
+		lda = 1
+	}
+	// Generate a random general matrix A.
+	a := randomGeneral(n, n, lda, rnd)
+	// ipiv and jpiv are outputs.
+	ipiv := make([]int, n)
+	jpiv := make([]int, n)
+	for i := 0; i < n; i++ {
+		ipiv[i], jpiv[i] = -1, -1 // Set to non-indices.
+	}
+	// Copy to store output (LU decomposition).
+	lu := cloneGeneral(a)
+	k := impl.Dgetc2(n, lu.Data, lu.Stride, ipiv, jpiv)
+	if k >= 0 {
+		t.Logf("%v: matrix was perturbed at %d", name, k)
+	}
+
+	// Verify all indices are set.
+	for i := 0; i < n; i++ {
+		if ipiv[i] < 0 {
+			t.Errorf("%v: ipiv[%d] is negative", name, i)
+		}
+		if jpiv[i] < 0 {
+			t.Errorf("%v: jpiv[%d] is negative", name, i)
+		}
+	}
+	bi := blas64.Implementation()
+	// Construct L and U triangular matrices from Dgetc2 output.
+	L := zeros(n, n, lda)
+	U := zeros(n, n, lda)
+	for i := 0; i < n; i++ {
+		for j := 0; j < n; j++ {
+			idx := i*lda + j
+			if j >= i { // On upper triangle and setting of L's unit diagonal elements.
+				U.Data[idx] = lu.Data[idx]
+				if j == i {
+					L.Data[idx] = 1.0
+				}
+			} else if i > j { // On diagonal or lower triangle.
+				L.Data[idx] = lu.Data[idx]
+			}
+		}
+	}
+	work := zeros(n, n, lda)
+	bi.Dgemm(blas.NoTrans, blas.NoTrans, n, n, n, 1, L.Data, L.Stride, U.Data, U.Stride, 0, work.Data, work.Stride)
+
+	// Apply Permutations P and Q to L*U.
+	for i := n - 1; i >= 0; i-- {
+		ipv, jpv := ipiv[i], jpiv[i]
+		if ipv != i {
+			bi.Dswap(n, work.Data[i*lda:], 1, work.Data[ipv*lda:], 1)
+		}
+		if jpv != i {
+			bi.Dswap(n, work.Data[i:], work.Stride, work.Data[jpv:], work.Stride)
+		}
+	}
+
+	// A should be reconstructed by now.
+	for i := range work.Data {
+		if math.Abs(work.Data[i]-a.Data[i]) > tol {
+			t.Errorf("%v: matrix %d idx not equal after reconstruction. got %g, expected %g", name, i, work.Data[i], a.Data[i])
+		}
+	}
+}
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)
+	}
+}
diff --git a/stat/distuv/distribution_test.go b/stat/distuv/distribution_test.go
index 777aa4ab..a63fad8b 100644
--- a/stat/distuv/distribution_test.go
+++ b/stat/distuv/distribution_test.go
@@ -70,43 +70,48 @@ type cumulantProber interface {
 	probLogprober
 }
 
-func checkMean(t *testing.T, i int, x []float64, m meaner, tol float64) {
+func checkMean(t *testing.T, cas int, x []float64, m meaner, tol float64) {
+	t.Helper()
 	mean := stat.Mean(x, nil)
 	if !scalar.EqualWithinAbsOrRel(mean, m.Mean(), tol, tol) {
-		t.Errorf("Mean mismatch case %v: want: %v, got: %v", i, mean, m.Mean())
+		t.Errorf("Mean mismatch case %v: want: %v, got: %v", cas, mean, m.Mean())
 	}
 }
 
-func checkMedian(t *testing.T, i int, x []float64, m medianer, tol float64) {
+func checkMedian(t *testing.T, cas int, x []float64, m medianer, tol float64) {
+	t.Helper()
 	median := stat.Quantile(0.5, stat.Empirical, x, nil)
 	if !scalar.EqualWithinAbsOrRel(median, m.Median(), tol, tol) {
-		t.Errorf("Median mismatch case %v: want: %v, got: %v", i, median, m.Median())
+		t.Errorf("Median mismatch case %v: want: %v, got: %v", cas, median, m.Median())
 	}
 }
 
-func checkVarAndStd(t *testing.T, i int, x []float64, v varStder, tol float64) {
+func checkVarAndStd(t *testing.T, cas int, x []float64, v varStder, tol float64) {
+	t.Helper()
 	variance := stat.Variance(x, nil)
 	if !scalar.EqualWithinAbsOrRel(variance, v.Variance(), tol, tol) {
-		t.Errorf("Variance mismatch case %v: want: %v, got: %v", i, variance, v.Variance())
+		t.Errorf("Variance mismatch case %v: want: %v, got: %v", cas, variance, v.Variance())
 	}
 	std := math.Sqrt(variance)
 	if !scalar.EqualWithinAbsOrRel(std, v.StdDev(), tol, tol) {
-		t.Errorf("StdDev mismatch case %v: want: %v, got: %v", i, std, v.StdDev())
+		t.Errorf("StdDev mismatch case %v: want: %v, got: %v", cas, std, v.StdDev())
 	}
 }
 
-func checkEntropy(t *testing.T, i int, x []float64, e entropyer, tol float64) {
+func checkEntropy(t *testing.T, cas int, x []float64, e entropyer, tol float64) {
+	t.Helper()
 	tmp := make([]float64, len(x))
 	for i, v := range x {
 		tmp[i] = -e.LogProb(v)
 	}
 	entropy := stat.Mean(tmp, nil)
 	if !scalar.EqualWithinAbsOrRel(entropy, e.Entropy(), tol, tol) {
-		t.Errorf("Entropy mismatch case %v: want: %v, got: %v", i, entropy, e.Entropy())
+		t.Errorf("Entropy mismatch case %v: want: %v, got: %v", cas, entropy, e.Entropy())
 	}
 }
 
-func checkExKurtosis(t *testing.T, i int, x []float64, e exKurtosiser, tol float64) {
+func checkExKurtosis(t *testing.T, cas int, x []float64, e exKurtosiser, tol float64) {
+	t.Helper()
 	mean := e.Mean()
 	tmp := make([]float64, len(x))
 	for i, x := range x {
@@ -116,11 +121,12 @@ func checkExKurtosis(t *testing.T, i int, x []float64, e exKurtosiser, tol float
 	mu4 := stat.Mean(tmp, nil)
 	kurtosis := mu4/(variance*variance) - 3
 	if !scalar.EqualWithinAbsOrRel(kurtosis, e.ExKurtosis(), tol, tol) {
-		t.Errorf("ExKurtosis mismatch case %v: want: %v, got: %v", i, kurtosis, e.ExKurtosis())
+		t.Errorf("ExKurtosis mismatch case %v: want: %v, got: %v", cas, kurtosis, e.ExKurtosis())
 	}
 }
 
-func checkSkewness(t *testing.T, i int, x []float64, s skewnesser, tol float64) {
+func checkSkewness(t *testing.T, cas int, x []float64, s skewnesser, tol float64) {
+	t.Helper()
 	mean := s.Mean()
 	std := s.StdDev()
 	tmp := make([]float64, len(x))
@@ -130,24 +136,25 @@ func checkSkewness(t *testing.T, i int, x []float64, s skewnesser, tol float64)
 	mu3 := stat.Mean(tmp, nil)
 	skewness := mu3 / math.Pow(std, 3)
 	if !scalar.EqualWithinAbsOrRel(skewness, s.Skewness(), tol, tol) {
-		t.Errorf("Skewness mismatch case %v: want: %v, got: %v", i, skewness, s.Skewness())
+		t.Errorf("Skewness mismatch case %v: want: %v, got: %v", cas, skewness, s.Skewness())
 	}
 }
 
-func checkQuantileCDFSurvival(t *testing.T, i int, xs []float64, c cumulanter, tol float64) {
+func checkQuantileCDFSurvival(t *testing.T, cas int, xs []float64, c cumulanter, tol float64) {
+	t.Helper()
 	// Quantile, CDF, and survival check.
-	for i, p := range []float64{0.1, 0.25, 0.5, 0.75, 0.9} {
+	for _, p := range []float64{0.1, 0.25, 0.5, 0.75, 0.9} {
 		x := c.Quantile(p)
 		cdf := c.CDF(x)
 		estCDF := stat.CDF(x, stat.Empirical, xs, nil)
 		if !scalar.EqualWithinAbsOrRel(cdf, estCDF, tol, tol) {
-			t.Errorf("CDF mismatch case %v: want: %v, got: %v", i, estCDF, cdf)
+			t.Errorf("CDF mismatch case %v: want: %v, got: %v", cas, estCDF, cdf)
 		}
 		if !scalar.EqualWithinAbsOrRel(cdf, p, tol, tol) {
-			t.Errorf("Quantile/CDF mismatch case %v: want: %v, got: %v", i, p, cdf)
+			t.Errorf("Quantile/CDF mismatch case %v: want: %v, got: %v", cas, p, cdf)
 		}
 		if math.Abs(1-cdf-c.Survival(x)) > 1e-14 {
-			t.Errorf("Survival/CDF mismatch case %v: want: %v, got: %v", i, 1-cdf, c.Survival(x))
+			t.Errorf("Survival/CDF mismatch case %v: want: %v, got: %v", cas, 1-cdf, c.Survival(x))
 		}
 	}
 	if !panics(func() { c.Quantile(-0.0001) }) {
@@ -160,16 +167,17 @@ func checkQuantileCDFSurvival(t *testing.T, i int, xs []float64, c cumulanter, t
 
 // checkProbContinuous checks that the PDF is consistent with LogPDF
 // and integrates to 1 from the lower to upper bound.
-func checkProbContinuous(t *testing.T, i int, x []float64, lower float64, upper float64, p probLogprober, tol float64) {
+func checkProbContinuous(t *testing.T, cas int, x []float64, lower float64, upper float64, p probLogprober, tol float64) {
+	t.Helper()
 	q := quad.Fixed(p.Prob, lower, upper, 1000000, nil, 0)
 	if math.Abs(q-1) > tol {
-		t.Errorf("Probability distribution doesn't integrate to 1. Case %v: Got %v", i, q)
+		t.Errorf("Probability distribution doesn't integrate to 1. Case %v: Got %v", cas, q)
 	}
 
 	// Check that PDF and LogPDF are consistent.
-	for i, v := range x {
+	for _, v := range x {
 		if math.Abs(math.Log(p.Prob(v))-p.LogProb(v)) > 1e-14 {
-			t.Errorf("Prob and LogProb mismatch case %v at %v: want %v, got %v", i, v, math.Log(v), p.LogProb(v))
+			t.Errorf("Prob and LogProb mismatch case %v at %v: want %v, got %v", cas, v, math.Log(v), p.LogProb(v))
 			break
 		}
 	}
@@ -180,7 +188,8 @@ func checkProbContinuous(t *testing.T, i int, x []float64, lower float64, upper
 // to 1 and greater than 0). However, this is also true if the PDF of a different
 // distribution is used. This checks that PDF is also consistent with the
 // CDF implementation and the random samples.
-func checkProbQuantContinuous(t *testing.T, i int, xs []float64, c cumulantProber, tol float64) {
+func checkProbQuantContinuous(t *testing.T, cas int, xs []float64, c cumulantProber, tol float64) {
+	t.Helper()
 	ps := make([]float64, 101)
 	floats.Span(ps, 0, 1)
 
@@ -190,14 +199,14 @@ func checkProbQuantContinuous(t *testing.T, i int, xs []float64, c cumulantProbe
 		if p == 0 {
 			xp = x
 			if floats.Min(xs) < x {
-				t.Errorf("Sample of x less than Quantile(0). Case %v.", i)
+				t.Errorf("Sample of x less than Quantile(0). Case %v.", cas)
 				break
 			}
 			continue
 		}
 		if p == 1 {
 			if floats.Max(xs) > x {
-				t.Errorf("Sample of x greater than Quantile(1). Case %v.", i)
+				t.Errorf("Sample of x greater than Quantile(1). Case %v.", cas)
 				break
 			}
 		}
@@ -206,13 +215,13 @@ func checkProbQuantContinuous(t *testing.T, i int, xs []float64, c cumulantProbe
 		// the quantiles.
 		q := quad.Fixed(c.Prob, xp, x, 1000, nil, 0)
 		if math.Abs(q-(p-ps[i-1])) > 1e-5 {
-			t.Errorf("Integral of PDF doesn't match quantile. Case %v. Want %v, got %v.", i, p-ps[i-1], q)
+			t.Errorf("Integral of PDF doesn't match quantile. Case %v. Want %v, got %v.", cas, p-ps[i-1], q)
 			break
 		}
 
 		pEst := stat.CDF(x, stat.Empirical, xs, nil)
 		if math.Abs(pEst-p) > tol {
-			t.Errorf("Empirical CDF doesn't match quantile. Case %v.", i)
+			t.Errorf("Empirical CDF doesn't match quantile. Case %v.", cas)
 		}
 		xp = x
 	}
@@ -222,7 +231,8 @@ type moder interface {
 	Mode() float64
 }
 
-func checkMode(t *testing.T, i int, xs []float64, m moder, dx float64, tol float64) {
+func checkMode(t *testing.T, cas int, xs []float64, m moder, dx float64, tol float64) {
+	t.Helper()
 	rXs := make([]float64, len(xs))
 	for j, x := range xs {
 		rXs[j] = math.RoundToEven(x/dx) * dx
@@ -230,12 +240,13 @@ func checkMode(t *testing.T, i int, xs []float64, m moder, dx float64, tol float
 	want, _ := stat.Mode(rXs, nil)
 	got := m.Mode()
 	if !scalar.EqualWithinAbs(want, got, tol) {
-		t.Errorf("Mode mismatch case %d: want %g, got %v", i, want, got)
+		t.Errorf("Mode mismatch case %d: want %g, got %v", cas, want, got)
 	}
 }
 
 // checkProbDiscrete confirms that PDF and Rand are consistent for discrete distributions.
-func checkProbDiscrete(t *testing.T, i int, xs []float64, p probLogprober, tol float64) {
+func checkProbDiscrete(t *testing.T, cas int, xs []float64, p probLogprober, tol float64) {
+	t.Helper()
 	// Make a map of all of the unique samples.
 	m := make(map[float64]int)
 	for _, v := range xs {
@@ -244,17 +255,18 @@ func checkProbDiscrete(t *testing.T, i int, xs []float64, p probLogprober, tol f
 	for x, count := range m {
 		prob := float64(count) / float64(len(xs))
 		if math.Abs(prob-p.Prob(x)) > tol {
-			t.Errorf("PDF mismatch case %v at %v: want %v, got %v", i, x, prob, p.Prob(x))
+			t.Errorf("PDF mismatch case %v at %v: want %v, got %v", cas, x, prob, p.Prob(x))
 		}
 		if math.Abs(math.Log(p.Prob(x))-p.LogProb(x)) > 1e-14 {
-			t.Errorf("Prob and LogProb mismatch case %v at %v: want %v, got %v", i, x, math.Log(x), p.LogProb(x))
+			t.Errorf("Prob and LogProb mismatch case %v at %v: want %v, got %v", cas, x, math.Log(x), p.LogProb(x))
 		}
 	}
 }
 
 // testRandLogProb tests that LogProb and Rand give consistent results. This
 // can be used when the distribution does not implement CDF.
-func testRandLogProbContinuous(t *testing.T, i int, min float64, x []float64, f LogProber, tol float64, bins int) {
+func testRandLogProbContinuous(t *testing.T, cas int, min float64, x []float64, f LogProber, tol float64, bins int) {
+	t.Helper()
 	for cdf := 1 / float64(bins); cdf <= 1-1/float64(bins); cdf += 1 / float64(bins) {
 		// Get the estimated CDF from the samples
 		pt := stat.Quantile(cdf, stat.Empirical, x, nil)
@@ -265,7 +277,7 @@ func testRandLogProbContinuous(t *testing.T, i int, min float64, x []float64, f
 		// Integrate the PDF to find the CDF
 		estCDF := quad.Fixed(prob, min, pt, 10000, nil, 0)
 		if !scalar.EqualWithinAbsOrRel(cdf, estCDF, tol, tol) {
-			t.Errorf("Mismatch between integral of PDF and empirical CDF. Case %v. Want %v, got %v", i, cdf, estCDF)
+			t.Errorf("Mismatch between integral of PDF and empirical CDF. Case %v. Want %v, got %v", cas, cdf, estCDF)
 		}
 	}
 }