[Parva-26]

Follow-up to #14 and #48. In #14, @jefish003 preferred adding very small noise to break ties. This PR does that and finishes the Chebyshev default for the cases #48 didn't cover.

1. Tie-breaking noise

When a point has k or more exact copies, the distance to its k-th neighbour is 0, the strict < counts give -1 and digamma(0) is -inf. knn_mutual_information then returned inf and knn_conditional_mutual_information returned nan, for example on Poisson(3) counts with k=3 or on bootstrap resamples. A new helper, _add_tie_breaking_noise, adds independent noise of relative size 1e-10 to X, Y and Z before the neighbour search, as Tigramite's CMIknn and scikit-learn's mutual_info_regression do. The seed is fixed, so results stay deterministic. On continuous data the estimates don't change (the numbers in [1] below are identical with and without the noise).

2. Chebyshev when Z is None

With metric=None, the Z=None branch still used Minkowski (p=2) and knn_mutual_information defaulted to Euclidean, so the marginal case kept the Euclidean bias: with N=1000 and a true MI of 0.144, the mean estimate was -0.118 for dependent data and -0.264 for independent data, against +0.146 and -0.001 with Chebyshev (script below). Both defaults are now Chebyshev. On a small kNN discovery run (n=4, T=300, 6 seeds) this made little difference to the graph (TP/FP/FN 20/1/5 before, 19/1/6 with Chebyshev), but direct MI calls returned the wrong sign.

3. Docstring

The knn_conditional_mutual_information docstring gave the two incorrect identities noted in #14 and described a decomposition the code doesn't use. It now gives the Frenzel-Pompe / Vejmelka-Paluš formula that the code implements, with references.

4. Bootstrap intervals with k-NN estimators

bootstrap_cmi_confidence_interval(method="knn") used to raise a ValueError because of NaN replicates (see #40). With the tie fix it runs, but repeated rows in the resamples still bias the estimates: over 20 datasets the 95% interval covered the true CMI 7/20 times for knn, 13/20 for geometric_knn and 19/20 for gaussian. The function now emits a UserWarning for knn and geometric_knn. Its metric default is now None, so method="knn" uses Chebyshev like discover_network.

Tests

New causationentropy/tests/core/information/test_knn_ties_and_defaults.py: finite estimates with repeated values and count data, determinism, inputs not modified, the Chebyshev default and its accuracy for Z=None, and the bootstrap warning. 6 of the 8 tests fail on current main.

Validation

• Default suite: 471 passed, 17 deselected.
• black --check, isort --check-only and the flake8 hard gate pass.

Not changed: the discover_network docstring still says metric defaults to 'euclidean' (out of date since #48), and the internal functions in discovery.py keep metric="euclidean" defaults. I left discovery.py alone to avoid conflicts with #33 and #46.

Script for the numbers above (run from the repo root of this branch; part [2] takes about 4 minutes)

# Numbers for the k-NN follow-up. Run from the repo root of this branch.
# Part [2] takes about 4 minutes because of geometric_knn.
import warnings

import numpy as np

from causationentropy.core.information.conditional_mutual_information import (
    knn_conditional_mutual_information,
)
from causationentropy.core.stats import bootstrap_cmi_confidence_interval

warnings.filterwarnings("ignore")

print("[1] Z=None, k=5, 100 datasets per row, Y = 0.5 X + noise (true MI 0.144)")
rng = np.random.default_rng(1)
r = 0.5
for N in (300, 1000):
    res = {"minkowski": ([], []), None: ([], [])}
    for _ in range(100):
        X = rng.standard_normal((N, 1))
        Y_indep = rng.standard_normal((N, 1))
        Y_dep = r * X + np.sqrt(1 - r**2) * rng.standard_normal((N, 1))
        for m in res:
            res[m][0].append(knn_conditional_mutual_information(X, Y_indep, None, metric=m, k=5))
            res[m][1].append(knn_conditional_mutual_information(X, Y_dep, None, metric=m, k=5))
    for m, (indep, dep) in res.items():
        label = "old default (minkowski)" if m else "new default (chebyshev)"
        print(f"    N={N} {label}: independent {np.mean(indep):+.3f} | dependent {np.mean(dep):+.3f}")

print("[2] bootstrap_cmi_confidence_interval, 20 datasets, N=300, true CMI 0.154")
truth = 0.5 * np.log(1 + 0.6**2)
for method in ("knn", "geometric_knn", "gaussian"):
    cover = 0
    for s in range(20):
        rng = np.random.default_rng(100 + s)
        X = rng.standard_normal((300, 1))
        Z = rng.standard_normal((300, 1))
        Y = 0.6 * X + 0.8 * Z + rng.standard_normal((300, 1))
        _, lo, hi, _ = bootstrap_cmi_confidence_interval(
            X, Y, Z, method=method, k=5, n_bootstraps=100, seed=s
        )
        cover += lo <= truth <= hi
    print(f"    {method:14s} 95% interval covered the true value in {cover}/20 datasets")