Utility Functions¶

def bh(p_vals: [float], alpha: float) -> (int, float):
    """Static Benjamini-Hochberg (BH) procedure for FDR control.

    The Benjamini-Hochberg procedure is the fundamental method for controlling
    the False Discovery Rate in multiple testing. It provides FDR control at
    level α under independence of p-values or positive dependence (PRDS).

    The algorithm sorts p-values and finds the largest index i such that
    P(i) ≤ (i/n) × α, where P(i) is the i-th smallest p-value and n is the
    total number of tests.

    Args:
        p_vals: List of p-values to test. Should be in [0,1].
        alpha: Target FDR level. Must be in (0,1).

    Returns:
        Tuple of (number_of_rejections, rejection_threshold).
        - number_of_rejections: Number of hypotheses rejected
        - rejection_threshold: The p-value threshold used for rejection

    Examples:
        >>> # Basic usage
        >>> p_values = [0.01, 0.02, 0.5, 0.8]
        >>> num_rej, threshold = bh(p_values, alpha=0.05)
        >>> print(f"Rejected {num_rej} hypotheses at threshold {threshold:.4f}")

        >>> # Check which p-values are rejected
        >>> rejected = [p <= threshold for p in p_values]
        >>> print(f"Rejected p-values: {[p for p, r in zip(p_values, rejected) if r]}")

    References:
        Benjamini, Y., and Y. Hochberg (1995). "Controlling the False Discovery Rate:
        A Practical and Powerful Approach to Multiple Testing." Journal of the Royal
        Statistical Society: Series B, 57(1):289-300.
    """
    n = len(p_vals)
    sorted_p_vals = sorted(p_vals)

    def condition(i):
        return sorted_p_vals[i] <= alpha * (i + 1) / n

    left, right = 0, n
    while left < right:
        mid = (left + right) // 2
        if condition(mid):
            left = mid + 1
        else:
            right = mid

    return left, alpha * left / n if left else 0

def storey_bh(p_vals: [float], alpha: float, lambda_: float) -> (int, float):
    """Static Storey-Benjamini-Hochberg procedure with π₀ estimation.

    The Storey-BH procedure extends the classical Benjamini-Hochberg method by
    estimating π₀ (the proportion of true null hypotheses) and incorporating
    this estimate into the rejection threshold. This typically provides higher
    power when π₀ < 1, which is common in genomic and other high-dimensional
    applications.

    The algorithm estimates π₀ using Storey's method:
    π̂₀ = min(1, (1 + #{p > λ}) / (n × (1 - λ)))

    Then applies the BH procedure with adjusted significance level α/π̂₀.

    Args:
        p_vals: List of p-values to test. Should be in [0,1].
        alpha: Target FDR level. Must be in (0,1).
        lambda_: Threshold parameter for π₀ estimation. Must be in [0,1).
                Common choice: λ = 0.5. Higher values give more conservative
                π₀ estimates.

    Returns:
        Tuple of (number_of_rejections, rejection_threshold).
        - number_of_rejections: Number of hypotheses rejected
        - rejection_threshold: The p-value threshold used for rejection

    Raises:
        ValueError: If lambda_ >= 1.0.

    Examples:
        >>> # Standard usage with λ = 0.5
        >>> p_values = [0.001, 0.02, 0.5, 0.8, 0.9]
        >>> num_rej, threshold = storey_bh(p_values, alpha=0.05, lambda_=0.5)
        >>> print(f"Storey-BH rejected {num_rej} at threshold {threshold:.4f}")

        >>> # Compare with standard BH
        >>> num_rej_bh, threshold_bh = bh(p_values, alpha=0.05)
        >>> print(f"Standard BH rejected {num_rej_bh} at threshold {threshold_bh:.4f}")
        >>> # Storey-BH typically has higher power when π₀ < 1

    References:
        Storey, J. D. (2002). "A direct approach to false discovery rates."
        Journal of the Royal Statistical Society: Series B, 64(3):479-498.

        Storey, J. D. (2003). "The positive false discovery rate: a Bayesian
        interpretation and the q-value." Annals of Statistics, 31(6):2013-2035.
    """
    if not p_vals:
        return 0, 0.0

    n = len(p_vals)

    # Estimate π₀ using Storey's method
    # π₀ = (1 + #{p_i > λ}) / (n(1-λ))
    if lambda_ >= 1.0:
        raise ValueError("lambda_ must be less than 1.0")

    num_above_lambda = sum(1 for p in p_vals if p > lambda_)
    pi0 = min(1.0, (1 + num_above_lambda) / (n * (1 - lambda_)))

    # Apply BH procedure with adjusted alpha
    sorted_p_vals = sorted(p_vals)

    # Find the largest i such that P(i) ≤ (i/n) * (alpha/π₀)
    num_reject = 0
    threshold = 0.0

    for i in range(n):
        if sorted_p_vals[i] <= (i + 1) * alpha / (n * pi0):
            num_reject = i + 1
            threshold = sorted_p_vals[i]
        else:
            break

    return num_reject, threshold

def by(p_vals: [float], alpha: float) -> (int, float):
    """Static Benjamini-Yekutieli (BY) procedure for FDR control under dependence.

    The Benjamini-Yekutieli procedure extends the Benjamini-Hochberg method to
    provide FDR control even under arbitrary dependence among p-values. It uses
    harmonic weights to maintain FDR control at the cost of reduced power
    compared to the standard BH procedure.

    The algorithm finds the largest index i such that:
    P(i) ≤ (i/n) × α / H_n

    where H_n = Σ(1/j) for j=1 to n is the n-th harmonic number.

    Args:
        p_vals: List of p-values to test. Should be in [0,1].
        alpha: Target FDR level. Must be in (0,1).

    Returns:
        Tuple of (number_of_rejections, rejection_threshold).
        - number_of_rejections: Number of hypotheses rejected
        - rejection_threshold: The p-value threshold used for rejection

    Examples:
        >>> # Usage under arbitrary dependence
        >>> p_values = [0.01, 0.02, 0.5, 0.8]  # May be dependent
        >>> num_rej, threshold = by(p_values, alpha=0.05)
        >>> print(f"BY rejected {num_rej} at threshold {threshold:.4f}")

        >>> # Compare with BH (more conservative under dependence)
        >>> num_rej_bh, threshold_bh = bh(p_values, alpha=0.05)
        >>> print(f"BH rejected {num_rej_bh} (may not control FDR under dependence)")

    Notes:
        The BY procedure is recommended when:
        - P-values may be arbitrarily dependent
        - Conservative FDR control is required
        - The exact dependence structure is unknown

        It maintains valid FDR control under any dependence structure but
        typically has lower power than BH under independence.

    References:
        Benjamini, Y., and D. Yekutieli (2001). "The control of the false discovery
        rate in multiple testing under dependency." Annals of Statistics, 29(4):1165-1188.
    """
    n = len(p_vals)
    sorted_p_vals = sorted(p_vals)
    harmonic_sum = sum(1 / (i + 1) for i in range(n))

    def condition(i):
        return sorted_p_vals[i] <= alpha * (i + 1) / (n * harmonic_sum)

    left, right = 0, n
    while left < right:
        mid = (left + right) // 2
        if condition(mid):
            left = mid + 1
        else:
            right = mid

    return left, alpha * left / (n * harmonic_sum) if left else 0

def check_p_val(p_val: float) -> None:
    """Validate that a p-value is in the valid range [0, 1].

    Args:
        p_val: The p-value to validate.

    Raises:
        ValueError: If p_val is not in [0, 1].

    Examples:
        >>> check_p_val(0.05)  # Valid - no exception
        >>> check_p_val(1.5)   # Raises ValueError
        Traceback (most recent call last):
        ValueError: Given p-value must be between [0,1].
    """
    if not 0 <= p_val <= 1:
        raise ValueError(
            """
            Given p-value must be between [0,1].
            """
        )

def check_alpha(p_val: float) -> None:
    """Validate that an alpha value is in the valid range (0, 1).

    Args:
        p_val: The alpha/significance level to validate.

    Raises:
        ValueError: If p_val is not in (0, 1).

    Examples:
        >>> check_alpha(0.05)  # Valid - no exception
        >>> check_alpha(0.0)   # Raises ValueError (not > 0)
        >>> check_alpha(1.0)   # Raises ValueError (not < 1)
    """
    if not 0 < p_val < 1:
        raise ValueError(
            """
            Given alpha value must be between (0,1).
            """
        )

def check_initial_wealth(initial_wealth: float, alpha: float) -> None:
    """Validate that initial wealth is properly bounded relative to alpha.

    In alpha-wealth procedures (LORD, SAFFRON, etc.), the initial wealth
    must be positive but less than the target FDR level to ensure proper
    wealth dynamics and FDR control.

    Args:
        initial_wealth: The initial alpha wealth.
        alpha: The target FDR level.

    Raises:
        ValueError: If initial_wealth is not in (0, alpha).

    Examples:
        >>> check_initial_wealth(0.025, 0.05)  # Valid - no exception
        >>> check_initial_wealth(0.1, 0.05)    # Raises ValueError (> alpha)
        >>> check_initial_wealth(0.0, 0.05)    # Raises ValueError (not > 0)
    """
    if not 0 < initial_wealth < alpha:
        raise ValueError(
            """
            The initial wealth should be between (0, alpha).
            """
        )

def check_candidate_threshold(lambda_: float) -> None:
    """Validate that a candidate threshold λ is in the valid range (0, 1).

    The candidate threshold λ is used in SAFFRON and ADDIS algorithms to
    determine which p-values are considered "candidates" for rejection.
    It must be strictly between 0 and 1.

    Args:
        lambda_: The candidate threshold to validate.

    Raises:
        ValueError: If lambda_ is not in (0, 1).

    Examples:
        >>> check_candidate_threshold(0.5)  # Valid - no exception  
        >>> check_candidate_threshold(0.0)  # Raises ValueError
        >>> check_candidate_threshold(1.0)  # Raises ValueError
    """
    if not 0 < lambda_ < 1:
        raise ValueError(
            """
            The candidate threshold (lambda) should be between (0, 1).
            """
        )

def check_wealth(wealth: float) -> None:
    """Validate that alpha wealth is positive (not depleted).

    In alpha-wealth procedures, when wealth reaches zero or below,
    no more tests can be performed while maintaining FDR control.
    This function checks for wealth depletion.

    Args:
        wealth: Current alpha wealth level.

    Raises:
        ValueError: If wealth <= 0.

    Examples:
        >>> check_wealth(0.01)  # Valid - no exception
        >>> check_wealth(0.0)   # Raises ValueError (wealth depleted)
        >>> check_wealth(-0.1)  # Raises ValueError (negative wealth)
    """
    if wealth <= 0:
        raise ValueError(
            """
            Alpha wealth depleted. Test execution stopped.
            """
        )

def check_decay_factor(decay_factor: float) -> None:
    """Validate that a decay factor is in the valid range (0, 1).

    Decay factors are used in various algorithms to control the rate
    at which parameters decrease over time. They must be strictly
    between 0 and 1 for proper convergence behavior.

    Args:
        decay_factor: The decay factor to validate.

    Raises:
        ValueError: If decay_factor is not in (0, 1).

    Examples:
        >>> check_decay_factor(0.9)   # Valid - no exception
        >>> check_decay_factor(0.0)   # Raises ValueError (not > 0)
        >>> check_decay_factor(1.0)   # Raises ValueError (not < 1)
        >>> check_decay_factor(1.5)   # Raises ValueError (not < 1)
    """
    if not 0 < decay_factor < 1:
        raise ValueError(
            """
            Decay factor must be between (0,1).
            """
        )