Blum's test for independence

blums.py

def blums(x, y=None, sig_lev=0.01):
    """
    Performs Blum's test for independence.
    
    Blum, J. R., et al. “Distribution Free Tests of Independence Based on the 
    Sample Distribution Function.” The Annals of Mathematical Statistics, vol. 
    32, no. 2, 1961, pp. 485–498. JSTOR, www.jstor.org/stable/2237758.

    Parameters
    ----------
    x : numpy.ndarray
        First vector
    y : numpy.ndarray, optional
        Second vector. If None, use lag-1 of first vector.
    sig_lev : TYPE, optional
        Significance level; must be [0.1, 0.05, 0.02, 0.01, 0.005]. The default is 0.01.

    Returns
    -------
    h0 : bool
        Null hypothesis of independence
    B : float
        Blum's test statistic
    Bcr : float
        Critical value for B at the specified significace level

    """
    
    if y is None:
        y = x[1:]
        x = x[:-1]
    
    assert isinstance(x, np.ndarray)
    assert isinstance(y, np.ndarray)
    assert len(x) == len(y)
    
    Bcrs = {0.1: 2.286,
           0.05: 2.2844,
           0.02: 3.622,
           0.01: 4.230,
           0.005: 4.851}
    assert sig_lev in list(Bcrs.keys()), 'argumemt ''sig_lev'' must be one of following: [0.1, 0.05, 0.02, 0.01, 0.005], received {:}'.format(sig_lev)
    Bcr = Bcrs[sig_lev]
    
    N = len(x)
    N1 = np.zeros_like(x)
    N2 = np.zeros_like(x)
    N3 = np.zeros_like(x)
    N4 = np.zeros_like(x)
    for idx, (xi, yi) in enumerate(zip(x,y)):
        N1[idx] = sum([xj <= xi and yj <= yi for (xj,yj) in zip(x,y)])
        N2[idx] = sum([xj > xi and yj <= yi for (xj,yj) in zip(x,y)])
        N3[idx] = sum([xj <= xi and yj > yi for (xj,yj) in zip(x,y)])
        N4[idx] = sum([xj > xi and yj > yi for (xj,yj) in zip(x,y)])
    
    B = 0.5*np.pi**4*N*np.sum((N1*N4 - N2*N3)**2/N**5)
    if B > Bcr:
        h0 = False
    else:
        h0 = True
    return h0, B, Bcr