LOND: Levels based On Number of Discoveries¶

from online_fdr.investing.lond.lond import Lond

# Create LOND instance
lond = Lond(alpha=0.05)

# Test individual p-values
p_values = [0.001, 0.15, 0.03, 0.8, 0.02, 0.45, 0.006]

print("LOND Online Testing:")
discoveries = []

for i, p_value in enumerate(p_values):
    decision = lond.test_one(p_value)

    if decision:
        discoveries.append(i + 1)
        print(f"✓ Test {i+1}: p={p_value:.3f} → DISCOVERY! (total: {lond.num_reject})")
    else:
        print(f"  Test {i+1}: p={p_value:.3f} → no rejection (threshold: {lond.alpha:.6f})")

print(f"\nTotal discoveries: {len(discoveries)}")
print(f"Discovery indices: {discoveries}")

def demonstrate_discovery_momentum():
    """Show how LOND's power depends on early discoveries."""

    # Scenario 1: Early discoveries
    print("Scenario 1: Early Discoveries")
    print("=" * 35)

    lond1 = Lond(alpha=0.1)  # Higher alpha for better visibility
    early_discoveries = [0.001, 0.005, 0.02, 0.8, 0.3, 0.04, 0.2]

    for i, p_val in enumerate(early_discoveries, 1):
        decision = lond1.test_one(p_val)
        print(f"Test {i}: p={p_val:.3f} → {'REJECT' if decision else 'ACCEPT'} "
              f"(threshold: {lond1.alpha:.6f}, discoveries: {lond1.num_reject})")

    print(f"Final discoveries: {lond1.num_reject}\n")

    # Scenario 2: No early discoveries  
    print("Scenario 2: No Early Discoveries")
    print("=" * 37)

    lond2 = Lond(alpha=0.1)
    no_early = [0.8, 0.9, 0.7, 0.001, 0.005, 0.02, 0.04]  # Same p-values, reordered

    for i, p_val in enumerate(no_early, 1):
        decision = lond2.test_one(p_val)
        print(f"Test {i}: p={p_val:.3f} → {'REJECT' if decision else 'ACCEPT'} "
              f"(threshold: {lond2.alpha:.6f}, discoveries: {lond2.num_reject})")

    print(f"Final discoveries: {lond2.num_reject}")
    print(f"\nPower difference due to ordering: {lond1.num_reject - lond2.num_reject} discoveries")

demonstrate_discovery_momentum()

def lond_with_dependence():
    """Compare LOND variants for different dependence assumptions."""

    print("LOND Dependence Handling:")
    print("=" * 30)

    # Test p-values with some correlation structure
    p_values = [0.01, 0.02, 0.015, 0.8, 0.03, 0.9, 0.025, 0.7]

    # Independent version
    lond_indep = Lond(alpha=0.05, dependent=False)
    indep_decisions = [lond_indep.test_one(p) for p in p_values]

    # Dependent version (with harmonic correction)
    lond_dep = Lond(alpha=0.05, dependent=True)
    dep_decisions = [lond_dep.test_one(p) for p in p_values]

    print("Results:")
    print(f"Independent assumption: {sum(indep_decisions)} discoveries")
    print(f"Dependent correction: {sum(dep_decisions)} discoveries")
    print(f"Dependent version is more conservative as expected")

    # Show threshold differences
    print(f"\nFinal thresholds:")
    print(f"Independent: {lond_indep.alpha:.6f}")  
    print(f"Dependent: {lond_dep.alpha:.6f}")

lond_with_dependence()

def compare_lond_variants():
    """Compare original and modified LOND formulations."""

    print("LOND Formulation Comparison:")
    print("=" * 32)

    p_values = [0.002, 0.8, 0.01, 0.9, 0.005, 0.7, 0.015]

    # Original: uses (num_reject + 1) in threshold
    lond_orig = Lond(alpha=0.05, original=True)

    # Modified: uses max(num_reject, 1) in threshold
    lond_mod = Lond(alpha=0.05, original=False)

    print("Test | P-value | Original | Modified")
    print("-" * 35)

    for i, p_val in enumerate(p_values, 1):
        orig_decision = lond_orig.test_one(p_val)
        mod_decision = lond_mod.test_one(p_val)

        print(f"{i:4d} | {p_val:7.3f} | {'REJECT' if orig_decision else 'ACCEPT':>8} | "
              f"{'REJECT' if mod_decision else 'ACCEPT':>8}")

    print(f"\nOriginal formulation: {lond_orig.num_reject} discoveries")
    print(f"Modified formulation: {lond_mod.num_reject} discoveries")

compare_lond_variants()

from online_fdr.utils.generation import DataGenerator, GaussianLocationModel

def evaluate_lond_performance():
    """Evaluate LOND on simulated data."""

    print("LOND Performance Evaluation:")
    print("=" * 32)

    # Generate realistic test scenario
    dgp = GaussianLocationModel(alt_mean=2.0, alt_std=1.0, one_sided=True)
    generator = DataGenerator(n=100, pi0=0.9, dgp=dgp)  # 90% nulls

    # Create LOND instance
    lond = Lond(alpha=0.05)

    # Simulate testing
    true_positives = 0
    false_positives = 0
    total_tests = 50

    print(f"Testing {total_tests} hypotheses (10% alternatives expected):")
    print()

    for i in range(total_tests):
        p_value, is_alternative = generator.sample_one()
        decision = lond.test_one(p_value)

        if decision:
            if is_alternative:
                true_positives += 1
                result = "TRUE discovery ✓"
            else:
                false_positives += 1
                result = "FALSE discovery ✗"

            truth = "ALT" if is_alternative else "NULL"
            print(f"Test {i+1:2d}: p={p_value:.3f} ({truth}) → REJECT ({result})")

    # Calculate metrics
    total_discoveries = true_positives + false_positives
    empirical_fdr = false_positives / max(total_discoveries, 1)

    print(f"\nPerformance Summary:")
    print(f"Total discoveries: {total_discoveries}")
    print(f"True positives: {true_positives}")
    print(f"False positives: {false_positives}")
    print(f"Empirical FDR: {empirical_fdr:.3f}")
    print(f"Target FDR: {lond.alpha0}")
    print(f"FDR controlled: {'✓' if empirical_fdr <= lond.alpha0 else '✗'}")

evaluate_lond_performance()