diff --git a/cmd/internal-tools/analyze/censoring_test.go b/cmd/internal-tools/analyze/censoring_test.go index 5d7b49b..6a9f2ec 100644 --- a/cmd/internal-tools/analyze/censoring_test.go +++ b/cmd/internal-tools/analyze/censoring_test.go @@ -96,3 +96,23 @@ func TestRun_EffectSizeFollowsWhatCensoringDetermines(t *testing.T) { "and none of the early arm's twenty had violated by step 12", pair.A12) } } + +// The seed-matched contrast reads the same censored runs and reaches the same +// conclusion or it is not measuring the same thing. +func TestRun_PairedContrastFollowsWhatCensoringDetermines(t *testing.T) { + late := append(violatedAt(6, 5), violatedAt(14, 100)...) + earlyDirectory, lateDirectory := wallClockArms(t, stoppedShort(20, 12), late) + result := analyseCampaigns(t, "--paired", earlyDirectory, lateDirectory) + + paired := *result.Paired + if paired.First != "early" || paired.Second != "late" { + t.Fatalf("paired %s minus %s, want early minus late", paired.First, paired.Second) + } + if paired.Sign != 1 { + t.Errorf("sign %+d, want +1: the late arm is the one seen to violate first, in the six pairs "+ + "where the order is determined at all", paired.Sign) + } + if paired.A12 <= 0.5 { + t.Errorf("a12 within pairs %.4f, want above 0.5", paired.A12) + } +} diff --git a/cmd/internal-tools/analyze/paired.go b/cmd/internal-tools/analyze/paired.go index 16fbd57..f45eb15 100644 --- a/cmd/internal-tools/analyze/paired.go +++ b/cmd/internal-tools/analyze/paired.go @@ -6,178 +6,73 @@ import ( "slices" ) -type signedRankResult struct { - Pairs int `json:"pairs"` - NonZero int `json:"non_zero_pairs"` - Statistic float64 `json:"signed_rank_v"` - PValue float64 `json:"p_value"` - Exact bool `json:"exact"` -} - -// exactSignedRankLimit matches R's wilcox.test: the exact null distribution is -// used only below this many non-zero differences, and only when nothing is tied. -const exactSignedRankLimit = 50 - -// signedRank is the two-sided Wilcoxon signed-rank test over paired -// differences. Zero differences are dropped before ranking and the statistic is -// the sum of the ranks carried by the positive differences, which is the -// quantity R's wilcox.test calls V. Wilcoxon (1945), "Individual Comparisons by -// Ranking Methods", Biometrics Bulletin 1(6), 80-83. -func signedRank(differences []float64) signedRankResult { - result := signedRankResult{ - Pairs: len(differences), - Statistic: math.NaN(), - PValue: math.NaN(), - } - var magnitudes []float64 - var positive []bool - for _, difference := range differences { - if difference == 0 { - continue - } - magnitudes = append(magnitudes, math.Abs(difference)) - positive = append(positive, difference > 0) - } - result.NonZero = len(magnitudes) - if result.NonZero == 0 { - return result - } - ranks, tieGroups := midRanks(magnitudes) - statistic := 0.0 - for index, rank := range ranks { - if positive[index] { - statistic += rank - } - } - result.Statistic = statistic - - droppedZeros := len(differences) != result.NonZero - if len(tieGroups) == 0 && !droppedZeros && result.NonZero < exactSignedRankLimit { - result.Exact = true - result.PValue = exactSignedRankTwoSided(statistic, result.NonZero) - return result - } - result.PValue = normalSignedRankTwoSided(statistic, result.NonZero, tieGroups) - return result -} - -// normalSignedRankTwoSided follows the large-sample branch of R's wilcox.test: +// signTest is the exact two-sided sign test over matched pairs. Under the null +// that neither arm reaches its first violation sooner, a pair whose order the +// censoring determines falls either way with probability one half, so the count +// is binomial and the two-sided p-value doubles the smaller tail. Pairs left +// with no order carry no information and are not trials. // -// mean = n(n+1)/4 -// variance = n(n+1)(2n+1)/24 - sum(t^3 - t)/48 -// -// where t runs over the sizes of the groups tied on the absolute difference. -// The 0.5 shift toward the null mean is the continuity correction. -func normalSignedRankTwoSided(statistic float64, count int, tieGroups []int) float64 { - variance := signedRankVariance(count, tieGroups) - if variance <= 0 { - return 1 +// It is the seed-matched form of the comparison the unpaired test makes, and +// it is what the log-rank stratified by seed reduces to with one run per arm in +// each stratum. The magnitude-based alternatives are not available: a +// difference in steps needs both runs to have violated, and a paired test built +// on scores of censored times, the paired Prentice-Wilcoxon among them, is +// centred at zero under the null only when the two arms censor alike, which is +// exactly what the wall clock stops them from doing. +func signTest(favouringFirst, favouringSecond int) float64 { + trials := favouringFirst + favouringSecond + if trials == 0 { + return math.NaN() } - size := float64(count) - centered := statistic - size*(size+1)/4 - correction := 0.0 - switch { - case centered > 0: - correction = 0.5 - case centered < 0: - correction = -0.5 + tail := 0.0 + for count := 0; count <= min(favouringFirst, favouringSecond); count++ { + tail += math.Exp(logBinomialCoefficient(trials, count) - float64(trials)*math.Ln2) } - z := (centered - correction) / math.Sqrt(variance) - tail := math.Min(standardNormalUpperTail(z), standardNormalUpperTail(-z)) return math.Min(2*tail, 1) } -func signedRankVariance(count int, tieGroups []int) float64 { - size := float64(count) - tieAdjustment := 0.0 - for _, group := range tieGroups { - tied := float64(group) - tieAdjustment += tied*tied*tied - tied - } - return size*(size+1)*(2*size+1)/24 - tieAdjustment/48 -} - -// exactSignedRankTwoSided doubles the smaller exact tail, as R's wilcox.test -// does. -func exactSignedRankTwoSided(statistic float64, count int) float64 { - if statistic > float64(count)*float64(count+1)/4 { - return math.Min(2*exactSignedRankUpperTail(statistic, count), 1) - } - return math.Min(2*exactSignedRankLowerTail(statistic, count), 1) -} - -// exactSignedRankUpperTail is P(V >= statistic) under the null with no ties. -func exactSignedRankUpperTail(statistic float64, count int) float64 { - counts := exactSignedRankCounts(count) - total, tail := 0.0, 0.0 - for value, weight := range counts { - total += weight - if float64(value) >= statistic { - tail += weight - } - } - return tail / total -} - -func exactSignedRankLowerTail(statistic float64, count int) float64 { - counts := exactSignedRankCounts(count) - total, tail := 0.0, 0.0 - for value, weight := range counts { - total += weight - if float64(value) <= statistic { - tail += weight - } - } - return tail / total -} - -// exactSignedRankCounts returns the number of sign assignments producing each -// value of V from 0 to n(n+1)/2. V is the sum of the ranks held by the positive -// differences, so the count is a subset-sum tally over the ranks 1 to n. -func exactSignedRankCounts(count int) []float64 { - high := count * (count + 1) / 2 - table := make([]float64, high+1) - table[0] = 1 - for rank := 1; rank <= count; rank++ { - for sum := high; sum >= rank; sum-- { - if table[sum-rank] != 0 { - table[sum] += table[sum-rank] - } - } - } - return table +func logBinomialCoefficient(trials, chosen int) float64 { + all, _ := math.Lgamma(float64(trials + 1)) + picked, _ := math.Lgamma(float64(chosen + 1)) + rest, _ := math.Lgamma(float64(trials-chosen) + 1) + return all - picked - rest } // pairedComparison is the seed-matched contrast the actuation ablation reports. -// The difference is the first arm's steps to first violation less the second's, -// so a positive median means the second arm reached its first violation sooner, -// and Sign carries that direction as a number the decision rule can read. +// A pair is scored the way the unpaired comparison scores one, by which run +// outlived the other, so Sign is +1 when the second arm is the one seen to +// violate sooner across the pairs whose order censoring determines. type pairedComparison struct { - First string `json:"first"` - Second string `json:"second"` - Pairs int `json:"pairs"` - UnpairedSeeds []int64 `json:"unpaired_seeds,omitempty"` - MedianDifference float64 `json:"median_step_difference"` - Sign int `json:"sign"` - FirstSooner int `json:"first_sooner"` - SecondSooner int `json:"second_sooner"` - Tied int `json:"tied"` + First string `json:"first"` + Second string `json:"second"` + Pairs int `json:"pairs"` + UnpairedSeeds []int64 `json:"unpaired_seeds,omitempty"` + Sign int `json:"sign"` + FirstSooner int `json:"first_sooner"` + SecondSooner int `json:"second_sooner"` + // Unordered is the pairs the censoring leaves in no order, either because + // both runs ended clean or because the run that stopped first stopped before + // the other violated. They are not evidence either way and are not trials. + Unordered int `json:"unordered_pairs"` + // MedianDifference is in steps and is undefined unless some pair has both + // runs violating, which is the only shape a difference in steps can be read + // off. BothViolated says how many pairs it summarizes, because it describes + // those pairs and not the sample. + MedianDifference *float64 `json:"median_step_difference"` + BothViolated int `json:"both_violated_pairs"` // A12 is the within-pair form of the Vargha-Delaney effect size, the share - // of matched seeds on which the first arm took more steps, counting a tie as - // half. A matched design has no reason to compare the two arms as pooled - // bags of runs when each seed has a partner. + // of matched seeds on which the first arm took more steps, an unordered pair + // counting as half. A matched design has no reason to compare the two arms + // as pooled bags of runs when each seed has a partner. A12 float64 `json:"a12_within_pairs"` - Statistic float64 `json:"signed_rank_v"` PValue float64 `json:"p_value"` HolmPValue float64 `json:"holm_p_value"` - Exact bool `json:"exact"` } -// pairArms matches the two arms by seed and contrasts them pair by pair. -// Censored runs enter at the steps they ran, the same convention the unpaired -// comparison uses. A seed usable in one arm and not the other is named rather -// than dropped silently, because that is a host that lost a run and it is what -// the campaign manifest exists to make visible. +// pairArms matches the two arms by seed and contrasts them pair by pair. A seed +// usable in one arm and not the other is named rather than dropped silently, +// because that is a host that lost a run and it is what the campaign manifest +// exists to make visible. func pairArms(first, second arm) (pairedComparison, error) { firstBySeed, err := usableBySeed(first) if err != nil { @@ -192,7 +87,6 @@ func pairArms(first, second arm) (pairedComparison, error) { First: first.Name, Second: second.Name, A12: math.NaN(), - Statistic: math.NaN(), PValue: math.NaN(), HolmPValue: math.NaN(), } @@ -204,34 +98,38 @@ func pairArms(first, second arm) (pairedComparison, error) { comparison.UnpairedSeeds = append(comparison.UnpairedSeeds, seed) continue } - difference := observationOf(left, first.Budget).Steps - observationOf(right, second.Budget).Steps - differences = append(differences, difference) - switch { - case difference < 0: - comparison.FirstSooner++ - case difference > 0: + leftRun := observationOf(left, first.Budget) + rightRun := observationOf(right, second.Budget) + comparison.Pairs++ + switch outlives(leftRun, rightRun) { + case 1: comparison.SecondSooner++ + case -1: + comparison.FirstSooner++ default: - comparison.Tied++ + comparison.Unordered++ + } + if leftRun.Event && rightRun.Event { + comparison.BothViolated++ + differences = append(differences, leftRun.Steps-rightRun.Steps) } } - comparison.Pairs = len(differences) if comparison.Pairs == 0 { return comparison, nil } - comparison.MedianDifference = medianOf(differences) + if len(differences) > 0 { + median := medianOf(differences) + comparison.MedianDifference = &median + } switch { - case comparison.MedianDifference > 0: + case comparison.SecondSooner > comparison.FirstSooner: comparison.Sign = 1 - case comparison.MedianDifference < 0: + case comparison.FirstSooner > comparison.SecondSooner: comparison.Sign = -1 } - comparison.A12 = (float64(comparison.SecondSooner) + 0.5*float64(comparison.Tied)) / float64(comparison.Pairs) - test := signedRank(differences) - comparison.Statistic = test.Statistic - comparison.PValue = test.PValue - comparison.Exact = test.Exact + comparison.A12 = (float64(comparison.SecondSooner) + 0.5*float64(comparison.Unordered)) / float64(comparison.Pairs) + comparison.PValue = signTest(comparison.FirstSooner, comparison.SecondSooner) return comparison, nil } diff --git a/cmd/internal-tools/analyze/paired_test.go b/cmd/internal-tools/analyze/paired_test.go index 56d8ac6..b6950df 100644 --- a/cmd/internal-tools/analyze/paired_test.go +++ b/cmd/internal-tools/analyze/paired_test.go @@ -5,249 +5,52 @@ import ( "testing" ) -// Hollander and Wolfe (1973), 29f: Hamilton depression scale factor -// measurements on nine patients, first at admission and again after tranquilizer -// treatment. R's wilcox.test help page uses exactly these vectors as its paired -// example and reports -// -// wilcox.test(x, y, paired = TRUE, alternative = "greater") -// ## V = 40, p-value = 0.01953 -var ( - depressionAtAdmission = []float64{1.83, 0.50, 1.62, 2.48, 1.68, 1.88, 1.55, 3.06, 1.30} - depressionAfterOneWeek = []float64{0.878, 0.647, 0.598, 2.05, 1.06, 1.29, 1.06, 3.14, 1.29} -) - -func differencesOf(first, second []float64) []float64 { - differences := make([]float64, len(first)) - for index := range first { - differences[index] = first[index] - second[index] +// The two-sided sign test is R's binom.test(k, n) at p = 0.5, which is the +// doubled tail of a symmetric binomial and can be worked out by hand from the +// coefficients: 2 * sum(C(n, i), i <= min(k, n-k)) / 2^n. +func TestSignTest_MatchesTheBinomialTail(t *testing.T) { + cases := []struct { + first, second int + want float64 + }{ + {0, 10, 2.0 / 1024}, + {1, 9, 2 * 11.0 / 1024}, + {3, 7, 2 * 176.0 / 1024}, + {5, 5, 1}, + {0, 1, 1}, + {2, 0, 0.5}, } - return differences -} - -func TestSignedRank_MatchesPublishedDepressionResult(t *testing.T) { - result := signedRank(differencesOf(depressionAtAdmission, depressionAfterOneWeek)) - - if result.Statistic != 40 { - t.Errorf("statistic %v, want 40", result.Statistic) - } - if !result.Exact { - t.Error("expected the exact null distribution for nine untied differences") - } - upper := exactSignedRankUpperTail(40, 9) - if math.Abs(upper-0.01953) > 5e-6 { - t.Errorf("one-sided p-value %.6f, want 0.01953", upper) - } - if math.Abs(result.PValue-2*0.01953125) > 1e-9 { - t.Errorf("two-sided p-value %.6f, want %.6f", result.PValue, 2*0.01953125) - } -} - -// Reversing the pairs mirrors the statistic about n(n+1)/2 and leaves the -// two-sided p-value alone, which R reports as V = 5 on the same data. -func TestSignedRank_ReversedPairsMirrorTheStatistic(t *testing.T) { - forward := signedRank(differencesOf(depressionAtAdmission, depressionAfterOneWeek)) - reversed := signedRank(differencesOf(depressionAfterOneWeek, depressionAtAdmission)) - if reversed.Statistic != 5 { - t.Errorf("reversed statistic %v, want 5", reversed.Statistic) - } - if math.Abs(reversed.PValue-forward.PValue) > 1e-12 { - t.Errorf("reversed p-value %v, want %v", reversed.PValue, forward.PValue) - } -} - -// The exact null distribution must be a proper distribution: 2^n sign -// assignments in total, symmetric about n(n+1)/4. -func TestExactSignedRankCounts_FormAProperSymmetricDistribution(t *testing.T) { - counts := exactSignedRankCounts(8) - total := 0.0 - for _, count := range counts { - total += count - } - if total != 256 { - t.Errorf("counts sum to %v, want 2^8 = 256", total) - } - for index := range counts { - if counts[index] != counts[len(counts)-1-index] { - t.Errorf("count at %d is %v but %v at the mirrored point", index, counts[index], counts[len(counts)-1-index]) + for _, test := range cases { + got := signTest(test.first, test.second) + if math.Abs(got-test.want) > 1e-12 { + t.Errorf("sign test on %d against %d gives %v, want %v", test.first, test.second, got, test.want) + } + if reversed := signTest(test.second, test.first); math.Abs(reversed-got) > 1e-12 { + t.Errorf("sign test on %d against %d gives %v reversed and %v forward", + test.first, test.second, reversed, got) } } } -// The tie-corrected variance is checked against the exact permutation variance -// of the statistic, computed here by enumerating every sign assignment over the -// observed midranks. That is an independent calculation rather than a second -// call into the implementation under test. -func TestSignedRankVariance_MatchesExactPermutationVariance(t *testing.T) { - cases := [][]float64{ - {1, -2, 3, -4, 5, 6, -7, 8}, - {12, -12, 12, 12, -5, 5, 30, -30}, - {-400, 400, 400, -400, 400, 400, 400, 400}, - {3, 3, 3, -3, -3, 7, 7, 9, 9}, +// A campaign runs tens of seeds, not tens of thousands, but the tail is summed +// through log-gamma rather than through factorials so that a lopsided family +// stays a number rather than becoming an overflow. +func TestSignTest_LargeCountsStayFinite(t *testing.T) { + if got := signTest(0, 200); got <= 0 || got > 1e-59 { + t.Errorf("sign test on 0 against 200 gives %v, want a positive value around 2^-199", got) } - for _, differences := range cases { - magnitudes := make([]float64, len(differences)) - for index, difference := range differences { - magnitudes[index] = math.Abs(difference) - } - ranks, tieGroups := midRanks(magnitudes) - mean, variance := permutationMomentsOfSignedRank(ranks) - size := float64(len(ranks)) - if expected := size * (size + 1) / 4; math.Abs(mean-expected) > 1e-9 { - t.Errorf("permutation mean %v for %v, want %v", mean, differences, expected) - } - if got := signedRankVariance(len(ranks), tieGroups); math.Abs(got-variance) > 1e-9 { - t.Errorf("variance %v for %v, want the permutation variance %v", got, differences, variance) - } + if got := signTest(100, 100); math.Abs(got-1) > 1e-12 { + t.Errorf("sign test on an even split gives %v, want 1", got) } } -// With ties present the normal approximation is the only branch available, so -// it is checked against the exact permutation p-value of the same statistic on -// the same data. -func TestSignedRank_TiedDifferencesTrackTheExactPermutationPValue(t *testing.T) { - cases := [][]float64{ - {40, 40, 40, -12, 33, 40, 40, -3, 40, 21, 40, 40}, - {-5, -5, -5, -5, 9, 9, 2, 2, -1, -1, 40, 40}, - {40, 40, -12, 33, 40, -3, 21, 40, 15, -9}, - } - for _, differences := range cases { - result := signedRank(differences) - if result.Exact { - t.Errorf("%v used the exact null distribution despite ties", differences) - } - exact := permutationSignedRankTwoSided(differences) - if math.Abs(result.PValue-exact) > 0.03 { - t.Errorf("normal approximation p %.4f for %v, want near the permutation p %.4f", - result.PValue, differences, exact) - } +func TestSignTest_NoOrderedPairHasNoTest(t *testing.T) { + if got := signTest(0, 0); !math.IsNaN(got) { + t.Errorf("sign test with nothing to test gives %v, want undefined", got) } } -// Every difference the same size is the degenerate end of the tie correction, -// and the normal approximation is genuinely poor there: the exact randomization -// p-value on this sample is 0.3438 against the approximation's 0.2273. The tool -// keeps R's formula rather than the randomization p-value so that a reviewer -// running wilcox.test on the same differences reads the same number, and the -// expected value here is that published formula worked through by hand: -// -// n = 10, one tied group of 10, V = 7 * 5.5 = 38.5 -// mean = 10 * 11 / 4 = 27.5 -// variance = 10 * 11 * 21 / 24 - (10^3 - 10) / 48 = 96.25 - 20.625 = 75.625 -// z = (38.5 - 27.5 - 0.5) / sqrt(75.625) -func TestSignedRank_EveryDifferenceTheSameSizeFollowsTheDocumentedFormula(t *testing.T) { - differences := []float64{7, 7, 7, 7, 7, 7, 7, -7, -7, -7} - result := signedRank(differences) - - if result.Statistic != 38.5 { - t.Errorf("statistic %v, want 38.5", result.Statistic) - } - if got := signedRankVariance(10, []int{10}); math.Abs(got-75.625) > 1e-12 { - t.Errorf("variance %v, want 75.625", got) - } - expected := 2 * standardNormalUpperTail(10.5/math.Sqrt(75.625)) - if math.Abs(result.PValue-expected) > 1e-12 { - t.Errorf("p-value %v, want %v", result.PValue, expected) - } - if randomization := permutationSignedRankTwoSided(differences); math.Abs(randomization-0.3438) > 5e-4 { - t.Errorf("randomization p-value %.4f, want 0.3438", randomization) - } -} - -// R drops zero differences before ranking and tests what is left, so a pair -// where both arms took the same number of steps carries no direction and must -// not be ranked as though it did. -func TestSignedRank_ZeroDifferencesAreDropped(t *testing.T) { - result := signedRank([]float64{0, 0, 3, -1, 2}) - if result.Pairs != 5 || result.NonZero != 3 { - t.Errorf("pairs %d non-zero %d, want 5 and 3", result.Pairs, result.NonZero) - } - // Ranks over |{3, 1, 2}| are 3, 1, 2, and the positive differences hold 3 - // and 2. - if result.Statistic != 5 { - t.Errorf("statistic %v, want 5", result.Statistic) - } - if result.Exact { - t.Error("used the exact null distribution despite dropped zeros") - } -} - -func TestSignedRank_EveryDifferenceZero(t *testing.T) { - result := signedRank([]float64{0, 0, 0}) - if result.NonZero != 0 { - t.Errorf("non-zero pairs %d, want 0", result.NonZero) - } - if !math.IsNaN(result.PValue) || !math.IsNaN(result.Statistic) { - t.Errorf("result %+v, want everything undefined", result) - } -} - -// permutationMomentsOfSignedRank enumerates every sign assignment and returns -// the mean and variance of the statistic over them. -func permutationMomentsOfSignedRank(ranks []float64) (float64, float64) { - values := signedRankPermutationValues(ranks) - mean := 0.0 - for _, value := range values { - mean += value - } - mean /= float64(len(values)) - variance := 0.0 - for _, value := range values { - variance += (value - mean) * (value - mean) - } - return mean, variance / float64(len(values)) -} - -// permutationSignedRankTwoSided is the exact randomization p-value: the -// proportion of sign assignments whose statistic is at least as far from the -// null mean as the observed one. -func permutationSignedRankTwoSided(differences []float64) float64 { - magnitudes := make([]float64, 0, len(differences)) - observed := 0.0 - for _, difference := range differences { - if difference == 0 { - continue - } - magnitudes = append(magnitudes, math.Abs(difference)) - } - ranks, _ := midRanks(magnitudes) - position := 0 - for _, difference := range differences { - if difference == 0 { - continue - } - if difference > 0 { - observed += ranks[position] - } - position++ - } - size := float64(len(ranks)) - mean := size * (size + 1) / 4 - values := signedRankPermutationValues(ranks) - extreme := 0 - for _, value := range values { - if math.Abs(value-mean) >= math.Abs(observed-mean)-1e-9 { - extreme++ - } - } - return float64(extreme) / float64(len(values)) -} - -func signedRankPermutationValues(ranks []float64) []float64 { - values := make([]float64, 0, 1< 1e-12 { t.Errorf("a12 within pairs %v, want %v", comparison.A12, 2.5/3) } + // Only seeds 1 and 3 have a difference in steps to take a median of, 20 and + // 0: the pair holding a clean run has no difference either arm supports. + if comparison.BothViolated != 2 || comparison.MedianDifference == nil || *comparison.MedianDifference != 10 { + t.Errorf("median difference %v over %d pair(s), want 10 over 2", + comparison.MedianDifference, comparison.BothViolated) + } + if want := signTest(0, 2); comparison.PValue != want { + t.Errorf("p %v, want the sign test's %v over the two ordered pairs", comparison.PValue, want) + } +} + +// Two clean runs are two runs that were still going when they stopped, whatever +// step each stopped on, so the pair says nothing and is not a trial. +func TestPairArms_PairsOfCleanRunsAreNotEvidence(t *testing.T) { + early := arm{Name: "early", Budget: 400, Runs: []classifiedRun{cleanRun(1, 12), cleanRun(2, 14)}} + late := arm{Name: "late", Budget: 400, Runs: []classifiedRun{cleanRun(1, 400), cleanRun(2, 380)}} + + comparison, err := pairArms(early, late) + if err != nil { + t.Fatal(err) + } + if comparison.Unordered != 2 || comparison.Sign != 0 { + t.Errorf("comparison %+v, want both pairs unordered and no direction", comparison) + } + if !math.IsNaN(comparison.PValue) { + t.Errorf("p %v, want undefined with no ordered pair", comparison.PValue) + } + if comparison.MedianDifference != nil { + t.Errorf("median difference %v, want undefined where no pair has two violations", + *comparison.MedianDifference) + } + if comparison.A12 != 0.5 { + t.Errorf("a12 within pairs %v, want 0.5", comparison.A12) + } } // A run excluded as missing data cannot be paired against anything, and the @@ -344,8 +177,9 @@ func TestPairArms_DirectionReversesWithTheArms(t *testing.T) { if forward.Sign != 1 || reversed.Sign != -1 { t.Errorf("signs %+d and %+d, want +1 then -1", forward.Sign, reversed.Sign) } - if forward.MedianDifference != -reversed.MedianDifference { - t.Errorf("median differences %v and %v, want opposites", forward.MedianDifference, reversed.MedianDifference) + if *forward.MedianDifference != -*reversed.MedianDifference { + t.Errorf("median differences %v and %v, want opposites", + *forward.MedianDifference, *reversed.MedianDifference) } if math.Abs(forward.PValue-reversed.PValue) > 1e-12 { t.Errorf("p-values %v and %v, want the same two-sided value", forward.PValue, reversed.PValue) diff --git a/cmd/internal-tools/analyze/planted_test.go b/cmd/internal-tools/analyze/planted_test.go index 2dde2fd..58c8484 100644 --- a/cmd/internal-tools/analyze/planted_test.go +++ b/cmd/internal-tools/analyze/planted_test.go @@ -581,7 +581,12 @@ func TestPlanted_PairedComparisonRecoversTheShiftAndItsSign(t *testing.T) { source := rand.New(rand.NewSource(5150)) var post, pre []plantedRun - var differences []float64 + // Arms are contrasted in the order their labels sort, so the plant is stated + // the same way: post-repair against pre-repair. A seed where the unshifted + // run violated is a pair the shifted run is known to have outlived, whether + // it violated later or ran on clean; a seed where neither violated is a pair + // with no order. + firstSooner, bothViolated, unordered := 0, 0, 0 for seed := int64(1); seed <= 30; seed++ { fast := base.draw(seed, source) slow := plantedRun{seed: seed, steps: fast.steps + shift, violated: fast.violated} @@ -590,13 +595,18 @@ func TestPlanted_PairedComparisonRecoversTheShiftAndItsSign(t *testing.T) { } post = append(post, fast) pre = append(pre, slow) - // Arms are contrasted in the order their labels sort, so the planted - // difference is stated the same way: post-repair less pre-repair. - differences = append(differences, float64(recordedStep(fast, budget)-recordedStep(slow, budget))) + switch { + case !fast.violated: + unordered++ + case slow.violated: + bothViolated++ + firstSooner++ + default: + firstSooner++ + } } - expected := medianOf(differences) - if expected >= 0 { - t.Fatalf("planted median difference %v, want the shifted arm to violate later", expected) + if bothViolated == 0 { + t.Fatal("no pair has two violations, so the planted shift is nowhere the analysis can read it") } root := t.TempDir() @@ -618,17 +628,24 @@ func TestPlanted_PairedComparisonRecoversTheShiftAndItsSign(t *testing.T) { if paired.Pairs != 30 { t.Errorf("%d pairs, want 30", paired.Pairs) } - if paired.MedianDifference != expected { - t.Errorf("median difference %v, want the planted %v", paired.MedianDifference, expected) + // Every pair where both runs violated was shifted by exactly the plant, so + // the median over them is the plant itself rather than a mixture of it with + // the step counts censored runs never reached. + if paired.BothViolated != bothViolated { + t.Errorf("%d pair(s) with two violations, want %d", paired.BothViolated, bothViolated) + } + if paired.MedianDifference == nil || *paired.MedianDifference != -shift { + t.Errorf("median difference %v, want the planted %d", paired.MedianDifference, -shift) } if paired.Sign != -1 { t.Errorf("sign %+d, want -1 for the arm that violated sooner", paired.Sign) } - if paired.SecondSooner != 0 { - t.Errorf("%d pairs favour the shifted arm, want none: every pair was shifted the same way", paired.SecondSooner) + if paired.FirstSooner != firstSooner || paired.SecondSooner != 0 || paired.Unordered != unordered { + t.Errorf("counts %+v, want %d favouring the shifted arm, none the other way and %d unordered", + paired, firstSooner, unordered) } - if paired.A12 != 0 { - t.Errorf("a12 within pairs %v, want 0 where no pair favours the first arm", paired.A12) + if want := 0.5 * float64(unordered) / 30; paired.A12 != want { + t.Errorf("a12 within pairs %v, want %v where no pair favours the first arm", paired.A12, want) } if paired.PValue > 0.001 { t.Errorf("p-value %v for a shift planted in every pair", paired.PValue) @@ -638,13 +655,6 @@ func TestPlanted_PairedComparisonRecoversTheShiftAndItsSign(t *testing.T) { } } -func recordedStep(run plantedRun, budget int) int { - if run.violated { - return run.steps - } - return budget -} - // The paired test is the ablation's decision rule, so a null there has to stay a // null: two arms drawn from the same model, matched by seed, must not report a // difference more often than the level allows. diff --git a/cmd/internal-tools/analyze/report.go b/cmd/internal-tools/analyze/report.go index 3eecd95..d915791 100644 --- a/cmd/internal-tools/analyze/report.go +++ b/cmd/internal-tools/analyze/report.go @@ -140,15 +140,16 @@ func writeReport(result analysis, out io.Writer) { } func writePaired(out io.Writer, comparison pairedComparison) { - fmt.Fprintf(out, "\npaired per-seed difference, %s minus %s, censored runs held at the steps they ran\n", + fmt.Fprintf(out, "\npaired per-seed contrast, %s against %s, each pair scored by which run outlived the other\n", comparison.First, comparison.Second) - fmt.Fprintf(out, "%d seed pair(s): %s sooner in %d, %s sooner in %d, tied in %d\n", + fmt.Fprintf(out, "%d seed pair(s): %s sooner in %d, %s sooner in %d, left in no order by censoring in %d\n", comparison.Pairs, comparison.First, comparison.FirstSooner, - comparison.Second, comparison.SecondSooner, comparison.Tied) - fmt.Fprintf(out, "median difference %+.1f steps, sign %+d, a12 within pairs %.3f\n", - comparison.MedianDifference, comparison.Sign, comparison.A12) - fmt.Fprintf(out, "wilcoxon signed-rank v %.1f, p %s, holm p %s\n", - comparison.Statistic, formatPValue(comparison.PValue), formatPValue(comparison.HolmPValue)) + comparison.Second, comparison.SecondSooner, comparison.Unordered) + fmt.Fprintf(out, "median difference %s over the %d pair(s) where both runs violated, sign %+d, a12 within pairs %.3f\n", + formatStepDifference(comparison.MedianDifference), comparison.BothViolated, comparison.Sign, comparison.A12) + fmt.Fprintf(out, "sign test over the %d ordered pair(s), p %s, holm p %s\n", + comparison.FirstSooner+comparison.SecondSooner, + formatPValue(comparison.PValue), formatPValue(comparison.HolmPValue)) if len(comparison.UnpairedSeeds) > 0 { fmt.Fprintf(out, "%d seed(s) usable in one arm only and left out of the pairing: %v\n", len(comparison.UnpairedSeeds), comparison.UnpairedSeeds) @@ -177,6 +178,13 @@ func formatMedian(value *float64) string { return strconv.FormatFloat(*value, 'f', -1, 64) } +func formatStepDifference(value *float64) string { + if value == nil { + return "undefined" + } + return fmt.Sprintf("%+.1f steps", *value) +} + // formatActions marks a denominator with actions whose producer nothing names, // because the rate beside it then divides by a count that may include the // login the spec's setup drove.