mirror of
https://github.com/priyanshujain/sanderling.git
synced 2026-10-02 19:17:10 +00:00
pj
019d608f65
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
84 lines
2.1 KiB
Go
84 lines
2.1 KiB
Go
package main
|
|
|
|
import "math"
|
|
|
|
// standardNormalUpperTail is P(Z > z) for a standard normal Z.
|
|
func standardNormalUpperTail(z float64) float64 {
|
|
return 0.5 * math.Erfc(z/math.Sqrt2)
|
|
}
|
|
|
|
// chiSquareUpperTail is P(X > x) for a chi-square variate with the given
|
|
// degrees of freedom, which is the regularized upper incomplete gamma
|
|
// Q(degreesOfFreedom/2, x/2).
|
|
func chiSquareUpperTail(x float64, degreesOfFreedom int) float64 {
|
|
if degreesOfFreedom <= 0 || math.IsNaN(x) {
|
|
return math.NaN()
|
|
}
|
|
if x <= 0 {
|
|
return 1
|
|
}
|
|
return regularizedUpperGamma(float64(degreesOfFreedom)/2, x/2)
|
|
}
|
|
|
|
const (
|
|
gammaIterationLimit = 2000
|
|
gammaTolerance = 1e-15
|
|
gammaTiny = 1e-300
|
|
)
|
|
|
|
// regularizedUpperGamma is Q(shape, x). The series is used below the crossover
|
|
// and the continued fraction above it, as in Numerical Recipes in C, 2nd ed.,
|
|
// section 6.2 (gammp/gammq).
|
|
func regularizedUpperGamma(shape, x float64) float64 {
|
|
if x < shape+1 {
|
|
return 1 - lowerGammaSeries(shape, x)
|
|
}
|
|
return upperGammaContinuedFraction(shape, x)
|
|
}
|
|
|
|
func lowerGammaSeries(shape, x float64) float64 {
|
|
term := 1 / shape
|
|
sum := term
|
|
for iteration := 1; iteration < gammaIterationLimit; iteration++ {
|
|
term *= x / (shape + float64(iteration))
|
|
sum += term
|
|
if math.Abs(term) < math.Abs(sum)*gammaTolerance {
|
|
break
|
|
}
|
|
}
|
|
return sum * math.Exp(-x+shape*math.Log(x)-logGamma(shape))
|
|
}
|
|
|
|
// upperGammaContinuedFraction evaluates Q(shape, x) with the modified Lentz
|
|
// algorithm, Numerical Recipes in C, 2nd ed., section 5.2.
|
|
func upperGammaContinuedFraction(shape, x float64) float64 {
|
|
b := x + 1 - shape
|
|
c := 1 / gammaTiny
|
|
d := 1 / b
|
|
h := d
|
|
for iteration := 1; iteration < gammaIterationLimit; iteration++ {
|
|
numerator := -float64(iteration) * (float64(iteration) - shape)
|
|
b += 2
|
|
d = numerator*d + b
|
|
if math.Abs(d) < gammaTiny {
|
|
d = gammaTiny
|
|
}
|
|
c = b + numerator/c
|
|
if math.Abs(c) < gammaTiny {
|
|
c = gammaTiny
|
|
}
|
|
d = 1 / d
|
|
delta := d * c
|
|
h *= delta
|
|
if math.Abs(delta-1) < gammaTolerance {
|
|
break
|
|
}
|
|
}
|
|
return h * math.Exp(-x+shape*math.Log(x)-logGamma(shape))
|
|
}
|
|
|
|
func logGamma(x float64) float64 {
|
|
value, _ := math.Lgamma(x)
|
|
return value
|
|
}
|