1007 lines
40 KiB
Go
1007 lines
40 KiB
Go
// Copyright (c) 2026 Petr Balvín <opensource@petrbalvin.org> (https://petrbalvin.org)
|
|
// SPDX-License-Identifier: MIT
|
|
|
|
// Behaviour pins for the package's entry points: every test here
|
|
// holds a corrected defect to the answer it must keep, among them the
|
|
// NaN quantile hang, the rolling dtype panics, the Quantile NaN
|
|
// panic, the Histogram2D complex panic and unbounded bin counts, the
|
|
// lower-half NormalQuantile bracket, the saturated GLM likelihood,
|
|
// the weighted regression intercept statistics, the missing
|
|
// finiteness and symmetry refusals, the BinomialDraws count guard,
|
|
// the Kolmogorov series' truncation, the NaN-parameter and
|
|
// complex-input sweeps, the KolmogorovSmirnovTest NaN loop and the
|
|
// degenerate-response F statistic. Each one fails with the old
|
|
// behaviour restored.
|
|
|
|
package stats
|
|
|
|
import (
|
|
"math"
|
|
"math/big"
|
|
"testing"
|
|
"time"
|
|
|
|
"sourcedock.dev/petrbalvin/tensor/internal/core"
|
|
)
|
|
|
|
// pinVector builds a rank-1 float array.
|
|
func pinVector(t *testing.T, vals ...float64) *core.Array {
|
|
t.Helper()
|
|
a, err := core.FromFloats(vals, len(vals))
|
|
if err != nil {
|
|
t.Fatalf("FromFloats(%v): %v", vals, err)
|
|
}
|
|
return a
|
|
}
|
|
|
|
// pinWatchdog runs f in its own goroutine and fails the test if it
|
|
// has not returned within limit. It is what keeps the hang regression a
|
|
// failure rather than a stuck test binary.
|
|
func pinWatchdog(t *testing.T, limit time.Duration, what string, f func()) {
|
|
t.Helper()
|
|
done := make(chan struct{})
|
|
go func() {
|
|
defer close(done)
|
|
f()
|
|
}()
|
|
select {
|
|
case <-done:
|
|
case <-time.After(limit):
|
|
t.Fatalf("%s has not returned after %v: the loop does not terminate", what, limit)
|
|
}
|
|
}
|
|
|
|
// TestDiscreteQuantileNaNQuantileDoesNotHang pins the P0: q = NaN passed
|
|
// the written-out range guard (`q < 0 || q > 1` is false for NaN), every
|
|
// bracket comparison below was then false too, and the doubling loop
|
|
// never terminated. Both discrete quantiles must refuse NaN instead.
|
|
func TestDiscreteQuantileNaNQuantileDoesNotHang(t *testing.T) {
|
|
for _, tc := range []struct {
|
|
name string
|
|
call func() (float64, error)
|
|
}{
|
|
{"PoissonQuantile", func() (float64, error) { return PoissonQuantile(math.NaN(), 3) }},
|
|
{"BinomialQuantile", func() (float64, error) { return BinomialQuantile(math.NaN(), 0.5, 5) }},
|
|
} {
|
|
t.Run(tc.name, func(t *testing.T) {
|
|
var err error
|
|
pinWatchdog(t, 10*time.Second, tc.name+"(NaN)", func() {
|
|
_, err = tc.call()
|
|
})
|
|
if err == nil {
|
|
t.Fatalf("%s(NaN) returned a value, want a refusal", tc.name)
|
|
}
|
|
})
|
|
}
|
|
}
|
|
|
|
// bigSqrt2 is sqrt(2) at the working precision.
|
|
func bigSqrt2() *big.Float {
|
|
const prec = 256
|
|
return new(big.Float).SetPrec(prec).Sqrt(new(big.Float).SetPrec(prec).SetInt64(2))
|
|
}
|
|
|
|
// pinIntArray builds a rank-1 int array.
|
|
func pinIntArray(t *testing.T, vals ...int64) *core.Array {
|
|
t.Helper()
|
|
a, err := core.FromInts(vals, len(vals))
|
|
if err != nil {
|
|
t.Fatalf("FromInts(%v): %v", vals, err)
|
|
}
|
|
return a
|
|
}
|
|
|
|
// pinFloat32Array builds a rank-1 float32 array.
|
|
func pinFloat32Array(t *testing.T, vals ...float32) *core.Array {
|
|
t.Helper()
|
|
a, err := core.FromFloat32s(vals, len(vals))
|
|
if err != nil {
|
|
t.Fatalf("FromFloat32s(%v): %v", vals, err)
|
|
}
|
|
return a
|
|
}
|
|
|
|
// pinSameFloats pins two float slices to the bit for the dtype tests.
|
|
func pinSameFloats(t *testing.T, label string, got, want []float64) {
|
|
t.Helper()
|
|
if len(got) != len(want) {
|
|
t.Fatalf("%s has length %d, want %d", label, len(got), len(want))
|
|
}
|
|
for i := range got {
|
|
if got[i] != want[i] {
|
|
t.Fatalf("%s[%d] = %.17g, want %.17g", label, i, got[i], want[i])
|
|
}
|
|
}
|
|
}
|
|
|
|
// TestRollingFamilyAcceptsIntAndFloat32 pins the P0: rollingCheck
|
|
// returned RawFloats()[:n], which is a nil slice with capacity 0 for an
|
|
// int or a float32 payload, so all four reductions panicked. The int,
|
|
// float32 and float routes must agree value for value.
|
|
func TestRollingFamilyAcceptsIntAndFloat32(t *testing.T) {
|
|
vals := []float64{1, 4, 2, 8, 5, 7}
|
|
floats := pinVector(t, vals...)
|
|
ints := pinIntArray(t, 1, 4, 2, 8, 5, 7)
|
|
f32 := pinFloat32Array(t, 1, 4, 2, 8, 5, 7)
|
|
cases := []struct {
|
|
name string
|
|
call func(a *core.Array) (*core.Array, error)
|
|
}{
|
|
{"RollingMean", func(a *core.Array) (*core.Array, error) { return RollingMean(a, 3) }},
|
|
{"RollingSum", func(a *core.Array) (*core.Array, error) { return RollingSum(a, 2) }},
|
|
{"RollingMax", func(a *core.Array) (*core.Array, error) { return RollingMax(a, 3) }},
|
|
{"RollingMin", func(a *core.Array) (*core.Array, error) { return RollingMin(a, 2) }},
|
|
}
|
|
for _, tc := range cases {
|
|
t.Run(tc.name, func(t *testing.T) {
|
|
want, err := tc.call(floats)
|
|
if err != nil {
|
|
t.Fatalf("float route: %v", err)
|
|
}
|
|
wantVals := make([]float64, want.Len())
|
|
for i := range wantVals {
|
|
wantVals[i] = want.FloatAt(i)
|
|
}
|
|
for _, dt := range []struct {
|
|
name string
|
|
a *core.Array
|
|
}{{"int", ints}, {"float32", f32}} {
|
|
got, err := tc.call(dt.a)
|
|
if err != nil {
|
|
t.Fatalf("%s route: %v", dt.name, err)
|
|
}
|
|
gotVals := make([]float64, got.Len())
|
|
for i := range gotVals {
|
|
gotVals[i] = got.FloatAt(i)
|
|
}
|
|
pinSameFloats(t, dt.name, gotVals, wantVals)
|
|
}
|
|
})
|
|
}
|
|
}
|
|
|
|
// TestRollingExtremesAllNaNWindow pins the P2: the window fold seeded
|
|
// its running best with ±Inf and only replaced it on a strict
|
|
// comparison, so an all-NaN window answered -Inf/+Inf where core.Min
|
|
// and core.Max answer NaN. A NaN that shares a window with a number
|
|
// must still lose.
|
|
func TestRollingExtremesAllNaNWindow(t *testing.T) {
|
|
allNaN := pinVector(t, math.NaN(), math.NaN())
|
|
max, err := RollingMax(allNaN, 2)
|
|
if err != nil {
|
|
t.Fatalf("RollingMax: %v", err)
|
|
}
|
|
if !math.IsNaN(max.FloatAt(0)) {
|
|
t.Fatalf("RollingMax of an all-NaN window = %v, want NaN as core.Max gives", max.FloatAt(0))
|
|
}
|
|
min, err := RollingMin(allNaN, 2)
|
|
if err != nil {
|
|
t.Fatalf("RollingMin: %v", err)
|
|
}
|
|
if !math.IsNaN(min.FloatAt(0)) {
|
|
t.Fatalf("RollingMin of an all-NaN window = %v, want NaN as core.Min gives", min.FloatAt(0))
|
|
}
|
|
// A NaN never wins against a number.
|
|
mixed := pinVector(t, math.NaN(), 5, 2, math.NaN())
|
|
wantMax := []float64{5, 5, 2}
|
|
wantMin := []float64{5, 2, 2}
|
|
gotMax, err := RollingMax(mixed, 2)
|
|
if err != nil {
|
|
t.Fatalf("RollingMax: %v", err)
|
|
}
|
|
gotMin, err := RollingMin(mixed, 2)
|
|
if err != nil {
|
|
t.Fatalf("RollingMin: %v", err)
|
|
}
|
|
for i := range wantMax {
|
|
if gotMax.FloatAt(i) != wantMax[i] {
|
|
t.Fatalf("RollingMax[%d] = %v, want %v (NaN never wins)", i, gotMax.FloatAt(i), wantMax[i])
|
|
}
|
|
if gotMin.FloatAt(i) != wantMin[i] {
|
|
t.Fatalf("RollingMin[%d] = %v, want %v (NaN never wins)", i, gotMin.FloatAt(i), wantMin[i])
|
|
}
|
|
}
|
|
}
|
|
|
|
// TestQuantileRejectsNaNQuantile pins the P0: `q < 0 || q > 1` is
|
|
// false for NaN, so the quantile reached int(math.Floor(NaN)) and
|
|
// indexed the sorted values with it, panicking instead of refusing.
|
|
func TestQuantileRejectsNaNQuantile(t *testing.T) {
|
|
a := pinVector(t, 1, 2, 3, 4)
|
|
for _, qs := range [][]float64{{math.NaN()}, {0.5, math.NaN()}, {math.Inf(1)}, {math.Inf(-1)}} {
|
|
out, err := Quantile(a, qs)
|
|
if err == nil {
|
|
t.Fatalf("Quantile(%v) = %v, want a refusal", qs, out)
|
|
}
|
|
}
|
|
// The documented domain still works.
|
|
out, err := Quantile(a, []float64{0.25, 0.5, 0.75})
|
|
if err != nil {
|
|
t.Fatalf("Quantile: %v", err)
|
|
}
|
|
pinSameFloats(t, "Quantile", []float64{out.FloatAt(0), out.FloatAt(1), out.FloatAt(2)},
|
|
[]float64{1.75, 2.5, 3.25})
|
|
}
|
|
|
|
// TestHistogram2DRejectsComplexInput pins the P0: every other entry
|
|
// point of the package refuses complex input by name, while Histogram2D
|
|
// walked its float accessor into a nil int payload and panicked.
|
|
func TestHistogram2DRejectsComplexInput(t *testing.T) {
|
|
real := pinVector(t, 1, 2, 3)
|
|
cx, err := core.FromComplexes([]complex128{1, 2, 3}, 3)
|
|
if err != nil {
|
|
t.Fatalf("FromComplexes: %v", err)
|
|
}
|
|
if _, _, _, err := Histogram2D(cx, real, 2, 2); err == nil {
|
|
t.Fatal("complex x accepted")
|
|
}
|
|
if _, _, _, err := Histogram2D(real, cx, 2, 2); err == nil {
|
|
t.Fatal("complex y accepted")
|
|
}
|
|
// The real path is untouched.
|
|
counts, _, _, err := Histogram2D(real, real, 2, 2)
|
|
if err != nil {
|
|
t.Fatalf("Histogram2D: %v", err)
|
|
}
|
|
if counts.Shape()[0] != 2 || counts.Shape()[1] != 2 {
|
|
t.Fatalf("shape %v, want [2 2]", counts.Shape())
|
|
}
|
|
}
|
|
|
|
// TestHistogramBinCountsAreBounded pins the P1: neither Histogram nor
|
|
// Histogram2D bounded its bin count, so a hostile request either asked
|
|
// the allocator for an unbounded buffer or, with a product that wraps,
|
|
// produced an empty count payload the first write then indexed. The
|
|
// guard must fire before anything is allocated.
|
|
func TestHistogramBinCountsAreBounded(t *testing.T) {
|
|
a := pinVector(t, 0, 1)
|
|
if _, _, err := Histogram(a, maxHistBins+1); err == nil {
|
|
t.Fatal("a bin count past the limit was accepted")
|
|
}
|
|
if _, _, err := Histogram(a, 1<<21); err == nil {
|
|
t.Fatal("a 2^21-bin request was accepted")
|
|
}
|
|
counts, edges, err := Histogram(a, 4)
|
|
if err != nil {
|
|
t.Fatalf("a modest histogram: %v", err)
|
|
}
|
|
if counts.Len() != 4 || edges.Len() != 5 {
|
|
t.Fatalf("lengths %d and %d, want 4 and 5", counts.Len(), edges.Len())
|
|
}
|
|
if _, _, _, err := Histogram2D(a, a, maxHistBins+1, 1); err == nil {
|
|
t.Fatal("a per-axis bin count past the limit was accepted")
|
|
}
|
|
// 1025 * 1024 = 1049600 is past the total, each axis alone is not.
|
|
if _, _, _, err := Histogram2D(a, a, 1025, 1024); err == nil {
|
|
t.Fatal("a bin product past the limit was accepted")
|
|
}
|
|
// A wide but legal request goes through: the guard refuses what is
|
|
// past the limit, not what is merely large.
|
|
if _, _, _, err := Histogram2D(a, a, 512, 512); err != nil {
|
|
t.Fatalf("a 512 x 512 histogram: %v", err)
|
|
}
|
|
}
|
|
|
|
// TestNormalQuantileLowerTail pins the P1: continuousQuantile brackets
|
|
// upwards from the seed only, and Phi(0) = 0.5 is the infimum it can
|
|
// reach, so every q < 0.5 failed to bracket. The mirrored tail must
|
|
// reproduce an independent high-precision inverse of Phi.
|
|
func TestNormalQuantileLowerTail(t *testing.T) {
|
|
qs := []float64{1e-6, 1e-5, 1e-4, 1e-3, 0.01, 0.025, 0.05, 0.1, 0.25, 0.4, 0.4999}
|
|
worst := 0.0
|
|
for _, q := range qs {
|
|
z, err := NormalQuantile(q)
|
|
if err != nil {
|
|
t.Fatalf("NormalQuantile(%g): %v", q, err)
|
|
}
|
|
if z >= 0 {
|
|
t.Fatalf("NormalQuantile(%g) = %v, want a negative quantile", q, z)
|
|
}
|
|
ref := bigNormalQuantile(q)
|
|
dev := math.Abs(z - ref)
|
|
if dev > worst {
|
|
worst = dev
|
|
}
|
|
if dev > 1e-9 {
|
|
t.Fatalf("NormalQuantile(%g) = %.17g, the 256-bit bisection reference is %.17g (off by %.3g)",
|
|
q, z, ref, dev)
|
|
}
|
|
// The symmetry is exact in the implementation and must hold to
|
|
// the last bit: the mirrored call returns the negated value.
|
|
mirror, err := NormalQuantile(1 - q)
|
|
if err != nil {
|
|
t.Fatalf("NormalQuantile(%g): %v", 1-q, err)
|
|
}
|
|
if z+mirror != 0 {
|
|
t.Fatalf("NormalQuantile(%g) + NormalQuantile(%g) = %.17g, want exactly 0", q, 1-q, z+mirror)
|
|
}
|
|
// And the library's own CDF inverts the result.
|
|
if dev := math.Abs(NormalCDF(z) - q); dev > 1e-12 {
|
|
t.Fatalf("NormalCDF(NormalQuantile(%g)) is off by %.3g", q, dev)
|
|
}
|
|
}
|
|
t.Logf("worst deviation from the 256-bit reference: %.3g", worst)
|
|
// The upper half was already correct and must stay so.
|
|
upper, err := NormalQuantile(0.975)
|
|
if err != nil {
|
|
t.Fatalf("NormalQuantile(0.975): %v", err)
|
|
}
|
|
if math.Abs(upper-1.9599639845400532) > 1e-12 {
|
|
t.Fatalf("NormalQuantile(0.975) = %.17g, want 1.9599639845400532", upper)
|
|
}
|
|
}
|
|
|
|
// TestNormalQuantileNaNRejected pins the P2: NaN passed the written-out
|
|
// range guard and every comparison of the bisection returned false, so
|
|
// the seed came back as the answer (v = 1 with a nil error).
|
|
func TestNormalQuantileNaNRejected(t *testing.T) {
|
|
for _, q := range []float64{math.NaN(), math.Inf(1), math.Inf(-1), -0.5, 1.5, 0, 1} {
|
|
v, err := NormalQuantile(q)
|
|
if err == nil {
|
|
t.Fatalf("NormalQuantile(%v) = %v, want a refusal", q, v)
|
|
}
|
|
}
|
|
if v, err := StudentTQuantile(math.NaN(), 5); err == nil {
|
|
t.Fatalf("StudentTQuantile(NaN, 5) = %v, want a refusal", v)
|
|
}
|
|
}
|
|
|
|
// TestLogisticRegressionSaturatedLikelihoodFinite pins the P1: the
|
|
// convergence exit recomputed the sigmoid without the clamp the loop
|
|
// relies on, so one far-out covariate saturated Fitted to exactly 0/1
|
|
// and turned LogLikelihood into 0·log(0) = NaN while Converged was
|
|
// true. The reported values must be the ones the loop maximised.
|
|
func TestLogisticRegressionSaturatedLikelihoodFinite(t *testing.T) {
|
|
design, err := core.FromFloats([]float64{
|
|
1, 0, 1, 1, 1, 2, 1, 3, 1, 4, 1, 5, 1, 1e3,
|
|
}, 7, 2)
|
|
if err != nil {
|
|
t.Fatalf("FromFloats: %v", err)
|
|
}
|
|
y := pinVector(t, 0, 1, 0, 1, 1, 0, 1)
|
|
res, err := LogisticRegression(design, y)
|
|
if err != nil {
|
|
t.Fatalf("LogisticRegression: %v", err)
|
|
}
|
|
if !res.Converged {
|
|
t.Fatal("the fit did not converge")
|
|
}
|
|
manual := 0.0
|
|
for i, p := range res.Fitted {
|
|
if p <= 0 || p >= 1 {
|
|
t.Fatalf("Fitted[%d] = %v, want a probability strictly inside (0, 1)", i, p)
|
|
}
|
|
yv := y.FloatAt(i)
|
|
manual += yv*math.Log(p) + (1-yv)*math.Log(1-p)
|
|
}
|
|
if math.IsNaN(res.LogLikelihood) || math.IsInf(res.LogLikelihood, 0) {
|
|
t.Fatalf("LogLikelihood = %v, want the finite value the loop maximised", res.LogLikelihood)
|
|
}
|
|
if math.Abs(res.LogLikelihood-manual) > 1e-12 {
|
|
t.Fatalf("LogLikelihood = %.17g, recomputing from Fitted gives %.17g",
|
|
res.LogLikelihood, manual)
|
|
}
|
|
// The inference must stay finite as well.
|
|
for j := range res.Coefficients {
|
|
mustFinite := []float64{res.Coefficients[j], res.StandardErrors[j],
|
|
res.ZStatistics[j], res.PValues[j]}
|
|
for _, v := range mustFinite {
|
|
if math.IsNaN(v) || math.IsInf(v, 0) {
|
|
t.Fatalf("coefficient %d has a non-finite statistic: %v", j, mustFinite)
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// TestWeightedLinearRegressionInterceptStatistics pins the P1: the
|
|
// intercept column of the unweighted design stops being constant once
|
|
// the sqrt weights are applied, so the delegated fit reported the
|
|
// no-intercept conventions (uncentred TSS, DModel = p, an F test
|
|
// against the zero model) for a regression the caller described with an
|
|
// intercept. The statistics must be the weighted-theory ones, against
|
|
// the weighted mean.
|
|
func TestWeightedLinearRegressionInterceptStatistics(t *testing.T) {
|
|
design, err := core.FromFloats([]float64{1, 0, 1, 1, 1, 2, 1, 3, 1, 4, 1, 5}, 6, 2)
|
|
if err != nil {
|
|
t.Fatalf("FromFloats: %v", err)
|
|
}
|
|
yVals := []float64{10, 10.1, 9.9, 10, 10.2, 9.8}
|
|
y := pinVector(t, yVals...)
|
|
plain, err := LinearRegression(design, y)
|
|
if err != nil {
|
|
t.Fatalf("LinearRegression: %v", err)
|
|
}
|
|
// Uniform weights are the ordinary fit, statistic for statistic.
|
|
uniform, err := WeightedLinearRegression(design, y, pinVector(t, 1, 1, 1, 1, 1, 1))
|
|
if err != nil {
|
|
t.Fatalf("WeightedLinearRegression: %v", err)
|
|
}
|
|
if uniform.DModel != plain.DModel || uniform.RSquared != plain.RSquared ||
|
|
uniform.FStatistic != plain.FStatistic || uniform.FPValue != plain.FPValue {
|
|
t.Fatalf("uniform weights report DModel=%d R2=%.17g F=%.17g p=%.17g, the ordinary fit DModel=%d R2=%.17g F=%.17g p=%.17g",
|
|
uniform.DModel, uniform.RSquared, uniform.FStatistic, uniform.FPValue,
|
|
plain.DModel, plain.RSquared, plain.FStatistic, plain.FPValue)
|
|
}
|
|
for _, w6 := range []float64{1.0000001, 2} {
|
|
w := pinVector(t, 1, 1, 1, 1, 1, w6)
|
|
res, err := WeightedLinearRegression(design, y, w)
|
|
if err != nil {
|
|
t.Fatalf("WeightedLinearRegression(w6=%g): %v", w6, err)
|
|
}
|
|
// The reference from first principles, on the residuals the call
|
|
// itself reports: sumW, the weighted mean, the weighted centred
|
|
// total and the weighted residual sum.
|
|
sumW, sumWY, rss := 0.0, 0.0, 0.0
|
|
weights := []float64{1, 1, 1, 1, 1, w6}
|
|
for i, wi := range weights {
|
|
sumW += wi
|
|
sumWY += wi * yVals[i]
|
|
rss += wi * res.Residuals[i] * res.Residuals[i]
|
|
}
|
|
meanW := sumWY / sumW
|
|
tss := 0.0
|
|
for i, wi := range weights {
|
|
d := yVals[i] - meanW
|
|
tss += wi * d * d
|
|
}
|
|
wantR2 := 1 - rss/tss
|
|
wantF := (tss - rss) / 1 / (rss / 4)
|
|
if res.DModel != 1 {
|
|
t.Fatalf("w6=%g: DModel = %d, want 1 (the design carries the intercept)", w6, res.DModel)
|
|
}
|
|
if math.Abs(res.RSquared-wantR2) > 1e-12 {
|
|
t.Fatalf("w6=%g: R2 = %.17g, the weighted-centred reference is %.17g", w6, res.RSquared, wantR2)
|
|
}
|
|
if math.Abs(res.FStatistic-wantF) > 1e-9*wantF {
|
|
t.Fatalf("w6=%g: F = %.17g, the weighted-centred reference is %.17g", w6, res.FStatistic, wantF)
|
|
}
|
|
// For a single slope F = t^2, so the model p-value is the
|
|
// slope's own two-sided tail.
|
|
if math.Abs(res.FPValue-res.PValues[1]) > 1e-9 {
|
|
t.Fatalf("w6=%g: model p = %.17g, the slope's two-sided tail is %.17g",
|
|
w6, res.FPValue, res.PValues[1])
|
|
}
|
|
}
|
|
// The report's exact acceptance values on this data with w6 = 2, the
|
|
// weighted-centred values (R2 0.172938829787234, F 0.836401640003216,
|
|
// p 0.4121703772248073 from a numerical integration of F(1,4)):
|
|
// the defect reported 0.9998404595, 12534.0045019697 and
|
|
// 2.54531602909285e-08 instead.
|
|
res, err := WeightedLinearRegression(design, y, pinVector(t, 1, 1, 1, 1, 1, 2))
|
|
if err != nil {
|
|
t.Fatalf("WeightedLinearRegression: %v", err)
|
|
}
|
|
if math.Abs(res.RSquared-0.172938829787234) > 1e-9 {
|
|
t.Fatalf("R2 = %.17g, the weighted-centred value is 0.172938829787234", res.RSquared)
|
|
}
|
|
if math.Abs(res.FStatistic-0.836401640003216) > 1e-9 {
|
|
t.Fatalf("F = %.17g, the weighted-centred value is 0.836401640003216", res.FStatistic)
|
|
}
|
|
if math.Abs(res.FPValue-0.4121703772248073) > 1e-9 {
|
|
t.Fatalf("model p = %.17g, the weighted-centred value is 0.4121703772248073", res.FPValue)
|
|
}
|
|
if res.DModel != 1 || res.DResidual != 4 {
|
|
t.Fatalf("degrees of freedom (%d, %d), want (1, 4)", res.DModel, res.DResidual)
|
|
}
|
|
// A single weight a hair away from uniform must not move the
|
|
// statistics: the flip the defect produced was 0.05 to 0.9998.
|
|
hair, err := WeightedLinearRegression(design, y, pinVector(t, 1, 1, 1, 1, 1, 1.0000001))
|
|
if err != nil {
|
|
t.Fatalf("WeightedLinearRegression: %v", err)
|
|
}
|
|
if math.Abs(hair.RSquared-uniform.RSquared) > 1e-6 {
|
|
t.Fatalf("a 1e-7 weight change moved R2 from %.17g to %.17g", uniform.RSquared, hair.RSquared)
|
|
}
|
|
if math.Abs(hair.FPValue-uniform.FPValue) > 1e-6 {
|
|
t.Fatalf("a 1e-7 weight change moved the model p from %.17g to %.17g",
|
|
uniform.FPValue, hair.FPValue)
|
|
}
|
|
}
|
|
|
|
// TestBinomialDrawsRejectsNonPositiveCount pins the P2: the only draw
|
|
// function in the file without the n < 1 guard returned a nil array
|
|
// with a nil error, which a caller checking only the error dereferences.
|
|
func TestBinomialDrawsRejectsNonPositiveCount(t *testing.T) {
|
|
g := core.NewGenerator(1)
|
|
for _, n := range []int{-1, 0} {
|
|
out, err := BinomialDraws(g, n, 4, 0.5)
|
|
if err == nil {
|
|
t.Fatalf("BinomialDraws(n=%d) = %v, want a refusal", n, out)
|
|
}
|
|
if out != nil {
|
|
t.Fatalf("BinomialDraws(n=%d) returned a non-nil array with the refusal", n)
|
|
}
|
|
}
|
|
out, err := BinomialDraws(g, 3, 4, 0.5)
|
|
if err != nil {
|
|
t.Fatalf("BinomialDraws: %v", err)
|
|
}
|
|
if out.Len() != 3 {
|
|
t.Fatalf("length %d, want 3", out.Len())
|
|
}
|
|
}
|
|
|
|
// TestKolmogorovSmirnovRejectsNonFiniteSamples pins an unbounded loop
|
|
// found while verifying the package's finiteness policy: the merge walk
|
|
// in KolmogorovSmirnovTest advances only past values that compare true
|
|
// against the current support point, and a NaN never does, so the walk
|
|
// re-reads the same element for ever. A +Inf sample walks through but
|
|
// distorts the distance. The sibling tests (MannWhitneyU, ANOVAOneWay)
|
|
// refuse a non-finite sample by name and so must this one.
|
|
func TestKolmogorovSmirnovRejectsNonFiniteSamples(t *testing.T) {
|
|
clean := pinVector(t, 1, 2, 3, 4, 5)
|
|
nan := pinVector(t, 1, 2, math.NaN(), 4, 5)
|
|
inf := pinVector(t, 1, 2, math.Inf(1), 4, 5)
|
|
for _, tc := range []struct {
|
|
name string
|
|
a, b *core.Array
|
|
}{
|
|
{"NaN first sample", nan, clean},
|
|
{"NaN second sample", clean, nan},
|
|
{"+Inf first sample", inf, clean},
|
|
{"+Inf second sample", clean, inf},
|
|
} {
|
|
t.Run(tc.name, func(t *testing.T) {
|
|
var err error
|
|
pinWatchdog(t, 5*time.Second, "KolmogorovSmirnovTest("+tc.name+")", func() {
|
|
_, _, err = KolmogorovSmirnovTest(tc.a, tc.b)
|
|
})
|
|
if err == nil {
|
|
t.Fatal("a non-finite sample was accepted")
|
|
}
|
|
})
|
|
}
|
|
// The finite path is untouched: identical samples are distance zero
|
|
// with p = 1, and the report's {1,2,3} against {1.5,2.5,3.5} case
|
|
// gives d = 1/3 exactly.
|
|
d, p, err := KolmogorovSmirnovTest(clean, clean)
|
|
if err != nil {
|
|
t.Fatalf("KolmogorovSmirnovTest: %v", err)
|
|
}
|
|
if d != 0 || p != 1 {
|
|
t.Fatalf("identical samples give d = %g, p = %g, want 0 and 1", d, p)
|
|
}
|
|
d, _, err = KolmogorovSmirnovTest(pinVector(t, 1, 2, 3), pinVector(t, 1.5, 2.5, 3.5))
|
|
if err != nil {
|
|
t.Fatalf("KolmogorovSmirnovTest: %v", err)
|
|
}
|
|
// The largest gap is the first support point's 1/3, computed from
|
|
// integer counts: allow the couple of ulps that division costs.
|
|
if math.Abs(d-1.0/3.0) > 1e-15 {
|
|
t.Fatalf("d = %.17g, want 1/3", d)
|
|
}
|
|
}
|
|
|
|
// TestKernelDensityRejectsNonFiniteInput pins the first half of the
|
|
// finiteness gap: KernelDensity never scanned the sample or the
|
|
// evaluation points, so a NaN came back as a NaN estimate with a nil
|
|
// error and a +Inf silently dropped that sample from the average (the
|
|
// estimate was then the one of the remaining n-1 points). The clean
|
|
// path is pinned against the hand-written KDE definition.
|
|
func TestKernelDensityRejectsNonFiniteInput(t *testing.T) {
|
|
sample := pinVector(t, 0, 1, 2, 3)
|
|
points := pinVector(t, 0.5, 2.5)
|
|
for _, tc := range []struct {
|
|
name string
|
|
s, p *core.Array
|
|
}{
|
|
{"NaN sample", pinVector(t, 0, 1, math.NaN(), 3), points},
|
|
{"+Inf sample", pinVector(t, 0, 1, math.Inf(1), 3), points},
|
|
{"NaN point", sample, pinVector(t, 0.5, math.NaN())},
|
|
{"-Inf point", sample, pinVector(t, 0, math.Inf(-1))},
|
|
} {
|
|
t.Run(tc.name, func(t *testing.T) {
|
|
if out, err := KernelDensity(tc.s, 0.5, tc.p); err == nil {
|
|
t.Fatalf("a non-finite input was accepted, out = %v", out)
|
|
}
|
|
})
|
|
}
|
|
// The definition, written out: 1/(n h sqrt(2 pi)) * sum_k
|
|
// exp(-((x - sample[k])/h)^2 / 2) at x = 0.5 and x = 2.5, which are
|
|
// mirror images of each other in this sample and so agree exactly.
|
|
for _, x := range []float64{0.5, 2.5} {
|
|
want := 0.0
|
|
for _, s := range []float64{0, 1, 2, 3} {
|
|
z := (x - s) / 0.5
|
|
want += math.Exp(-0.5 * z * z)
|
|
}
|
|
want /= 4 * 0.5 * math.Sqrt(2*math.Pi)
|
|
out, err := KernelDensity(sample, 0.5, pinVector(t, x))
|
|
if err != nil {
|
|
t.Fatalf("KernelDensity(%g): %v", x, err)
|
|
}
|
|
if got := out.FloatAt(0); got != want {
|
|
t.Fatalf("KernelDensity at %g = %.17g, the hand-written KDE is %.17g", x, got, want)
|
|
}
|
|
}
|
|
}
|
|
|
|
// TestMultivariateNormalRejectsNonFiniteInput pins the second half of
|
|
// the finiteness gap: the covariance was scanned for non-finite
|
|
// entries, the mean and the point were not, so a NaN mean or point came
|
|
// back as a NaN log density and drew NaN vectors, all with a nil error.
|
|
func TestMultivariateNormalRejectsNonFiniteInput(t *testing.T) {
|
|
mean := pinVector(t, 0, 0)
|
|
cov, err := core.FromFloats([]float64{1, 0.5, 0.5, 1}, 2, 2)
|
|
if err != nil {
|
|
t.Fatalf("FromFloats: %v", err)
|
|
}
|
|
nanCov, err := core.FromFloats([]float64{1, math.NaN(), math.NaN(), 1}, 2, 2)
|
|
if err != nil {
|
|
t.Fatalf("FromFloats: %v", err)
|
|
}
|
|
x := pinVector(t, 1, 1)
|
|
nanMean := pinVector(t, math.NaN(), 0)
|
|
infMean := pinVector(t, math.Inf(1), 0)
|
|
nanPoint := pinVector(t, 1, math.NaN())
|
|
g := core.NewGenerator(7)
|
|
for _, tc := range []struct {
|
|
name string
|
|
call func() error
|
|
}{
|
|
{"density NaN mean", func() error { _, err := MultivariateNormalLogDensity(nanMean, cov, x); return err }},
|
|
{"density +Inf mean", func() error { _, err := MultivariateNormalLogDensity(infMean, cov, x); return err }},
|
|
{"density NaN covariance", func() error { _, err := MultivariateNormalLogDensity(mean, nanCov, x); return err }},
|
|
{"density NaN point", func() error { _, err := MultivariateNormalLogDensity(mean, cov, nanPoint); return err }},
|
|
{"draws NaN mean", func() error { _, err := MultivariateNormalDraws(g, 2, nanMean, cov); return err }},
|
|
{"draws NaN covariance", func() error { _, err := MultivariateNormalDraws(g, 2, mean, nanCov); return err }},
|
|
} {
|
|
t.Run(tc.name, func(t *testing.T) {
|
|
if err := tc.call(); err == nil {
|
|
t.Fatal("a non-finite input was accepted")
|
|
}
|
|
})
|
|
}
|
|
// The clean path is pinned by hand: with mean 0, cov [[1,.5],[.5,1]]
|
|
// and x = (1,1), the log density is -ln(2 pi) - 0.5 ln(0.75) - 2/3.
|
|
got, err := MultivariateNormalLogDensity(mean, cov, x)
|
|
if err != nil {
|
|
t.Fatalf("MultivariateNormalLogDensity: %v", err)
|
|
}
|
|
want := -math.Log(2*math.Pi) - 0.5*math.Log(0.75) - 2.0/3.0
|
|
if math.Abs(got-want) > 1e-12 {
|
|
t.Fatalf("log density = %.17g, hand value %.17g", got, want)
|
|
}
|
|
draws, err := MultivariateNormalDraws(g, 3, mean, cov)
|
|
if err != nil {
|
|
t.Fatalf("MultivariateNormalDraws: %v", err)
|
|
}
|
|
if draws.Len() != 6 {
|
|
t.Fatalf("draws length %d, want 6", draws.Len())
|
|
}
|
|
}
|
|
|
|
// TestMultivariateNormalRejectsAsymmetricCovariance pins the P2: the
|
|
// Cholesky factor read only the lower triangle, so an asymmetric matrix
|
|
// was silently symmetrised and the reported density described a
|
|
// different distribution from the one handed in. A mirror that differs
|
|
// beyond a relative tolerance is refused by name; a mirror that differs
|
|
// by rounding only is still accepted.
|
|
func TestMultivariateNormalRejectsAsymmetricCovariance(t *testing.T) {
|
|
mean := pinVector(t, 0, 0)
|
|
x := pinVector(t, 1, 1)
|
|
asym, err := core.FromFloats([]float64{1, 100, 0.5, 1}, 2, 2)
|
|
if err != nil {
|
|
t.Fatalf("FromFloats: %v", err)
|
|
}
|
|
for _, tc := range []struct {
|
|
name string
|
|
call func() error
|
|
}{
|
|
{"density", func() error { _, err := MultivariateNormalLogDensity(mean, asym, x); return err }},
|
|
{"draws", func() error { _, err := MultivariateNormalDraws(core.NewGenerator(3), 2, mean, asym); return err }},
|
|
} {
|
|
t.Run(tc.name, func(t *testing.T) {
|
|
if err := tc.call(); err == nil {
|
|
t.Fatal("the asymmetric covariance was accepted")
|
|
}
|
|
})
|
|
}
|
|
// Rounding-level asymmetry is not asymmetry: 0.5 + 1e-15 against 0.5
|
|
// is a mirror pair an A*A^T assembly can produce.
|
|
sym, err := core.FromFloats([]float64{1, 0.5 + 1e-15, 0.5, 1}, 2, 2)
|
|
if err != nil {
|
|
t.Fatalf("FromFloats: %v", err)
|
|
}
|
|
got, err := MultivariateNormalLogDensity(mean, sym, x)
|
|
if err != nil {
|
|
t.Fatalf("a rounding-level mirror difference was refused: %v", err)
|
|
}
|
|
want := -math.Log(2*math.Pi) - 0.5*math.Log(0.75) - 2.0/3.0
|
|
if math.Abs(got-want) > 1e-12 {
|
|
t.Fatalf("log density = %.17g, the symmetric hand value is %.17g", got, want)
|
|
}
|
|
}
|
|
|
|
// TestKolmogorovSeriesTruncationRule pins the inert defect: the early
|
|
// exit asked whether the current term was 1e-14 times the next, which is
|
|
// exp(2(2k+1)l^2) > 1 for every l > 0, so it could never fire and the
|
|
// comment described a termination check that was not there. The
|
|
// truncation is now the first term below the rounding level; the count
|
|
// must be exactly that term, and the summed series must agree with a
|
|
// directly summed one to the last bit.
|
|
func TestKolmogorovSeriesTruncationRule(t *testing.T) {
|
|
term := func(k int, lambda float64) float64 {
|
|
return 2 * math.Exp(-2*float64(k)*float64(k)*lambda*lambda)
|
|
}
|
|
for _, lambda := range []float64{0.2, 0.3, 0.5, 0.75, 1, 1.5, 2, 3, 5} {
|
|
n := kolmogorovTermCount(lambda)
|
|
if n < 1 || n > maxKolmogorovTerms {
|
|
t.Fatalf("lambda = %g: %d terms is outside [1, %d]", lambda, n, maxKolmogorovTerms)
|
|
}
|
|
if term(n, lambda) >= kolmogorovEps {
|
|
t.Fatalf("lambda = %g: term %d is %g, still above the %g rounding level: the series stops too early",
|
|
lambda, n, term(n, lambda), kolmogorovEps)
|
|
}
|
|
if n > 1 && term(n-1, lambda) < kolmogorovEps {
|
|
t.Fatalf("lambda = %g: term %d is already below the rounding level: the series adds terms that cannot move the sum",
|
|
lambda, n-1)
|
|
}
|
|
direct := 0.0
|
|
for k := 1; k <= 400; k++ {
|
|
if k%2 == 0 {
|
|
direct -= term(k, lambda)
|
|
} else {
|
|
direct += term(k, lambda)
|
|
}
|
|
}
|
|
want := min(1, max(0, direct))
|
|
if got := kolmogorovTail(lambda); got != want {
|
|
t.Fatalf("lambda = %g: tail = %.17g, a directly summed 400-term series gives %.17g", lambda, got, want)
|
|
}
|
|
}
|
|
// The documented short circuit below lambda = 0.2 is untouched.
|
|
if got := kolmogorovTail(0.1); got != 1 {
|
|
t.Fatalf("kolmogorovTail(0.1) = %g, want 1", got)
|
|
}
|
|
}
|
|
|
|
// bootstrapMean is the bootstrap statistic the sweeps below resample with.
|
|
func bootstrapMean(a *core.Array) (float64, error) {
|
|
return core.Mean(a)
|
|
}
|
|
|
|
// TestNaNParametersRejected pins the class sweep over the package's
|
|
// float parameters. Each guard used to be written out (`x < 0 || x > 1`
|
|
// and friends), which a NaN walks through because every comparison
|
|
// against a NaN is false: the parameter reached the algorithm and the
|
|
// call answered a value, a NaN, or nothing at all. The NaN-rejecting
|
|
// form `!(x >= 0 && x <= 1)` must turn every one of them into the
|
|
// refusal the caller deserves. The watchdog turns a regression to the
|
|
// unbounded discrete-quantile loop into a failure rather than a stuck
|
|
// test binary.
|
|
func TestNaNParametersRejected(t *testing.T) {
|
|
g := core.NewGenerator(1)
|
|
clean := pinVector(t, 1, 2, 3, 4, 5)
|
|
nanSample := pinVector(t, 1, 2, math.NaN(), 4, 5)
|
|
cases := []struct {
|
|
name string
|
|
call func() error
|
|
}{
|
|
{"GammaLower(NaN shape)", func() error { _, err := GammaLower(math.NaN(), 1); return err }},
|
|
{"GammaLower(NaN x)", func() error { _, err := GammaLower(1, math.NaN()); return err }},
|
|
{"GammaUpper(NaN shape)", func() error { _, err := GammaUpper(math.NaN(), 1); return err }},
|
|
{"GammaUpper(NaN x)", func() error { _, err := GammaUpper(1, math.NaN()); return err }},
|
|
{"BetaIncomplete(NaN x)", func() error { _, err := BetaIncomplete(math.NaN(), 2, 2); return err }},
|
|
{"BetaIncomplete(NaN a)", func() error { _, err := BetaIncomplete(0.5, math.NaN(), 2); return err }},
|
|
{"BetaIncomplete(NaN b)", func() error { _, err := BetaIncomplete(0.5, 2, math.NaN()); return err }},
|
|
{"ExponentialCDF(NaN x)", func() error { _, err := ExponentialCDF(math.NaN(), 2); return err }},
|
|
{"ExponentialCDF(NaN rate)", func() error { _, err := ExponentialCDF(1, math.NaN()); return err }},
|
|
{"GammaCDF(NaN x)", func() error { _, err := GammaCDF(math.NaN(), 2, 3); return err }},
|
|
{"GammaCDF(NaN shape)", func() error { _, err := GammaCDF(1, math.NaN(), 3); return err }},
|
|
{"GammaCDF(NaN rate)", func() error { _, err := GammaCDF(1, 2, math.NaN()); return err }},
|
|
{"ChiSquareCDF(NaN x)", func() error { _, err := ChiSquareCDF(math.NaN(), 2); return err }},
|
|
{"StudentTCDF(NaN t)", func() error { _, err := StudentTCDF(math.NaN(), 5); return err }},
|
|
{"PoissonCDF(NaN lambda)", func() error { _, err := PoissonCDF(3, math.NaN()); return err }},
|
|
{"BinomialCDF(NaN p)", func() error { _, err := BinomialCDF(3, 10, math.NaN()); return err }},
|
|
{"ExponentialQuantile(NaN q)", func() error { _, err := ExponentialQuantile(math.NaN(), 2); return err }},
|
|
{"ExponentialQuantile(NaN rate)", func() error { _, err := ExponentialQuantile(0.5, math.NaN()); return err }},
|
|
{"GammaQuantile(NaN q)", func() error { _, err := GammaQuantile(math.NaN(), 2, 3); return err }},
|
|
{"GammaQuantile(NaN shape)", func() error { _, err := GammaQuantile(0.5, math.NaN(), 3); return err }},
|
|
{"GammaQuantile(NaN rate)", func() error { _, err := GammaQuantile(0.5, 2, math.NaN()); return err }},
|
|
{"ChiSquareQuantile(NaN q)", func() error { _, err := ChiSquareQuantile(math.NaN(), 2); return err }},
|
|
{"StudentTQuantile(NaN q)", func() error { _, err := StudentTQuantile(math.NaN(), 5); return err }},
|
|
{"PoissonQuantile(NaN q)", func() error { _, err := PoissonQuantile(math.NaN(), 3); return err }},
|
|
{"PoissonQuantile(NaN lambda)", func() error { _, err := PoissonQuantile(0.5, math.NaN()); return err }},
|
|
{"BinomialQuantile(NaN q)", func() error { _, err := BinomialQuantile(math.NaN(), 0.5, 5); return err }},
|
|
{"BinomialQuantile(NaN p)", func() error { _, err := BinomialQuantile(0.5, math.NaN(), 5); return err }},
|
|
{"ExponentialDraws(NaN rate)", func() error { _, err := ExponentialDraws(g, 3, math.NaN()); return err }},
|
|
{"GammaDraws(NaN shape)", func() error { _, err := GammaDraws(g, 3, math.NaN(), 2); return err }},
|
|
{"GammaDraws(NaN rate)", func() error { _, err := GammaDraws(g, 3, 2, math.NaN()); return err }},
|
|
{"PoissonDraws(NaN lambda)", func() error { _, err := PoissonDraws(g, 3, math.NaN()); return err }},
|
|
{"BinomialDraws(NaN p)", func() error { _, err := BinomialDraws(g, 3, 4, math.NaN()); return err }},
|
|
{"ChiSquareGoodnessOfFit(NaN observed)", func() error {
|
|
_, _, _, err := ChiSquareGoodnessOfFit(nanSample, clean)
|
|
return err
|
|
}},
|
|
{"ChiSquareGoodnessOfFit(NaN expected)", func() error {
|
|
_, _, _, err := ChiSquareGoodnessOfFit(clean, nanSample)
|
|
return err
|
|
}},
|
|
{"BootstrapCI(NaN level)", func() error { _, _, err := BootstrapCI(clean, bootstrapMean, math.NaN(), 20, 1); return err }},
|
|
{"TrimmedMean(NaN fraction)", func() error { _, err := TrimmedMean(clean, math.NaN()); return err }},
|
|
{"KernelDensity(NaN bandwidth)", func() error { _, err := KernelDensity(clean, math.NaN(), clean); return err }},
|
|
{"WelchTTest(NaN sample)", func() error { _, _, _, err := WelchTTest(nanSample, clean); return err }},
|
|
}
|
|
for _, tc := range cases {
|
|
t.Run(tc.name, func(t *testing.T) {
|
|
var err error
|
|
pinWatchdog(t, 5*time.Second, tc.name, func() { err = tc.call() })
|
|
if err == nil {
|
|
t.Fatalf("%s answered with no error for a NaN parameter", tc.name)
|
|
}
|
|
})
|
|
}
|
|
// No guard overreaches: the legal domain still answers, and the two
|
|
// values below are hand-checkable. The incomplete beta is the
|
|
// continued fraction's own rounding level, a few ulps at 0.5.
|
|
if p, err := BetaIncomplete(0.5, 2, 2); err != nil || math.Abs(p-0.5) > 1e-15 {
|
|
t.Fatalf("BetaIncomplete(0.5, 2, 2) = %v, %v, want 0.5 by symmetry", p, err)
|
|
}
|
|
if p, err := GammaLower(1, 1); err != nil || math.Abs(p-(1-math.Exp(-1))) > 1e-15 {
|
|
t.Fatalf("GammaLower(1, 1) = %v, %v, want 1 - e^-1", p, err)
|
|
}
|
|
if q, err := ExponentialQuantile(0.5, 3); err != nil || math.Abs(q-math.Ln2/3) > 1e-15 {
|
|
t.Fatalf("ExponentialQuantile(0.5, 3) = %v, %v, want ln 2 / 3", q, err)
|
|
}
|
|
if p, err := PoissonQuantile(0.5, 3); err != nil || p != 3 {
|
|
t.Fatalf("PoissonQuantile(0.5, 3) = %v, %v, want 3", p, err)
|
|
}
|
|
}
|
|
|
|
// TestComplexInputsRejected pins the dtype sweep: every entry point in
|
|
// the package refuses a complex array by name, and the ones exercised
|
|
// here (the three MVN inputs, BootstrapCI's data, KernelDensity's
|
|
// sample and points) walked the float accessor into a nil payload and
|
|
// panicked instead.
|
|
func TestComplexInputsRejected(t *testing.T) {
|
|
cx2, err := core.FromComplexes([]complex128{1 + 1i, 2 - 1i}, 2)
|
|
if err != nil {
|
|
t.Fatalf("FromComplexes: %v", err)
|
|
}
|
|
cx5, err := core.FromComplexes([]complex128{1 + 1i, 2, 3, 4, 5}, 5)
|
|
if err != nil {
|
|
t.Fatalf("FromComplexes: %v", err)
|
|
}
|
|
cxMat, err := core.FromComplexes([]complex128{1, 0, 0, 1}, 2, 2)
|
|
if err != nil {
|
|
t.Fatalf("FromComplexes: %v", err)
|
|
}
|
|
mean := pinVector(t, 0, 0)
|
|
x := pinVector(t, 1, 1)
|
|
cov, err := core.FromFloats([]float64{1, 0.5, 0.5, 1}, 2, 2)
|
|
if err != nil {
|
|
t.Fatalf("FromFloats: %v", err)
|
|
}
|
|
g := core.NewGenerator(3)
|
|
for _, tc := range []struct {
|
|
name string
|
|
call func() error
|
|
}{
|
|
{"MVN density complex mean", func() error { _, err := MultivariateNormalLogDensity(cx2, cov, x); return err }},
|
|
{"MVN density complex covariance", func() error { _, err := MultivariateNormalLogDensity(mean, cxMat, x); return err }},
|
|
{"MVN density complex point", func() error { _, err := MultivariateNormalLogDensity(mean, cov, cx2); return err }},
|
|
{"MVN draws complex mean", func() error { _, err := MultivariateNormalDraws(g, 2, cx2, cov); return err }},
|
|
{"MVN draws complex covariance", func() error { _, err := MultivariateNormalDraws(g, 2, mean, cxMat); return err }},
|
|
{"BootstrapCI complex data", func() error { _, _, err := BootstrapCI(cx5, bootstrapMean, 0.95, 20, 1); return err }},
|
|
{"KernelDensity complex sample", func() error { _, err := KernelDensity(cx5, 0.5, cx5); return err }},
|
|
{"KernelDensity complex points", func() error { _, err := KernelDensity(pinVector(t, 1, 2, 3, 4, 5), 0.5, cx5); return err }},
|
|
{"Histogram2D complex sample", func() error { _, _, _, err := Histogram2D(cx5, cx5, 2, 2); return err }},
|
|
} {
|
|
t.Run(tc.name, func(t *testing.T) {
|
|
if err := tc.call(); err == nil {
|
|
t.Fatal("a complex input was accepted")
|
|
}
|
|
})
|
|
}
|
|
// The real path is untouched: a constant sample has a degenerate
|
|
// bootstrap interval at its own value.
|
|
lo, hi, err := BootstrapCI(pinVector(t, 5, 5, 5, 5), bootstrapMean, 0.95, 50, 1)
|
|
if err != nil {
|
|
t.Fatalf("BootstrapCI: %v", err)
|
|
}
|
|
if lo != 5 || hi != 5 {
|
|
t.Fatalf("the constant sample gives [%g, %g], want [5, 5]", lo, hi)
|
|
}
|
|
}
|
|
|
|
// TestLinearRegressionConstantResponseStatistics pins the NaN guard a
|
|
// degenerate response needs: with no variation in y the fit reproduces
|
|
// it exactly, so the explained and residual sums are both zero and the
|
|
// F statistic is 0/0. A NaN there would spread through every consumer
|
|
// of the result, where the finite zero is the value the tail reads as
|
|
// p = 1: there is no evidence of a model.
|
|
func TestLinearRegressionConstantResponseStatistics(t *testing.T) {
|
|
design, err := core.FromFloats([]float64{1, 0, 1, 1, 1, 2, 1, 3}, 4, 2)
|
|
if err != nil {
|
|
t.Fatalf("FromFloats: %v", err)
|
|
}
|
|
y := pinVector(t, 5, 5, 5, 5)
|
|
for _, tc := range []struct {
|
|
name string
|
|
fit func() (*LinearRegressionResult, error)
|
|
}{
|
|
{"ordinary", func() (*LinearRegressionResult, error) { return LinearRegression(design, y) }},
|
|
{"weighted", func() (*LinearRegressionResult, error) {
|
|
return WeightedLinearRegression(design, y, pinVector(t, 1, 1, 1, 2))
|
|
}},
|
|
} {
|
|
t.Run(tc.name, func(t *testing.T) {
|
|
res, err := tc.fit()
|
|
if err != nil {
|
|
t.Fatalf("the constant-response fit: %v", err)
|
|
}
|
|
if math.IsNaN(res.FStatistic) || res.FStatistic != 0 {
|
|
t.Fatalf("FStatistic = %v, want the finite 0", res.FStatistic)
|
|
}
|
|
if res.FPValue != 1 {
|
|
t.Fatalf("FPValue = %v, want 1", res.FPValue)
|
|
}
|
|
if res.DModel != 1 || res.DResidual != 2 {
|
|
t.Fatalf("degrees of freedom (%d, %d), want (1, 2)", res.DModel, res.DResidual)
|
|
}
|
|
})
|
|
}
|
|
}
|
|
|
|
// bigPhi returns Phi(x) in 256-bit arithmetic through the Taylor
|
|
// series of the error function at x/sqrt(2). Nothing here calls the
|
|
// library, so it is an independent reference for the quantile tests.
|
|
func bigPhi(x *big.Float) *big.Float {
|
|
const prec = 256
|
|
one := new(big.Float).SetPrec(prec).SetInt64(1)
|
|
y := new(big.Float).SetPrec(prec).Quo(new(big.Float).SetPrec(prec).Set(x), bigSqrt2())
|
|
neg := y.Sign() < 0
|
|
ay := new(big.Float).SetPrec(prec).Abs(y)
|
|
// erf(ay) = (2/sqrt(pi)) * sum_{n>=0} (-1)^n ay^(2n+1) / (n! (2n+1)).
|
|
term := new(big.Float).SetPrec(prec).Set(ay) // ay^(2n+1)/n! at n = 0
|
|
sum := new(big.Float).SetPrec(prec).Set(ay)
|
|
ay2 := new(big.Float).SetPrec(prec).Mul(ay, ay)
|
|
for n := 1; ; n++ {
|
|
term.Mul(term, ay2)
|
|
term.Quo(term, new(big.Float).SetPrec(prec).SetInt64(int64(n)))
|
|
term.Neg(term)
|
|
contrib := new(big.Float).SetPrec(prec).Quo(term,
|
|
new(big.Float).SetPrec(prec).SetInt64(int64(2*n+1)))
|
|
sum.Add(sum, contrib)
|
|
if contrib.Sign() == 0 || contrib.MantExp(nil) < -400 {
|
|
break
|
|
}
|
|
}
|
|
// 2/sqrt(pi).
|
|
twoOverSqrtPi := new(big.Float).SetPrec(prec).Quo(
|
|
new(big.Float).SetPrec(prec).SetInt64(2),
|
|
new(big.Float).SetPrec(prec).Sqrt(new(big.Float).SetPrec(prec).SetFloat64(math.Pi)))
|
|
erf := new(big.Float).SetPrec(prec).Mul(sum, twoOverSqrtPi)
|
|
if neg {
|
|
erf.Neg(erf)
|
|
}
|
|
// Phi(x) = (1 + erf(x/sqrt(2)))/2.
|
|
return new(big.Float).SetPrec(prec).Quo(
|
|
new(big.Float).SetPrec(prec).Add(one, erf),
|
|
new(big.Float).SetPrec(prec).SetInt64(2))
|
|
}
|
|
|
|
// bigNormalQuantile inverts bigPhi by bisection on [-12, 12], the
|
|
// independent reference the NormalQuantile test compares against.
|
|
func bigNormalQuantile(q float64) float64 {
|
|
const prec = 256
|
|
target := new(big.Float).SetPrec(prec).SetFloat64(q)
|
|
lo := new(big.Float).SetPrec(prec).SetFloat64(-12)
|
|
hi := new(big.Float).SetPrec(prec).SetFloat64(12)
|
|
mid := new(big.Float).SetPrec(prec)
|
|
for range 220 {
|
|
mid.Add(lo, hi)
|
|
mid.Quo(mid, new(big.Float).SetPrec(prec).SetInt64(2))
|
|
if bigPhi(mid).Cmp(target) < 0 {
|
|
lo.Set(mid)
|
|
} else {
|
|
hi.Set(mid)
|
|
}
|
|
}
|
|
mid.Add(lo, hi)
|
|
mid.Quo(mid, new(big.Float).SetPrec(prec).SetInt64(2))
|
|
out, _ := mid.Float64()
|
|
return out
|
|
}
|