Modified Mandel-Betensky simultaneous nonparametric confidence interval

modified_mandel_betensky.py

"""
Implements the nonparametric simultaneous confidence interval described
in Gao et al. [1], which modifies the Mandel-Betensky simultaneous CI [2]
such that it is admissible. Notation in the code below comes from [1].

Takes as input an array of bootstrap samples / replicates, returns lower
and upper bounds at 1-alpha confidence.

[1]: https://www.mdpi.com/2073-8994/13/7/1212
[2]: https://pubmed.ncbi.nlm.nih.gov/26452746/
"""
import numpy as np

def mmb_simul_ci(boots, alpha=0.05):
    # boots is n_boot x n_est
    B, p = boots.shape
    
    # step1: ranks over boot dim
    ranks = np.argsort(boots, axis=0)
    assert ranks.shape == boots.shape
    b_tilde = np.sort(boots, axis=0)
    max_ranks = ranks.max(axis=1)
    assert max_ranks.shape == (B,)
    
    # step2: 1-a/2 quantile
    upper_q = np.quantile(max_ranks, 1 - alpha / 2, interpolation="higher")
    
    # step3: cull outside upper to make Phi
    Phi = max_ranks <= upper_q
    
    # step4: lower ranks
    boots_Phi = boots[Phi]
    Phi_ranks = np.argsort(boots_Phi, axis=0)
    assert Phi_ranks.shape == boots_Phi.shape
    min_ranks = Phi_ranks.min(axis=1)
    assert min_ranks.shape == (Phi.sum(),)
                      
    # step5: lower end
    lower_q = np.quantile(min_ranks, alpha / (2 - alpha), interpolation="lower")
    Psi = Phi.copy()
    Psi[Phi][min_ranks < lower_q] = False
    
    # step6: get the ci lims
    ranks_Psi = ranks[Psi, :]
    arange_p = np.arange(p)
    w = ranks_Psi.min(axis=0)
    low = b_tilde[w, arange_p]
    t = ranks_Psi.max(axis=0)
    high = b_tilde[t, arange_p]
    
    return low, high