internal/testmath: add two-sample Welch's t-test for performance tests

This CL copies code from github.com/aclements/go-moremath/stats and
github.com/aclements/go-moremath/mathx for Welch's t-test. Several
existing tests in the Go repository check performance and scalability,
and this import is part of a move toward a more rigorous measurement of
both.

Note that the copied code is already licensed to Go Authors, so there's
no need to worry about additional licensing considerations.

For #32986.

Change-Id: I058630fab7216d1a589bb182b69fa2231e6f5475
Reviewed-on: https://go-review.googlesource.com/c/go/+/411395
Reviewed-by: Michael Pratt <mpratt@google.com>
This commit is contained in:
Michael Anthony Knyszek 2022-06-09 20:41:18 +00:00 committed by Michael Knyszek
parent 24b9039149
commit 6130461149
2 changed files with 217 additions and 1 deletions

View file

@ -543,7 +543,10 @@ var depsRules = `
internal/fuzz, internal/testlog, runtime/pprof, regexp
< testing/internal/testdeps;
OS, flag, testing, internal/cfg
MATH, errors, testing
< internal/testmath;
OS, flag, testing, internal/cfg, internal/testmath
< internal/testenv;
OS, encoding/base64

View file

@ -0,0 +1,213 @@
// Copyright 2022 The Go 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 testmath
import (
"errors"
"math"
)
// A TTestSample is a sample that can be used for a one or two sample
// t-test.
type TTestSample interface {
Weight() float64
Mean() float64
Variance() float64
}
var (
ErrSampleSize = errors.New("sample is too small")
ErrZeroVariance = errors.New("sample has zero variance")
ErrMismatchedSamples = errors.New("samples have different lengths")
)
// TwoSampleWelchTTest performs a two-sample (unpaired) Welch's t-test
// on samples x1 and x2. This t-test does not assume the distributions
// have equal variance.
func TwoSampleWelchTTest(x1, x2 TTestSample, alt LocationHypothesis) (*TTestResult, error) {
n1, n2 := x1.Weight(), x2.Weight()
if n1 <= 1 || n2 <= 1 {
// TODO: Can we still do this with n == 1?
return nil, ErrSampleSize
}
v1, v2 := x1.Variance(), x2.Variance()
if v1 == 0 && v2 == 0 {
return nil, ErrZeroVariance
}
dof := math.Pow(v1/n1+v2/n2, 2) /
(math.Pow(v1/n1, 2)/(n1-1) + math.Pow(v2/n2, 2)/(n2-1))
s := math.Sqrt(v1/n1 + v2/n2)
t := (x1.Mean() - x2.Mean()) / s
return newTTestResult(int(n1), int(n2), t, dof, alt), nil
}
// A TTestResult is the result of a t-test.
type TTestResult struct {
// N1 and N2 are the sizes of the input samples. For a
// one-sample t-test, N2 is 0.
N1, N2 int
// T is the value of the t-statistic for this t-test.
T float64
// DoF is the degrees of freedom for this t-test.
DoF float64
// AltHypothesis specifies the alternative hypothesis tested
// by this test against the null hypothesis that there is no
// difference in the means of the samples.
AltHypothesis LocationHypothesis
// P is p-value for this t-test for the given null hypothesis.
P float64
}
func newTTestResult(n1, n2 int, t, dof float64, alt LocationHypothesis) *TTestResult {
dist := TDist{dof}
var p float64
switch alt {
case LocationDiffers:
p = 2 * (1 - dist.CDF(math.Abs(t)))
case LocationLess:
p = dist.CDF(t)
case LocationGreater:
p = 1 - dist.CDF(t)
}
return &TTestResult{N1: n1, N2: n2, T: t, DoF: dof, AltHypothesis: alt, P: p}
}
// A LocationHypothesis specifies the alternative hypothesis of a
// location test such as a t-test or a Mann-Whitney U-test. The
// default (zero) value is to test against the alternative hypothesis
// that they differ.
type LocationHypothesis int
const (
// LocationLess specifies the alternative hypothesis that the
// location of the first sample is less than the second. This
// is a one-tailed test.
LocationLess LocationHypothesis = -1
// LocationDiffers specifies the alternative hypothesis that
// the locations of the two samples are not equal. This is a
// two-tailed test.
LocationDiffers LocationHypothesis = 0
// LocationGreater specifies the alternative hypothesis that
// the location of the first sample is greater than the
// second. This is a one-tailed test.
LocationGreater LocationHypothesis = 1
)
// A TDist is a Student's t-distribution with V degrees of freedom.
type TDist struct {
V float64
}
// PDF returns the value at x of the probability distribution function for the
// distribution.
func (t TDist) PDF(x float64) float64 {
return math.Exp(lgamma((t.V+1)/2)-lgamma(t.V/2)) /
math.Sqrt(t.V*math.Pi) * math.Pow(1+(x*x)/t.V, -(t.V+1)/2)
}
// CDF returns the value at x of the cumulative distribution function for the
// distribution.
func (t TDist) CDF(x float64) float64 {
if x == 0 {
return 0.5
} else if x > 0 {
return 1 - 0.5*betaInc(t.V/(t.V+x*x), t.V/2, 0.5)
} else if x < 0 {
return 1 - t.CDF(-x)
} else {
return math.NaN()
}
}
func (t TDist) Bounds() (float64, float64) {
return -4, 4
}
func lgamma(x float64) float64 {
y, _ := math.Lgamma(x)
return y
}
// betaInc returns the value of the regularized incomplete beta
// function Iₓ(a, b) = 1 / B(a, b) * ∫₀ˣ tᵃ⁻¹ (1-t)ᵇ⁻¹ dt.
//
// This is not to be confused with the "incomplete beta function",
// which can be computed as BetaInc(x, a, b)*Beta(a, b).
//
// If x < 0 or x > 1, returns NaN.
func betaInc(x, a, b float64) float64 {
// Based on Numerical Recipes in C, section 6.4. This uses the
// continued fraction definition of I:
//
// (xᵃ*(1-x)ᵇ)/(a*B(a,b)) * (1/(1+(d₁/(1+(d₂/(1+...))))))
//
// where B(a,b) is the beta function and
//
// d_{2m+1} = -(a+m)(a+b+m)x/((a+2m)(a+2m+1))
// d_{2m} = m(b-m)x/((a+2m-1)(a+2m))
if x < 0 || x > 1 {
return math.NaN()
}
bt := 0.0
if 0 < x && x < 1 {
// Compute the coefficient before the continued
// fraction.
bt = math.Exp(lgamma(a+b) - lgamma(a) - lgamma(b) +
a*math.Log(x) + b*math.Log(1-x))
}
if x < (a+1)/(a+b+2) {
// Compute continued fraction directly.
return bt * betacf(x, a, b) / a
} else {
// Compute continued fraction after symmetry transform.
return 1 - bt*betacf(1-x, b, a)/b
}
}
// betacf is the continued fraction component of the regularized
// incomplete beta function Iₓ(a, b).
func betacf(x, a, b float64) float64 {
const maxIterations = 200
const epsilon = 3e-14
raiseZero := func(z float64) float64 {
if math.Abs(z) < math.SmallestNonzeroFloat64 {
return math.SmallestNonzeroFloat64
}
return z
}
c := 1.0
d := 1 / raiseZero(1-(a+b)*x/(a+1))
h := d
for m := 1; m <= maxIterations; m++ {
mf := float64(m)
// Even step of the recurrence.
numer := mf * (b - mf) * x / ((a + 2*mf - 1) * (a + 2*mf))
d = 1 / raiseZero(1+numer*d)
c = raiseZero(1 + numer/c)
h *= d * c
// Odd step of the recurrence.
numer = -(a + mf) * (a + b + mf) * x / ((a + 2*mf) * (a + 2*mf + 1))
d = 1 / raiseZero(1+numer*d)
c = raiseZero(1 + numer/c)
hfac := d * c
h *= hfac
if math.Abs(hfac-1) < epsilon {
return h
}
}
panic("betainc: a or b too big; failed to converge")
}