tensor/stats/multipletest_test.go

// Copyright (c) 2026 Petr Balvín <opensource@petrbalvin.org> (https://petrbalvin.org)
// SPDX-License-Identifier: MIT

package stats

import (
	"math"
	"testing"
)

// The worked vector for all three corrections: sorted it reads
// 0.005, 0.01, 0.03, 0.04 and every adjusted value below is computed
// by hand from that order.
var multipleTestP = []float64{0.01, 0.04, 0.03, 0.005}

// TestBonferroniAdjusted pins the hand values m·p: (0.04, 0.16, 0.12,
// 0.02), and the clamp at 1.
func TestBonferroniAdjusted(t *testing.T) {
	got, err := Bonferroni(multipleTestP)
	if err != nil {
		t.Fatalf("Bonferroni: %v", err)
	}
	want := []float64{0.04, 0.16, 0.12, 0.02}
	for i := range want {
		if math.Abs(got[i]-want[i]) > 1e-15 {
			t.Fatalf("Bonferroni[%d] = %.17g, want %.17g", i, got[i], want[i])
		}
	}
	clamped, err := Bonferroni([]float64{0.5, 0.6})
	if err != nil || clamped[0] != 1 || clamped[1] != 1 {
		t.Fatalf("clamped Bonferroni = %v, want [1, 1]", clamped)
	}
}

// TestHolmAdjusted pins the step-down walk. Sorted, the multipliers
// 4, 3, 2, 1 give 0.02, 0.03, 0.06, 0.06 after the running maximum,
// mapped back through the original order to (0.03, 0.06, 0.06, 0.02).
func TestHolmAdjusted(t *testing.T) {
	got, err := Holm(multipleTestP)
	if err != nil {
		t.Fatalf("Holm: %v", err)
	}
	want := []float64{0.03, 0.06, 0.06, 0.02}
	for i := range want {
		if math.Abs(got[i]-want[i]) > 1e-15 {
			t.Fatalf("Holm[%d] = %.17g, want %.17g", i, got[i], want[i])
		}
	}
}

// TestBenjaminiHochbergAdjusted pins the step-up walk. Sorted, from
// the top: 1·0.04 = 0.04, min(0.04, 4/3·0.03) = 0.04,
// min(0.04, 2·0.01) = 0.02, min(0.02, 4·0.005) = 0.02, so the original
// order carries (0.02, 0.04, 0.04, 0.02). Evenly spaced p-values all
// receive the largest raw p as their q-value, the identity 5/j·j/100
// = 0.05 on (0.01 .. 0.05).
func TestBenjaminiHochbergAdjusted(t *testing.T) {
	got, err := BenjaminiHochberg(multipleTestP)
	if err != nil {
		t.Fatalf("BenjaminiHochberg: %v", err)
	}
	want := []float64{0.02, 0.04, 0.04, 0.02}
	for i := range want {
		if math.Abs(got[i]-want[i]) > 1e-15 {
			t.Fatalf("BenjaminiHochberg[%d] = %.17g, want %.17g", i, got[i], want[i])
		}
	}
	uniform, err := BenjaminiHochberg([]float64{0.01, 0.02, 0.03, 0.04, 0.05})
	if err != nil {
		t.Fatalf("BenjaminiHochberg uniform: %v", err)
	}
	for i, v := range uniform {
		if math.Abs(v-0.05) > 1e-15 {
			t.Fatalf("uniform q[%d] = %.17g, want 0.05", i, v)
		}
	}
}

// TestBenjaminiHochbergLiterature pins the fourteen p-values tabulated
// in Benjamini and Hochberg (1995), Controlling the false discovery
// rate, the classic worked example, with the q-values derived by hand
// from the sorted vector: q_i = min over j ≥ i of 14·p_j/j. From the
// top the raw products fall 0.7590, 14·0.6528/13 = 0.7030,
// 14·0.5719/12 = 0.6672, 14·0.4262/11 = 0.5424, 14·0.3240/10 = 0.4536
// and then 0.0714 on rank 9, so the running minimum holds 0.4536
// across ranks 10 and 11, 0.6672 on rank 12 and 0.7030 on rank 13;
// below, 0.0602 on rank 8 and the running minimum carries 0.0596 from
// rank 7 back over rank 6 (whose raw product is 0.0648).
func TestBenjaminiHochbergLiterature(t *testing.T) {
	p := []float64{
		0.0001, 0.0004, 0.0019, 0.0095, 0.0201, 0.0278, 0.0298,
		0.0344, 0.0459, 0.3240, 0.4262, 0.5719, 0.6528, 0.7590,
	}
	got, err := BenjaminiHochberg(p)
	if err != nil {
		t.Fatalf("BenjaminiHochberg: %v", err)
	}
	want := []float64{
		0.0014, 0.0028, 14 * 0.0019 / 3, 0.03325, 0.05628, 0.0596, 0.0596,
		0.0602, 0.0714, 0.4536, 14 * 0.4262 / 11, 14 * 0.5719 / 12, 14 * 0.6528 / 13, 0.7590,
	}
	for i := range want {
		if math.Abs(got[i]-want[i]) > 1e-14 {
			t.Fatalf("q[%d] = %.17g, want %.17g", i, got[i], want[i])
		}
	}
}

// TestMultipleTestProperties holds the contract every correction
// states: adjusted values never fall below their raw p-values, never
// leave [0, 1], and preserve the order of the raw p-values.
func TestMultipleTestProperties(t *testing.T) {
	raw := []float64{0.5, 0.002, 0.2, 0.04, 0.009, 0.7, 0.0001, 0.13}
	corrections := []struct {
		name  string
		apply func([]float64) ([]float64, error)
	}{
		{"Bonferroni", Bonferroni},
		{"Holm", Holm},
		{"BenjaminiHochberg", BenjaminiHochberg},
	}
	for _, c := range corrections {
		got, err := c.apply(raw)
		if err != nil {
			t.Fatalf("%s: %v", c.name, err)
		}
		for i := range raw {
			if got[i] < raw[i] {
				t.Fatalf("%s[%d] = %v sits below the raw %v", c.name, i, got[i], raw[i])
			}
			if got[i] < 0 || got[i] > 1 {
				t.Fatalf("%s[%d] = %v leaves [0, 1]", c.name, i, got[i])
			}
			for j := range raw {
				if raw[i] <= raw[j] && got[i] > got[j] {
					t.Fatalf("%s inverts the order at (%d, %d): %v raw %v vs %v raw %v",
						c.name, i, j, got[i], raw[i], got[j], raw[j])
				}
			}
		}
	}
	// Zeros stay zeros under every correction.
	zeroed, err := Holm([]float64{0, 0.2})
	if err != nil || zeroed[0] != 0 {
		t.Fatalf("Holm on a zero p-value = %v, %v, want 0 untouched", zeroed, err)
	}
	// Equal p-values receive equal adjustments through the running
	// extreme.
	tied, _ := Holm([]float64{0.05, 0.05})
	if tied[0] != tied[1] {
		t.Fatalf("tied Holm adjustments differ: %v", tied)
	}
}

// TestMultipleTestErrors pins the input contract on all three
// corrections.
func TestMultipleTestErrors(t *testing.T) {
	for _, apply := range []func([]float64) ([]float64, error){Bonferroni, Holm, BenjaminiHochberg} {
		if _, err := apply(nil); err == nil {
			t.Fatal("an empty vector: want an error")
		}
		if _, err := apply([]float64{0.1, math.NaN()}); err == nil {
			t.Fatal("a NaN p-value: want an error")
		}
		if _, err := apply([]float64{0.1, math.Inf(1)}); err == nil {
			t.Fatal("an infinite p-value: want an error")
		}
		if _, err := apply([]float64{1.5}); err == nil {
			t.Fatal("p above 1: want an error")
		}
		if _, err := apply([]float64{-0.1}); err == nil {
			t.Fatal("a negative p-value: want an error")
		}
	}
}