Cucconi.py

import numpy as np
import scipy.stats as stats

def Cucconi (x, y):


  # Calculates the test statistic for the Cucconi two-sample location-scale test
  # assume x is smallest, swap as needed
  
  m = len(x)
  n = len(y)
  N = m + n
  S = stats.rankdata(np.concatenate(x, y))[(m + 1):N]
  denom = sqrt(m * n * (N + 1) * (2 * N + 1) * (8 * N + 11) / 5)

  U = (6 * np.sum(S^2) - n * (N + 1) * (2 * N + 1)) / denom
  V = (6 * np.sum(( (N + 1) - S)^2) - n * (N + 1) * (2 * N + 1)) / denom
  rho = (2 * (N^2 - 4)) / ((2 * N + 1) * (8 * N + 11)) - 1
  C = (U^2 + V^2 - 2 * rho * U * V) / (2 * (1 - rho^2))

  return C



# Different pvalue, depending on size n,m
# start of range is 6,6 end is inf,inf
# use midpoint ?

C > 2.24 - 2.30  |  p-val ~10%   |  2.27
C > 2.71 - 2.99  |  p-val ~5%    |  2.85
C > 3.62 - 4.60  |  p-val ~1%    |  4.11

todo fit surface ?
use 1000, 10000, 100000 as infinite with same value