Files
tensor/stats/distrib2_test.go
petrbalvin af4ee19703
Release / gates (push) Successful in 4m38s
Test / test (push) Successful in 5m16s
Release / release (push) Successful in 35s
feat: initial release
Assisted-by: GLM 5.3 Flash
2026-09-03 10:00:00 +02:00

437 lines
18 KiB
Go
Raw Permalink Blame History

This file contains ambiguous Unicode characters
This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.
// Copyright (c) 2026 Petr Balvín <opensource@petrbalvin.org> (https://petrbalvin.org)
// SPDX-License-Identifier: MIT
package stats
import (
"math"
"strings"
"testing"
"sourcedock.dev/petrbalvin/tensor/internal/core"
)
// TestWeibullClosed pins the Weibull on its closed forms: with k = 2
// and λ = 1 the density is 2xe^{−x²}, the CDF 1 − e^{−x²} and the
// median √(ln 2); with k = 1 the law is the exponential of rate 1/λ.
func TestWeibullClosed(t *testing.T) {
d, err := WeibullDensity(1, 2, 1)
if err != nil || math.Abs(d-2/math.E) > 1e-15 {
t.Fatalf("WeibullDensity(1, 2, 1) = %v, %v, want 2/e", d, err)
}
v, err := WeibullCDF(1, 2, 1)
if err != nil || math.Abs(v-(1-1/math.E)) > 1e-15 {
t.Fatalf("WeibullCDF(1, 2, 1) = %v, %v, want 1 − 1/e", v, err)
}
v, err = WeibullCDF(3, 2, 1)
if err != nil || math.Abs(v-(1-math.Exp(-9))) > 1e-15 {
t.Fatalf("WeibullCDF(3, 2, 1) = %v, %v, want 1 − e⁻⁹", v, err)
}
q, err := WeibullQuantile(0.5, 2, 1)
if err != nil || math.Abs(q-math.Sqrt(math.Ln2)) > 1e-14 {
t.Fatalf("WeibullQuantile(0.5, 2, 1) = %v, %v, want √(ln 2)", q, err)
}
// k = 1 reduces to the exponential of rate 1/λ.
v, err = WeibullCDF(2, 1, 0.5)
exp, eerr := ExponentialCDF(2, 2)
if err != nil || eerr != nil || v != exp {
t.Fatalf("WeibullCDF(2, 1, 0.5) = %v vs ExponentialCDF %v (%v, %v)", v, exp, err, eerr)
}
// Support convention and round trips.
if v, _ = WeibullCDF(-1, 2, 1); v != 0 {
t.Fatalf("WeibullCDF below the support = %v, want 0", v)
}
if d, _ = WeibullDensity(-1, 2, 1); d != 0 {
t.Fatalf("WeibullDensity below the support = %v, want 0", d)
}
for _, q := range []float64{0.001, 0.1, 0.5, 0.9, 0.999} {
x, err := WeibullQuantile(q, 1.5, 2)
if err != nil {
t.Fatalf("WeibullQuantile(%g): %v", q, err)
}
back, err := WeibullCDF(x, 1.5, 2)
if err != nil || math.Abs(back-q) > 1e-13 {
t.Fatalf("round trip q = %g: CDF(quantile) = %v, %v", q, back, err)
}
}
}
// TestLognormalClosed pins the lognormal through the normal: the CDF
// at 1 with μ = 0 is exactly ½, the density there is 1/√(2π), and the
// 0.975 quantile is e to the normal quantile.
func TestLognormalClosed(t *testing.T) {
v, err := LognormalCDF(1, 0, 1)
if err != nil || v != 0.5 {
t.Fatalf("LognormalCDF(1, 0, 1) = %v, %v, want exactly 0.5", v, err)
}
v, err = LognormalCDF(math.E, 0, 1)
if err != nil || math.Abs(v-NormalCDF(1)) > 1e-15 {
t.Fatalf("LognormalCDF(e, 0, 1) = %v, %v, want Φ(1)", v, err)
}
d, err := LognormalDensity(1, 0, 1)
if want := 1 / (math.Sqrt2 * math.SqrtPi); err != nil || math.Abs(d-want) > 1e-15 {
t.Fatalf("LognormalDensity(1, 0, 1) = %v, %v, want %.16g", d, err, want)
}
z, _ := NormalQuantile(0.975)
q, err := LognormalQuantile(0.975, 0, 1)
if err != nil || math.Abs(q-math.Exp(z)) > 1e-12*math.Exp(z) {
t.Fatalf("LognormalQuantile(0.975, 0, 1) = %v, %v, want e^{%.16g}", q, err, z)
}
// The CDF is the normal CDF at ln x, pointwise.
for _, x := range []float64{0.05, 0.5, 2, 20} {
got, _ := LognormalCDF(x, 0.3, 0.8)
want := NormalCDF((math.Log(x) - 0.3) / 0.8)
if math.Abs(got-want) > 1e-15 {
t.Fatalf("LognormalCDF(%g, 0.3, 0.8) = %v, want %v", x, got, want)
}
}
if v, _ = LognormalCDF(0, 0, 1); v != 0 {
t.Fatalf("LognormalCDF at 0 = %v, want 0", v)
}
if d, _ = LognormalDensity(0, 0, 1); d != 0 {
t.Fatalf("LognormalDensity at 0 = %v, want 0", d)
}
// A subnormal x underflows the closed form's numerator and
// denominator together: the log-space tail must answer 0, not the
// NaN their division would report.
for _, sigma := range []float64{0.1, 0.01, 1} {
if d, err = LognormalDensity(5e-324, 0, sigma); err != nil || d != 0 {
t.Fatalf("LognormalDensity(5e-324, 0, %g) = %v, %v, want 0", sigma, d, err)
}
}
}
// TestParetoClosed pins the Pareto on its closed forms with x_m = 1
// and α = 3: the CDF at 2 is 1 − 2^{−3} = 7/8, the density there
// 3/16, and the 7/8 quantile exactly 2.
func TestParetoClosed(t *testing.T) {
v, err := ParetoCDF(2, 1, 3)
if err != nil || math.Abs(v-0.875) > 1e-15 {
t.Fatalf("ParetoCDF(2, 1, 3) = %v, %v, want 0.875", v, err)
}
d, err := ParetoDensity(2, 1, 3)
if err != nil || math.Abs(d-3.0/16) > 1e-15 {
t.Fatalf("ParetoDensity(2, 1, 3) = %v, %v, want 3/16", d, err)
}
q, err := ParetoQuantile(0.875, 1, 3)
if err != nil || math.Abs(q-2) > 1e-12 {
t.Fatalf("ParetoQuantile(0.875, 1, 3) = %v, %v, want 2", q, err)
}
if v, _ = ParetoCDF(0.5, 1, 3); v != 0 {
t.Fatalf("ParetoCDF below the support = %v, want 0", v)
}
if d, _ = ParetoDensity(0.5, 1, 3); d != 0 {
t.Fatalf("ParetoDensity below the support = %v, want 0", d)
}
// The Expm1 form keeps the digits just above the support. The
// exact probability of the representable input is 3u for the
// offset u the double actually carries.
input := 1 + 1e-12
u := input - 1
v, err = ParetoCDF(input, 1, 3)
if err != nil || math.Abs(v-3*u) > 1e-10*3*u {
t.Fatalf("ParetoCDF just above the support = %v, %v, want ≈ %.17g", v, err, 3*u)
}
}
// TestNegativeBinomialClosed pins the negative binomial on exact
// fractions with r = 3 and p = ½, and holds the summation CDF against
// the exact identity P(X ≤ k) = I_p(r, k+1).
func TestNegativeBinomialClosed(t *testing.T) {
d, err := NegativeBinomialPMF(0, 3, 0.5)
if err != nil || math.Abs(d-math.Pow(0.5, 3)) > 1e-16 {
t.Fatalf("NegativeBinomialPMF(0, 3, 0.5) = %v, %v, want 1/8", d, err)
}
d, err = NegativeBinomialPMF(1, 3, 0.5)
if err != nil || math.Abs(d-3*math.Pow(0.5, 4)) > 1e-16 {
t.Fatalf("NegativeBinomialPMF(1, 3, 0.5) = %v, %v, want 3/16", d, err)
}
v, err := NegativeBinomialCDF(1, 3, 0.5)
if err != nil || math.Abs(v-0.3125) > 1e-15 {
t.Fatalf("NegativeBinomialCDF(1, 3, 0.5) = %v, %v, want 5/16", v, err)
}
// Exact identity against the regularised beta.
for _, r := range []int{1, 2, 5, 12} {
for _, p := range []float64{0.2, 0.5, 0.8} {
for k := 0; k <= 25; k++ {
got, err := NegativeBinomialCDF(k, r, p)
if err != nil {
t.Fatalf("NegativeBinomialCDF(%d, %d, %g): %v", k, r, p, err)
}
want, err := BetaIncomplete(p, float64(r), float64(k+1))
if err != nil {
t.Fatalf("BetaIncomplete: %v", err)
}
if math.Abs(got-want) > 1e-14 {
t.Fatalf("NegativeBinomialCDF(%d, %d, %g) = %.16g, want I_p identity %.16g",
k, r, p, got, want)
}
}
}
}
// The CDF at k = 2 is exactly ½, so both sides of the smallest-k
// rule are pinned.
q, err := NegativeBinomialQuantile(0.5, 0.5, 3)
if err != nil || q != 2 {
t.Fatalf("NegativeBinomialQuantile(0.5, 0.5, 3) = %v, %v, want 2", q, err)
}
q, err = NegativeBinomialQuantile(0.3125, 0.5, 3)
if err != nil || q != 1 {
t.Fatalf("NegativeBinomialQuantile(0.3125, 0.5, 3) = %v, %v, want 1", q, err)
}
// q = 1 brackets once the summed tail underflows below a half ulp
// of 1, at a few hundred failures for these parameters.
qEnd, err := NegativeBinomialQuantile(1, 0.5, 3)
if err != nil || qEnd < 10 {
t.Fatalf("NegativeBinomialQuantile(1, 0.5, 3) = %v, %v", qEnd, err)
}
if d, err = NegativeBinomialPMF(-1, 3, 0.5); err != nil || d != 0 {
t.Fatalf("NegativeBinomialPMF below the support = %v, %v", d, err)
}
}
// TestDirichletDensityMoments pins the density on an exact rational
// value, the boundary conventions, and the mean and mode helpers.
// With α = (2, 3, 4) the normalising constant is
// B(α) = Γ2Γ3Γ4/Γ9 = 12/40320 = 1/3360, so the uniform point carries
// 3360·3^{−6} = 3360/729.
func TestDirichletDensityMoments(t *testing.T) {
third := 1.0 / 3
d, err := DirichletDensity([]float64{2, 3, 4}, []float64{third, third, third})
if err != nil || math.Abs(d-3360.0/729) > 1e-12*3360.0/729 {
t.Fatalf("DirichletDensity uniform = %v, %v, want %.16g", d, err, 3360.0/729)
}
// Boundary: α_i > 1 kills the density, α_i < 1 blows it up, and an
// α_i of exactly 1 contributes nothing.
if d, _ = DirichletDensity([]float64{2, 2}, []float64{0, 1}); d != 0 {
t.Fatalf("boundary density with α > 1 = %v, want 0", d)
}
if d, _ = DirichletDensity([]float64{0.5, 0.5}, []float64{0, 1}); !math.IsInf(d, 1) {
t.Fatalf("boundary density with α < 1 = %v, want +Inf", d)
}
d, err = DirichletDensity([]float64{1, 3}, []float64{0, 1})
if err != nil || math.Abs(d-3) > 1e-14 {
t.Fatalf("boundary density with α = 1 = %v, %v, want 3", d, err)
}
mean, err := DirichletMean([]float64{2, 3, 4})
if err != nil {
t.Fatalf("DirichletMean: %v", err)
}
for i, want := range []float64{2.0 / 9, 1.0 / 3, 4.0 / 9} {
if math.Abs(mean[i]-want) > 1e-15 {
t.Fatalf("DirichletMean[%d] = %.16g, want %.16g", i, mean[i], want)
}
}
mode, err := DirichletMode([]float64{2, 3, 4})
if err != nil {
t.Fatalf("DirichletMode: %v", err)
}
for i, want := range []float64{1.0 / 6, 1.0 / 3, 0.5} {
if math.Abs(mode[i]-want) > 1e-15 {
t.Fatalf("DirichletMode[%d] = %.16g, want %.16g", i, mode[i], want)
}
}
if _, err := DirichletMode([]float64{1, 2}); err == nil {
t.Fatal("an α of 1 leaves no interior mode: want an error")
}
}
// TestDirichletDraws checks the sampler statistically: every row is a
// probability vector, the column means track α/α₀, and the seed makes
// the run reproducible.
func TestDirichletDraws(t *testing.T) {
g := core.NewGenerator(7)
alpha := []float64{2, 3, 4}
const n = 200000
draws, err := DirichletDraws(g, n, alpha)
if err != nil {
t.Fatalf("DirichletDraws: %v", err)
}
if draws.Shape()[0] != n || draws.Shape()[1] != len(alpha) {
t.Fatalf("draw shape %v, want (%d, %d)", draws.Shape(), n, len(alpha))
}
colSum := make([]float64, len(alpha))
for r := range n {
sum := 0.0
for c := range alpha {
v := draws.FloatAt(r*len(alpha) + c)
if v < 0 {
t.Fatalf("draw (%d, %d) = %v is negative", r, c, v)
}
sum += v
colSum[c] += v
}
if math.Abs(sum-1) > 1e-12 {
t.Fatalf("row %d sums to %.17g, want 1", r, sum)
}
}
mean, _ := DirichletMean(alpha)
for c := range alpha {
m := colSum[c] / n
if math.Abs(m-mean[c]) > 0.005 {
t.Fatalf("column %d mean = %g, want ≈ %g", c, m, mean[c])
}
}
again, _ := DirichletDraws(core.NewGenerator(7), n, alpha)
for i := range n * len(alpha) {
if again.FloatAt(i) != draws.FloatAt(i) {
t.Fatalf("seeded run not deterministic at %d", i)
}
}
// Concentrations below 1 take the boosted gamma branch of the
// sampler; the mean still tracks α/α₀.
sparse, err := DirichletDraws(core.NewGenerator(3), n, []float64{0.5, 0.7})
if err != nil {
t.Fatalf("DirichletDraws: %v", err)
}
sums := []float64{0, 0}
for r := range n {
for c := range 2 {
v := sparse.FloatAt(r*2 + c)
sums[c] += v
}
}
for c, a := range []float64{0.5, 0.7} {
if m := sums[c] / n; math.Abs(m-a/1.2) > 0.005 {
t.Fatalf("sparse column %d mean = %g, want ≈ %g", c, m, a/1.2)
}
}
if _, err := DirichletDraws(core.NewGenerator(1), 1, []float64{2, -1}); err == nil {
t.Fatal("a negative concentration: want an error")
}
}
// TestDistribution2Errors pins the parameter contracts of the new
// distributions.
func TestDistribution2Errors(t *testing.T) {
if _, err := WeibullDensity(1, 0, 1); err == nil || !strings.Contains(err.Error(), "shape k") {
t.Fatalf("k = 0: got %v, want the shape refusal", err)
}
if _, err := WeibullDensity(1, 2, 0); err == nil || !strings.Contains(err.Error(), "scale λ") {
t.Fatalf("λ = 0: got %v, want the scale refusal", err)
}
if _, err := WeibullDensity(math.NaN(), 2, 1); err == nil || !strings.Contains(err.Error(), "must be a number") {
t.Fatalf("NaN x: got %v, want the x refusal", err)
}
if _, err := WeibullCDF(1, -1, 1); err == nil || !strings.Contains(err.Error(), "shape k") {
t.Fatalf("negative k: got %v, want the shape refusal", err)
}
if _, err := WeibullQuantile(1.5, 2, 1); err == nil || !strings.Contains(err.Error(), "q must lie") {
t.Fatalf("q outside [0, 1]: got %v, want the q refusal", err)
}
if _, err := WeibullQuantile(0, 2, 1); err == nil || !strings.Contains(err.Error(), "no finite quantile") {
t.Fatalf("q = 0: got %v, want the q refusal", err)
}
if _, err := WeibullCDF(1, 2, math.Inf(1)); err == nil || !strings.Contains(err.Error(), "scale λ") {
t.Fatalf("λ = +Inf: got %v, want the scale refusal", err)
}
if v, err := WeibullCDF(math.Inf(1), 2, 1); err != nil || v != 1 {
t.Fatalf("WeibullCDF(+Inf) = %v, %v, want 1", v, err)
}
if _, err := LognormalQuantile(0.5, math.NaN(), 1); err == nil || !strings.Contains(err.Error(), "location μ") {
t.Fatalf("μ = NaN in the quantile: got %v, want the location refusal", err)
}
if _, err := LognormalQuantile(0.5, 0, math.Inf(1)); err == nil || !strings.Contains(err.Error(), "log-scale σ") {
t.Fatalf("σ = +Inf in the quantile: got %v, want the log-scale refusal", err)
}
if _, err := LognormalCDF(math.NaN(), 0, 1); err == nil || !strings.Contains(err.Error(), "must be a number") {
t.Fatalf("NaN x: got %v, want the x refusal", err)
}
if _, err := ParetoQuantile(0.5, 0, 2); err == nil || !strings.Contains(err.Error(), "scale x_m") {
t.Fatalf("x_m = 0 in the quantile: got %v, want the scale refusal", err)
}
if _, err := ParetoQuantile(0.5, 1, math.NaN()); err == nil || !strings.Contains(err.Error(), "tail index α") {
t.Fatalf("NaN α: got %v, want the tail-index refusal", err)
}
if _, err := ParetoCDF(math.NaN(), 1, 2); err == nil || !strings.Contains(err.Error(), "must be a number") {
t.Fatalf("NaN x: got %v, want the x refusal", err)
}
if _, err := WeibullQuantile(0.5, math.Inf(1), 1); err == nil || !strings.Contains(err.Error(), "shape k") {
t.Fatalf("k = +Inf in the quantile: got %v, want the shape refusal", err)
}
if _, err := WeibullQuantile(0.5, 2, math.Inf(-1)); err == nil || !strings.Contains(err.Error(), "scale λ") {
t.Fatalf("λ = −Inf in the quantile: got %v, want the scale refusal", err)
}
if _, err := LognormalDensity(math.NaN(), 0, 1); err == nil || !strings.Contains(err.Error(), "must be a number") {
t.Fatalf("NaN x in the density: got %v, want the x refusal", err)
}
if _, err := ParetoDensity(math.NaN(), 1, 2); err == nil || !strings.Contains(err.Error(), "must be a number") {
t.Fatalf("NaN x in the Pareto density: got %v, want the x refusal", err)
}
if _, err := ParetoDensity(1, math.Inf(1), 2); err == nil || !strings.Contains(err.Error(), "scale x_m") {
t.Fatalf("x_m = +Inf: got %v, want the scale refusal", err)
}
if _, err := NegativeBinomialCDF(0, 0, 0.5); err == nil || !strings.Contains(err.Error(), "r must be") {
t.Fatalf("r = 0 in the CDF: got %v, want the r refusal", err)
}
if _, err := DirichletMode([]float64{2}); err == nil || !strings.Contains(err.Error(), "at least two components") {
t.Fatalf("one component in the mode: got %v, want the component floor refusal", err)
}
if _, err := DirichletMode([]float64{2, math.Inf(1)}); err == nil || !strings.Contains(err.Error(), "alpha[1]") {
t.Fatalf("α = +Inf in the mode: got %v, want the concentration refusal", err)
}
if _, err := DirichletDraws(core.NewGenerator(1), 1, []float64{2}); err == nil || !strings.Contains(err.Error(), "at least two components") {
t.Fatalf("one component in the draws: got %v, want the component floor refusal", err)
}
if _, err := DirichletMean([]float64{2, math.NaN()}); err == nil || !strings.Contains(err.Error(), "alpha[1]") {
t.Fatalf("NaN α in the mean: got %v, want the concentration refusal", err)
}
if _, err := LognormalDensity(1, math.Inf(1), 1); err == nil || !strings.Contains(err.Error(), "location μ") {
t.Fatalf("μ = +Inf: got %v, want the location refusal", err)
}
if _, err := LognormalCDF(1, 0, -1); err == nil || !strings.Contains(err.Error(), "log-scale σ") {
t.Fatalf("σ < 0: got %v, want the log-scale refusal", err)
}
if _, err := LognormalQuantile(0.5, 0, 0); err == nil || !strings.Contains(err.Error(), "log-scale σ") {
t.Fatalf("σ = 0: got %v, want the log-scale refusal", err)
}
if _, err := ParetoDensity(1, 0, 2); err == nil || !strings.Contains(err.Error(), "scale x_m") {
t.Fatalf("x_m = 0: got %v, want the scale refusal", err)
}
if _, err := ParetoCDF(1, 1, math.Inf(1)); err == nil || !strings.Contains(err.Error(), "tail index α") {
t.Fatalf("α = +Inf: got %v, want the tail-index refusal", err)
}
if _, err := ParetoQuantile(1, 1, 2); err == nil || !strings.Contains(err.Error(), "no finite quantile") {
t.Fatalf("q = 1: got %v, want the q refusal", err)
}
if _, err := NegativeBinomialPMF(0, 0, 0.5); err == nil || !strings.Contains(err.Error(), "r must be") {
t.Fatalf("r = 0: got %v, want the r refusal", err)
}
if _, err := NegativeBinomialPMF(0, 3, 1.5); err == nil || !strings.Contains(err.Error(), "p must lie") {
t.Fatalf("p above 1 in the PMF: got %v, want the p refusal", err)
}
if _, err := NegativeBinomialCDF(0, 3, 0); err == nil || !strings.Contains(err.Error(), "p must lie") {
t.Fatalf("p = 0: got %v, want the p refusal", err)
}
if _, err := NegativeBinomialCDF(0, 3, 1); err == nil || !strings.Contains(err.Error(), "p must lie") {
t.Fatalf("p = 1: got %v, want the p refusal", err)
}
if _, err := NegativeBinomialQuantile(0.5, 0.5, 0); err == nil || !strings.Contains(err.Error(), "r must be") {
t.Fatalf("r = 0 in the quantile: got %v, want the r refusal", err)
}
if _, err := NegativeBinomialQuantile(0.5, 1.5, 3); err == nil || !strings.Contains(err.Error(), "p must lie") {
t.Fatalf("p above 1 in the quantile: got %v, want the p refusal", err)
}
if _, err := NegativeBinomialQuantile(-0.1, 0.5, 3); err == nil || !strings.Contains(err.Error(), "q must lie") {
t.Fatalf("q < 0: got %v, want the q refusal", err)
}
if _, err := DirichletDensity([]float64{2}, []float64{1}); err == nil || !strings.Contains(err.Error(), "at least two components") {
t.Fatalf("one component: got %v, want the component floor refusal", err)
}
if _, err := DirichletDensity([]float64{2, 3}, []float64{0.5}); err == nil || !strings.Contains(err.Error(), "components, the point") {
t.Fatalf("length mismatch: got %v, want the length refusal", err)
}
if _, err := DirichletDensity([]float64{2, 3}, []float64{0.5, 0.2}); err == nil || !strings.Contains(err.Error(), "must sum to 1") {
t.Fatalf("off-simplex point: got %v, want the simplex refusal", err)
}
if _, err := DirichletDensity([]float64{-1, 3}, []float64{0.5, 0.5}); err == nil || !strings.Contains(err.Error(), "alpha[0]") {
t.Fatalf("negative α: got %v, want the concentration refusal", err)
}
if _, err := DirichletMean([]float64{0, 3}); err == nil || !strings.Contains(err.Error(), "alpha[0]") {
t.Fatalf("α = 0 in the mean: got %v, want the concentration refusal", err)
}
if _, err := DirichletDraws(core.NewGenerator(1), 0, []float64{2, 3}); err == nil || !strings.Contains(err.Error(), "n must be") {
t.Fatalf("n = 0: got %v, want the n refusal", err)
}
}