ENH: Add exact Cramér-von Mises two-sample p-value gufunc#118
ENH: Add exact Cramér-von Mises two-sample p-value gufunc#118fbourgey wants to merge 6 commits intoscipy:mainfrom
Conversation
…ue calculation for Cramér-von Mises two-sample test
| test_case{8862.0, 14, 8, 0.2679738562091503, 1e-10}, | ||
| test_case{3491.0000000000005, 14, 5, 0.34657722738218094, 1e-10}, | ||
| test_case{12559.0, 5, 26, 0.11812654860485784, 1e-10}, | ||
| test_case{8901.0, 23, 5, 0.9907610907610908, 1e-10} //, test_case{119376.0, 20, 21, 0.5716351061359124, 1e-10} |
There was a problem hiding this comment.
When m=20 and n=21, the test was extremely slow, likely due to it taking too much memory.
| // Device-safe greatest common divisor (gcd) for 64-bit integers | ||
| XSF_HOST_DEVICE inline int64_t gcd(int64_t a, int64_t b) { | ||
| a = (a < 0) ? -a : a; | ||
| b = (b < 0) ? -b : b; | ||
| while (b != 0) { | ||
| int64_t t = a % b; | ||
| a = b; | ||
| b = t; | ||
| } | ||
| return a; | ||
| } | ||
|
|
||
| // Device-safe least common multiple (lcm) for 64-bit integers | ||
| XSF_HOST_DEVICE inline int64_t lcm(int64_t a, int64_t b) { | ||
| if (a == 0 || b == 0) { | ||
| return 0; | ||
| } | ||
| int64_t g = gcd(a, b); | ||
| int64_t res = (a / g) * b; | ||
| return (res < 0) ? -res : res; | ||
| } |
There was a problem hiding this comment.
It does not seem like we can call std::lcm anymore. This implements gcd and lcm. We can move those somewhere else so other methods can use them in the future.
There was a problem hiding this comment.
It turns out thatcuda::std::numeric has it, so we're good on this front actually and don't need to reimplement.
| int64_t zeta = | ||
| static_cast<int64_t>(std::floor((lcm * lcm * (m + n) * (6.0 * s - mn * (4.0 * mn - 1))) / (6.0 * mn * mn))); | ||
|
|
||
| detail::cvm_freq_table_all(m, n, a, b, gs, next_gs); |
There was a problem hiding this comment.
This isn't the right idea. cramer_von_mises_exact should be taking in an already generated freq_table not generating a new one from scratch inside the kernel. Table generation has no s dependence, and we don't want to have to regenerate the table for each scalar value of s.
| using cuda::std::gcd; | ||
| using cuda::std::lcm; | ||
|
|
There was a problem hiding this comment.
Don't worry about the CUDA side. I still have to figure out cupy/cupy#9839 before we can even try these in CuPy. I'm pretty sure just adding using cuda::std::gcd won't work in all cases, and we actually need wrappers for stdlib functions like the other ones in this file. I recall I had suggested using using like this when Irwin first set this up, but there was a reason he had done things the way he did.
2 participants
Supersedes #108
#108 is a direct a C++ port of SciPy's
_pval_cvm_2samp_exactand contains a detailed review explaining the algorithm.The current PR tries to implement @steppi's comment for computing the exact p-value of the Cramér-von Mises two-sample test and making it CUDA compatible.