feat(analyze): survival analysis over campaign directories

Steps to first violation with clean runs right-censored at the budget, since
per-run yield is a binary at 11 to 45 percent and separating two arms on it
would need roughly 80 runs per arm. Kaplan-Meier, log-rank, Wilcoxon rank-sum
with Vargha-Delaney A12, Holm within each family.

A hand-rolled log-rank that is subtly wrong is a silent-wrong-number generator
and would be believed, so every statistic is validated against a published
worked example with the source named in the test: R survdiff on aml, Freireich
6-MP, Hollander and Wolfe 1973 for the rank sum, printed p.adjust output for
Holm. Two could not be: the k>2 log-rank, guarded by calibration instead, and
the tie-corrected variance, checked against an exact permutation variance.

Failed and timed-out runs are excluded as missing data and counted by reason,
never treated as censored observations, which would bias the result.

Claude-Session: https://claude.ai/code/session_01A5KmftdEJ49A9z5mF5ESrX

cmd/internal-tools/analyze/survival.go

@@ -0,0 +1,241 @@
+package main
+
+import (
+	"math"
+	"slices"
+)
+
+// observation is one run reduced to what the survival analysis needs: the step
+// count at which it left the risk set, and whether it left because a violation
+// was found (an event) or because it exhausted the step budget without one
+// (right-censored). A run that failed or timed out is neither and never
+// reaches this type.
+type observation struct {
+	Steps float64
+	Event bool
+}
+
+type survivalPoint struct {
+	Steps    float64 `json:"steps"`
+	AtRisk   int     `json:"at_risk"`
+	Events   int     `json:"events"`
+	Censored int     `json:"censored"`
+	Survival float64 `json:"survival"`
+}
+
+// kaplanMeier is the product-limit estimate of Kaplan and Meier (1958), one row
+// per distinct observed step count. Runs censored at a step count tied with an
+// event are counted in the risk set for that event, which is the standard
+// convention.
+func kaplanMeier(observations []observation) []survivalPoint {
+	if len(observations) == 0 {
+		return nil
+	}
+	remaining := len(observations)
+	survival := 1.0
+	var curve []survivalPoint
+	for _, steps := range distinctSteps(observations) {
+		events, censored := 0, 0
+		for _, item := range observations {
+			if item.Steps != steps {
+				continue
+			}
+			if item.Event {
+				events++
+			} else {
+				censored++
+			}
+		}
+		atRisk := remaining
+		if events > 0 {
+			survival *= 1 - float64(events)/float64(atRisk)
+		}
+		curve = append(curve, survivalPoint{
+			Steps:    steps,
+			AtRisk:   atRisk,
+			Events:   events,
+			Censored: censored,
+			Survival: survival,
+		})
+		remaining -= events + censored
+	}
+	return curve
+}
+
+// medianSurvival is the smallest step count at which the estimate falls to or
+// below one half. It is undefined whenever fewer than half the runs violate,
+// and the second return value says so: substituting a mean there would report a
+// number the data does not contain.
+func medianSurvival(curve []survivalPoint) (float64, bool) {
+	for _, point := range curve {
+		if point.Survival <= 0.5 {
+			return point.Steps, true
+		}
+	}
+	return 0, false
+}
+
+func distinctSteps(observations []observation) []float64 {
+	steps := make([]float64, 0, len(observations))
+	for _, item := range observations {
+		steps = append(steps, item.Steps)
+	}
+	slices.Sort(steps)
+	return slices.Compact(steps)
+}
+
+type logRankResult struct {
+	Groups           []string  `json:"groups"`
+	Sizes            []int     `json:"sizes"`
+	Observed         []float64 `json:"observed"`
+	Expected         []float64 `json:"expected"`
+	ChiSquare        float64   `json:"chi_square"`
+	DegreesOfFreedom int       `json:"degrees_of_freedom"`
+	PValue           float64   `json:"p_value"`
+}
+
+// logRank is the k-sample Mantel-Haenszel log-rank test. At every distinct
+// event time it contrasts observed with expected events under the null of equal
+// hazards, then combines the k-1 independent differences through their
+// covariance matrix: chi-square = U' V^-1 U on k-1 degrees of freedom.
+// Mantel (1966); Peto and Peto (1972); Klein and Moeschberger, Survival
+// Analysis, 2nd ed., section 7.3.
+func logRank(names []string, groups [][]observation) logRankResult {
+	var keptNames []string
+	var kept [][]observation
+	for index, group := range groups {
+		if len(group) == 0 {
+			continue
+		}
+		keptNames = append(keptNames, names[index])
+		kept = append(kept, group)
+	}
+	names, groups = keptNames, kept
+
+	count := len(groups)
+	result := logRankResult{
+		Groups:           names,
+		Sizes:            make([]int, count),
+		Observed:         make([]float64, count),
+		Expected:         make([]float64, count),
+		DegreesOfFreedom: count - 1,
+		PValue:           math.NaN(),
+	}
+	if count < 2 {
+		return result
+	}
+	var pooled []observation
+	for index, group := range groups {
+		result.Sizes[index] = len(group)
+		pooled = append(pooled, group...)
+	}
+
+	covariance := make([][]float64, count)
+	for index := range covariance {
+		covariance[index] = make([]float64, count)
+	}
+	atRisk := make([]float64, count)
+	deaths := make([]float64, count)
+
+	for _, steps := range distinctSteps(pooled) {
+		totalAtRisk, totalDeaths := 0.0, 0.0
+		for index, group := range groups {
+			atRisk[index], deaths[index] = 0, 0
+			for _, item := range group {
+				if item.Steps >= steps {
+					atRisk[index]++
+				}
+				if item.Steps == steps && item.Event {
+					deaths[index]++
+				}
+			}
+			totalAtRisk += atRisk[index]
+			totalDeaths += deaths[index]
+		}
+		if totalDeaths == 0 {
+			continue
+		}
+		for index := range groups {
+			result.Observed[index] += deaths[index]
+			result.Expected[index] += totalDeaths * atRisk[index] / totalAtRisk
+		}
+		if totalAtRisk <= 1 {
+			continue
+		}
+		scale := totalDeaths * (totalAtRisk - totalDeaths) / (totalAtRisk - 1)
+		for row := range groups {
+			share := atRisk[row] / totalAtRisk
+			covariance[row][row] += scale * share * (1 - share)
+			for column := range groups {
+				if column == row {
+					continue
+				}
+				covariance[row][column] -= scale * share * atRisk[column] / totalAtRisk
+			}
+		}
+	}
+
+	reduced := make([][]float64, count-1)
+	difference := make([]float64, count-1)
+	for row := 0; row < count-1; row++ {
+		reduced[row] = make([]float64, count-1)
+		copy(reduced[row], covariance[row][:count-1])
+		difference[row] = result.Observed[row] - result.Expected[row]
+	}
+	solution, ok := solveLinearSystem(reduced, difference)
+	if !ok {
+		result.ChiSquare = 0
+		result.PValue = 1
+		return result
+	}
+	statistic := 0.0
+	for index := range difference {
+		statistic += difference[index] * solution[index]
+	}
+	if statistic < 0 || math.IsNaN(statistic) {
+		statistic = 0
+	}
+	result.ChiSquare = statistic
+	result.PValue = chiSquareUpperTail(statistic, result.DegreesOfFreedom)
+	return result
+}
+
+// solveLinearSystem solves matrix*x = vector by Gaussian elimination with
+// partial pivoting, reporting failure rather than a value when the matrix is
+// singular, which is what a group with no events at all produces.
+func solveLinearSystem(matrix [][]float64, vector []float64) ([]float64, bool) {
+	size := len(vector)
+	work := make([][]float64, size)
+	for row := range work {
+		work[row] = make([]float64, size+1)
+		copy(work[row], matrix[row])
+		work[row][size] = vector[row]
+	}
+	for column := 0; column < size; column++ {
+		pivot := column
+		for row := column + 1; row < size; row++ {
+			if math.Abs(work[row][column]) > math.Abs(work[pivot][column]) {
+				pivot = row
+			}
+		}
+		if math.Abs(work[pivot][column]) < 1e-12 {
+			return nil, false
+		}
+		work[column], work[pivot] = work[pivot], work[column]
+		for row := column + 1; row < size; row++ {
+			factor := work[row][column] / work[column][column]
+			for next := column; next <= size; next++ {
+				work[row][next] -= factor * work[column][next]
+			}
+		}
+	}
+	solution := make([]float64, size)
+	for row := size - 1; row >= 0; row-- {
+		total := work[row][size]
+		for column := row + 1; column < size; column++ {
+			total -= work[row][column] * solution[column]
+		}
+		solution[row] = total / work[row][row]
+	}
+	return solution, true
+}