ENH: Implement the exact p-value of the Cramér-von Mises two-sample test

[fbourgey]

Reference issue

Toward #98

What does this implement/fix?

Implement the exact p-value of the Cramer-von Mises two-sample test for a given value of the test statistic from _pval_cvm_2samp_exact

https://en.wikipedia.org/wiki/Cram%C3%A9r%E2%80%93von_Mises_criterion

[mdhaber]

Can you provide additional information that helps the reader compare this to either the original code, which is pretty hard to follow, or the algorithm, which looks simple but has terms like $I[a^2u(u+1)(2u + 1)/6]$ and $I[b^2u(u+1)(2u + 1)/6]$ that don't obviously appear in either version? If it would be cleaner, feel free to express this in terms of Algorithm 1 from the paper rather than the Python code. There are a lot of aspects of the Python that don't look like a very efficient expression of the algorithm (lots of placeholding size-zero array, use of arrays instead of dicts, use of delete, etc.), so if it can't be translated directly, it is probably not worth trying to maintain a resemblance.

[fbourgey]

I agree; I've found the Python code quite difficult to follow. I’ll provide more detail and include screenshots of the relevant sections of the paper in the coming days.

[steppi]

@mdhaber, you want this to be able to work in CuPy right? std::map and std::vector aren't available in CUDA.

Also, rather than putting this whole thing into a scalar kernel, and making a ufunc out it, it seems like it may be better to take the idea of the original function and decompose it into the ufuncs that are needed to make it work. I think the key component here would be a gufunc to generate the frequency table. lcm and the binomial coefficients could be handled through delegation, and the rest looks like it could be made array-agnostic.

[mdhaber]

you want this to be able to work in CuPy right?

Yes. OK, good to know.

[fbourgey]

Page 3 in [1]: Definition of the variable lcm. This is $L$ in [1].

[fbourgey]

Page 4 in [1]: definitions of $a$ and $b$

[fbourgey]

Using consistent notations and following the Python comment:
$s$ is $U$ in [2] and $T_2$ in [1] same as $T$ in [2]

From [1, page 3]

Rewriting in terms of $\zeta$, we have

$$ \zeta = \frac{(n+n)^2 L^2}{mn} T_2 $$

and from [2, eq. 9] (and using the same notations as in [1])

$$ T_2 = \frac{s}{nm(n+m)} - \frac{4mn-1}{6(m+n)} $$

Simple algebra gives

$$ \zeta = \frac{(m+n)L^2}{6 m^2 n^2} (6 s - (4mn - 1) mn) $$

We then floor as $\zeta$ is an integer-valued statistic.

[fbourgey]

Initialize the frequency tables of $\zeta$ for paths reaching $(u,v)$.

Each table is a $2 \times k$ array where $k =$ number of distinct values of $\zeta$ at $(u,v)$ ie

g = [[zeta_1, zeta_2, ..., zeta_k],
     [count_1, count_2, ..., count_k]]

row 0 = distinct values of $\zeta$ and row 1 = how many paths give each value

Initially: $g_{0,0}$ = {(0,1)}

[fbourgey]

Compute p-value

$$ \mathbb{P}(\zeta \geq \zeta_{obs}) = \frac{\text{number of paths with } \zeta \geq \zeta_{obs}}{\binom{m+n}{m}} $$

[fbourgey]

Apply shift $S_{d^2}$. From Eq. 11 of [1]

[fbourgey]

Translation of the original Python code.

    vi, i0, i1 = np.intersect1d(tmp[0], g[0], return_indices=True)
    tmp = np.concatenate([
        np.stack([vi, tmp[1, i0] + g[1, i1]]),
        np.delete(tmp, i0, 1),
        np.delete(g, i1, 1)
    ], 1)

• First line: find common $\zeta$ values between $g_{u,v-1}$ and $g_{u-1,v}$
• Second line: this implements $g_{u,v-1} + g_{u-1,v}$

[fbourgey]

It appears that (12) is not needed, as this is produced automatically by how gs and tmp are initialized.

For $v=0$, we have $d^2 = (av-bu)^2 = b^2 u^2$. Hence,

$$ g_{u,0} = S_{b^2 u^2} (g_{u-1,0}). $$

Starting from $g_{0,0} = I[0]$, at step $u=1$, $g_{1,0} = I[b^2]$, ... and for a general,

$$ g_{u,0} = I\left[\sum_{j=1}^{u} (bj)^2\right] = b^2 \frac{u(u+1)(2u+1)}{6} $$

Here, we use that $S_d(I[c]) = I[c+d])$