andrey.core.ci_cmiknn.CMIknn

class andrey.core.ci_cmiknn.CMIknn(data, *, k_cmi=0.1, k_perm=5, n_permutations=500, seed=None)

Conditional mutual information test with a nearest-neighbor estimator (CMIknn).

Each query estimates the conditional mutual information of x and y given the conditioning set from nearest-neighbor distances, then compares it with estimates on copies of the data where x is shuffled only among rows close in the conditioning set. That keeps the dependence of x on the conditioning set and removes any other.

Parameters:

• data (ndarray) – Observed data, one row per sample. Each column is replaced by its ranks.

• k_cmi (float) – Nearest neighbors of the estimator: a fraction of the rows in (0, 1), or a count, an int from 1 to one less than the rows. More neighbors lower the estimate’s variance and raise its bias toward zero, which costs power as the conditioning set grows.

• k_perm (int) – Nearest neighbors in the conditioning set among which x is shuffled.

• n_permutations (int) – Shuffled copies that form the null distribution.

• seed (int | None) – Seed for the tie-breaking and the shuffles. None uses ANDREY_SEED, or 0 when it is unset.

Raises:

ValueError – If data is not 2-D, has fewer than 2 rows, or contains NaN or infinity, k_cmi is neither a fraction in (0, 1) nor an int from 1 to one less than the rows, or k_perm or n_permutations is below 1.

Notes

How it works. Each column is replaced by its ranks, ties ordered at random, and distances are in the maximum norm. With k the neighbor count (round(k_cmi n) for a fraction), eps_i the distance from row i to its k-th nearest row in the space of x, y, and the conditioning set z, and k_xz, k_yz, k_z the rows strictly closer than eps_i in each subspace (row i counted), the estimate is I = psi(k) + mean_i [psi(k_z) - psi(k_xz) - psi(k_yz)], psi the digamma function [1]. Each null copy maps every row to one of its k_perm nearest rows in z, avoiding a row already taken where it can. With b of the B copies estimating at least I, the p-value is (b + 1) / (B + 1): the observed data count as one more draw, so the p-value is never zero [3]. In the skeleton searches of PC, FCI, GFCI, and CD-NOD, a query’s shuffles stop once p > alpha or p <= alpha is certain whatever the remaining shuffles give, so the search decides as with every shuffle; where the search orients colliders by p-value (collider_order="p_value"), only p <= alpha stops them, so every p-value it ranks is the full one. An empty conditioning set shuffles x over all rows, and the estimate becomes the mutual information estimate of x and y [2]. The default k_cmi, a tenth of the rows, is the low end of the range [1] recommends; a larger share smooths the estimate, which costs power as the conditioning set grows.

statistic() returns the estimate I itself, in nats, for use as a dependence measure. It is biased on finite samples, so compare estimates at the same rows and k: on Gaussian pairs with 1,000 rows and k = 5 it reads 0.01 to 0.05 above the true value, independent pairs included (0.03 to 0.09 with 200 rows), and with k a tenth of the rows it reads near zero for independent pairs and low for strong dependence, 0.06 below at a correlation of 0.9 and 0.26 below given one more variable. The p-value compares I with shuffled copies, which share the bias.

Of x and y, the one with the lower tie_ranks() rank, a rank computed from each column’s own data, is shuffled. The tie-breaking is drawn once, from seed, when the test is built, and each column takes its keys by that rank. The shuffles come from a generator seeded by seed and the ranks of the conditioning set, so every query with the same conditioning set uses the same shuffles, and the test object keeps them for later queries. A p-value does not depend on the order of the queries, of the arguments, or of the columns.

The distances and neighbor counts run on a CUDA GPU under the cuda backend, or under auto from ANDREY_CMIKNN_GPU_THRESHOLD rows squared, a batch of shuffles at a time. The counts are integers, so the p-value equals the CPU’s.

Limitations. Measured in PC on simulated data from sparse random graphs (mean degree 2), each variable a weighted sum of z + sin(z), or of cos(z), over its parents plus Gaussian noise, with 20 datasets per size on 4 cores. At 5 variables with 500 to 2,000 rows, 10 with 500 and 1,000, and 20 with 500, its median errors (SHD) stayed below the empty graph’s at every size. The median fit took 47 seconds at 5 variables and 500 rows, and 12 to 14 minutes at 5 variables with 2,000 rows, 10 with 1,000, and 20 with 500. Each query computes n_permutations + 1 estimates, each over n x n distance matrices, so its time grows with the square of the number of rows times the permutations, and its memory with the square of the rows: hundreds of rows take seconds per query, and tens of thousands do not fit. The p-value moves with seed, is a multiple of 1 / (n_permutations + 1), and is at least that.

References

[1] (1,2)

Runge, J. (2018). Conditional independence testing based on a nearest-neighbor estimator of conditional mutual information. AISTATS, PMLR 84, 938-947.

[2]

Kraskov, A., Stogbauer, H., and Grassberger, P. (2004). Estimating mutual information. Physical Review E 69, 066138.

[3]

Phipson, B., and Smyth, G. K. (2010). Permutation p-values should never be zero: calculating exact p-values when permutations are randomly drawn. Statistical Applications in Genetics and Molecular Biology 9(1), Article 39.

statistic(x, y, condition_set=())

Return the estimate of the conditional mutual information of x and y given condition_set, in nats.

The estimate the p-value compares with its shuffles, from the test’s ranks (see Notes); with an empty condition_set, the mutual information estimate. It does not depend on the order of the arguments.

Return type:

float