Add symbolic (partition-based) entropy, MI and CMI estimators (#59)

CHANGELOG.md

@@ -23,6 +23,7 @@
 - `linear_gaussian_from_graph` for generating synthetic data from a chosen `networkx` graph (#39).
 - Block-bootstrap confidence intervals: `moving_block_bootstrap_indices`, `stationary_bootstrap_indices`, `bootstrap_confidence_interval` and `bootstrap_cmi_confidence_interval` (#40).
 - Optional futility early stopping in `shuffle_test` (`early_stop=True`). It is only available when calling `shuffle_test` directly (#41).
+- Symbolic (partition-based) estimators following Cafaro, Lord, Sun and Bollt (Chaos, 2015): `symbolize`, `symbolic_entropy`, `symbolic_mutual_information` and `symbolic_conditional_mutual_information`, available as `method='symbolic'` in `conditional_mutual_information` and `information='symbolic'` in `discover_network` (#59).
 - Multiple-testing corrections `bonferroni_correction`, `benjamini_hochberg_correction`, `benjamini_yekutieli_correction`, `adaptive_bh_correction` and `estimate_null_proportion`, plus `apply_test_correction` to correct the p-values of a discovered network after discovery (adds a `P_Adjusted` column to `network_to_dataframe`). Discovery itself does not apply a correction (#37).
 - Tutorials and expanded theory documentation (#30, #49, #50).
 

README.md

@@ -95,7 +95,7 @@ from causationentropy import discover_network
 network = discover_network(
     data,
     method='standard',          # 'standard', 'alternative', base lines: 'information_lasso', or 'lasso'
-    information='gaussian',     # 'gaussian', 'knn', 'kde', 'geometric_knn', or 'poisson'
+    information='gaussian',     # 'gaussian', 'knn', 'kde', 'geometric_knn', 'poisson', or 'symbolic'
     max_lag=5,                  # Maximum time lag to consider
     alpha_forward=0.05,         # Forward selection significance
     alpha_backward=0.05,        # Backward elimination significance

causationentropy/core/discovery.py

@@ -84,6 +84,8 @@ def discover_network(
         - 'kde': Kernel density estimation
         - 'geometric_knn': Geometric mean k-NN estimator
         - 'poisson': Poisson distribution assumption
+        - 'symbolic': Plug-in estimate on symbolized data (Cafaro et al., 2015);
+          each variable is split near its median (two symbols)
     max_lag : int, default=5
         Maximum time lag to consider in causal relationships. The algorithm examines
         lags from 1 to max_lag (inclusive).
@@ -176,7 +178,14 @@ def discover_network(
 
     if method not in ["standard", "alternative", "information_lasso", "lasso"]:
         raise NotImplementedError(f"discover_network: method={method} not supported.")
-    supported_information_types = ["gaussian", "knn", "kde", "geometric_knn", "poisson"]
+    supported_information_types = [
+        "gaussian",
+        "knn",
+        "kde",
+        "geometric_knn",
+        "poisson",
+        "symbolic",
+    ]
     if information not in supported_information_types:
         raise NotImplementedError(
             f"discover_network: information={information} not supported. "

causationentropy/core/information/conditional_mutual_information.py

@@ -1,3 +1,5 @@
+from typing import Optional
+
 import numpy as np
 from scipy.spatial.distance import cdist
 from scipy.special import digamma
@@ -7,12 +9,15 @@
     kde_entropy,
     poisson_entropy,
     poisson_joint_entropy,
+    symbolic_entropy,
+    symbolize,
 )
 from causationentropy.core.information.mutual_information import (
     gaussian_mutual_information,
     geometric_knn_mutual_information,
     kde_mutual_information,
     knn_mutual_information,
+    symbolic_mutual_information,
 )
 
 
@@ -371,6 +376,84 @@ def poisson_conditional_mutual_information(X, Y, Z):
         return cmi
 
 
+def symbolic_conditional_mutual_information(
+    X: np.ndarray,
+    Y: np.ndarray,
+    Z: Optional[np.ndarray] = None,
+    n_symbols: int = 2,
+    partition: str = "quantile",
+) -> float:
+    r"""
+    Plug-in conditional mutual information of symbolized data, in nats.
+
+    Every variable is symbolized with
+    :func:`~causationentropy.core.information.entropy.symbolize` (a partition of each
+    variable's range, as in Cafaro et al. [1]_), and the conditional mutual information
+    of Eq. (15) in [1]_ is estimated from empirical symbol frequencies:
+
+    .. math::
+
+        \hat I(X; Y | Z) = \hat H(X, Z) + \hat H(Y, Z) - \hat H(Z) - \hat H(X, Y, Z).
+
+    Parameters
+    ----------
+    X : np.ndarray, shape (n_samples,) or (n_samples, n_features_x)
+        First variable.
+    Y : np.ndarray, shape (n_samples,) or (n_samples, n_features_y)
+        Second variable.
+    Z : np.ndarray, shape (n_samples,) or (n_samples, n_features_z), or None
+        Conditioning variable. If None, returns the mutual information I(X; Y).
+    n_symbols : int, default=2
+        Number of symbols per variable.
+    partition : {'quantile', 'uniform', 'precomputed'}, default='quantile'
+        Partition type.
+
+    Returns
+    -------
+    float
+        Conditional mutual information estimate in nats.
+
+    Notes
+    -----
+    Cafaro et al. [1]_ show that the inferred causal structure can depend sensitively
+    on the partition, so ``n_symbols`` and ``partition`` should be chosen with care.
+    The plug-in estimate is biased upward when the number of occupied joint symbols is
+    not small compared with the sample size [2]_; in :func:`discover_network` the
+    permutation test uses the same estimator for its null distribution.
+    Time complexity is :math:`O(n \log n)` per call.
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> rng = np.random.default_rng(0)
+    >>> z = rng.integers(0, 2, 1000)
+    >>> x = z.copy(); y = z.copy()
+    >>> round(symbolic_conditional_mutual_information(x, y, z, partition="precomputed"), 6)
+    0.0
+
+    References
+    ----------
+    .. [1] C. Cafaro, W. M. Lord, J. Sun, and E. M. Bollt, "Causation entropy from
+           symbolic representations of dynamical systems," Chaos 25, 043106 (2015).
+           arXiv:1507.07262.
+    .. [2] L. Paninski, "Estimation of entropy and mutual information," Neural
+           Computation 15(6), 1191-1253 (2003).
+    """
+    if Z is None:
+        return symbolic_mutual_information(
+            X, Y, n_symbols=n_symbols, partition=partition
+        )
+
+    SX = symbolize(X, n_symbols=n_symbols, partition=partition)
+    SY = symbolize(Y, n_symbols=n_symbols, partition=partition)
+    SZ = symbolize(Z, n_symbols=n_symbols, partition=partition)
+
+    def H(*parts):
+        return symbolic_entropy(np.hstack(parts), partition="precomputed")
+
+    return H(SX, SZ) + H(SY, SZ) - H(SZ) - H(SX, SY, SZ)
+
+
 def conditional_mutual_information(
     X,
     Y,
@@ -380,6 +463,8 @@ def conditional_mutual_information(
     k=6,
     bandwidth="silverman",
     kernel="gaussian",
+    n_symbols: int = 2,
+    partition: str = "quantile",
 ):
     """
     Compute conditional mutual information using specified estimation method.
@@ -417,6 +502,7 @@ def conditional_mutual_information(
         - 'knn': k-nearest neighbor (KSG) estimator
         - 'geometric_knn': Geometric k-NN with manifold corrections
         - 'poisson': For discrete count data with Poisson assumptions
+        - 'symbolic': Plug-in estimate on symbolized data (Cafaro et al., 2015)
 
     metric : str or None, default=None
         Distance metric for k-NN based methods.
@@ -426,6 +512,10 @@ def conditional_mutual_information(
         Bandwidth parameter for KDE methods.
     kernel : str, default='gaussian'
         Kernel function for KDE methods.
+    n_symbols : int, default=2
+        Number of symbols per variable for the 'symbolic' method.
+    partition : {'quantile', 'uniform', 'precomputed'}, default='quantile'
+        Partition type for the 'symbolic' method.
 
     Returns
     -------
@@ -495,6 +585,11 @@ def conditional_mutual_information(
     elif method == "poisson":
         cmi = poisson_conditional_mutual_information(X, Y, Z)
 
+    elif method == "symbolic":
+        cmi = symbolic_conditional_mutual_information(
+            X, Y, Z, n_symbols=n_symbols, partition=partition
+        )
+
     else:
         supported_methods = [
             "gaussian",
@@ -503,6 +598,7 @@ def conditional_mutual_information(
             "knn",
             "geometric_knn",
             "poisson",
+            "symbolic",
         ]
         raise ValueError(
             f"Method '{method}' unavailable. Supported methods: {supported_methods}"

causationentropy/core/information/entropy.py

@@ -344,3 +344,169 @@ def poisson_joint_entropy(Cov):
     Ent1 = np.sum(poisson_entropy(U))
     Ent2 = np.sum(T)
     return Ent1 + Ent2
+
+
+def symbolize(
+    X: np.ndarray,
+    n_symbols: int = 2,
+    partition: str = "quantile",
+) -> np.ndarray:
+    r"""
+    Convert continuous data into symbols by partitioning each variable's range.
+
+    Following Cafaro et al. [1]_, a partition :math:`\mathcal{A} = \{A_0, \dots, A_m\}`
+    of the state space assigns to each observation :math:`x_t` the symbol
+
+    .. math::
+
+        s_t = \sum_{i=0}^{m} \chi_{A_i}(x_t) \cdot i,
+
+    where :math:`\chi_{A}` is the indicator function of the cell :math:`A`. Here the
+    partition is applied to each column (variable) separately, with ``n_symbols`` cells.
+
+    Parameters
+    ----------
+    X : np.ndarray, shape (n_samples,) or (n_samples, n_features)
+        Input data. A 1-D array is treated as a single column.
+    n_symbols : int, default=2
+        Number of cells (symbols) per variable. Must be at least 2. Ignored when
+        ``partition='precomputed'``.
+    partition : {'quantile', 'uniform', 'precomputed'}, default='quantile'
+        How the cells are chosen for each column:
+
+        - ``'quantile'``: approximately equal-frequency cells. Edges are placed at
+          midpoints between neighbouring distinct values, as close as possible to the
+          empirical quantiles, so tied values always share a cell (with
+          ``n_symbols=2`` this is a split near the median). Data with few distinct
+          values may therefore use fewer than ``n_symbols`` symbols.
+        - ``'uniform'``: equal-width cells between the column minimum and maximum.
+        - ``'precomputed'``: the data are already symbols (integer labels) and are
+          returned unchanged as integers.
+
+        For ``'uniform'``, a value equal to a cell edge is assigned to the upper cell.
+
+    Returns
+    -------
+    np.ndarray, shape (n_samples, n_features)
+        Integer symbols in ``{0, ..., n_symbols - 1}`` (or the given labels for
+        ``'precomputed'``).
+
+    Raises
+    ------
+    ValueError
+        If ``n_symbols < 2``, the partition type is unknown, the data contain NaN or
+        infinite values, or precomputed data are not integer-valued.
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> symbolize(np.array([0.1, 0.9, 0.4, 0.6]), n_symbols=2).ravel()
+    array([0, 1, 0, 1])
+
+    References
+    ----------
+    .. [1] C. Cafaro, W. M. Lord, J. Sun, and E. M. Bollt, "Causation entropy from
+           symbolic representations of dynamical systems," Chaos 25, 043106 (2015).
+           arXiv:1507.07262.
+    """
+    X = np.asarray(X)
+    if X.ndim == 1:
+        X = X.reshape(-1, 1)
+    X = X.astype(float)
+    if not np.all(np.isfinite(X)):
+        raise ValueError("symbolize requires finite values (no NaN or inf).")
+
+    if partition == "precomputed":
+        if not np.all(np.equal(np.mod(X, 1), 0)):
+            raise ValueError("partition='precomputed' requires integer-valued symbols.")
+        return X.astype(np.int64)
+
+    if n_symbols < 2:
+        raise ValueError(f"n_symbols must be at least 2, got {n_symbols}.")
+
+    symbols = np.empty(X.shape, dtype=np.int64)
+    for j in range(X.shape[1]):
+        col = X[:, j]
+        if partition == "quantile":
+            inner_edges = _quantile_edges(col, n_symbols)
+        elif partition == "uniform":
+            inner_edges = np.linspace(col.min(), col.max(), n_symbols + 1)[1:-1]
+        else:
+            raise ValueError(
+                f"Unknown partition '{partition}'. "
+                "Use 'quantile', 'uniform' or 'precomputed'."
+            )
+        symbols[:, j] = np.searchsorted(inner_edges, col, side="right")
+    return symbols
+
+
+def _quantile_edges(col: np.ndarray, n_symbols: int) -> np.ndarray:
+    """Equal-frequency cell edges placed between distinct values of ``col``."""
+    values, counts = np.unique(col, return_counts=True)
+    if values.size < 2:
+        return np.empty(0)
+    midpoints = (values[:-1] + values[1:]) / 2.0
+    share_below = np.cumsum(counts)[:-1] / col.size
+    targets = np.arange(1, n_symbols) / n_symbols
+    nearest = np.abs(share_below[None, :] - targets[:, None]).argmin(axis=1)
+    return np.unique(midpoints[nearest])
+
+
+def _joint_labels(S: np.ndarray) -> np.ndarray:
+    """Map each row of a 2-D symbol array to one integer label for its joint symbol."""
+    _, labels = np.unique(S, axis=0, return_inverse=True)
+    return labels.ravel()
+
+
+def symbolic_entropy(
+    X: np.ndarray,
+    n_symbols: int = 2,
+    partition: str = "quantile",
+) -> float:
+    r"""
+    Plug-in (maximum-likelihood) Shannon entropy of symbolized data, in nats.
+
+    The data are symbolized with :func:`symbolize`, and the joint symbol of each row is
+    counted. With empirical frequencies :math:`\hat P(s)`,
+
+    .. math::
+
+        \hat H(X) = -\sum_{s} \hat P(s) \log \hat P(s).
+
+    Parameters
+    ----------
+    X : np.ndarray, shape (n_samples,) or (n_samples, n_features)
+        Input data; several columns are treated as one joint variable.
+    n_symbols : int, default=2
+        Number of symbols per variable (see :func:`symbolize`).
+    partition : {'quantile', 'uniform', 'precomputed'}, default='quantile'
+        Partition type (see :func:`symbolize`).
+
+    Returns
+    -------
+    float
+        Entropy estimate in nats.
+
+    Notes
+    -----
+    The plug-in estimator is biased downward for small samples relative to the number
+    of occupied joint symbols [2]_. Time complexity is :math:`O(n \log n)`.
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> x = np.array([0, 1, 2, 3])
+    >>> round(symbolic_entropy(x, partition="precomputed"), 6)  # log(4)
+    1.386294
+
+    References
+    ----------
+    .. [1] C. Cafaro, W. M. Lord, J. Sun, and E. M. Bollt, "Causation entropy from
+           symbolic representations of dynamical systems," Chaos 25, 043106 (2015).
+    .. [2] L. Paninski, "Estimation of entropy and mutual information," Neural
+           Computation 15(6), 1191-1253 (2003).
+    """
+    labels = _joint_labels(symbolize(X, n_symbols=n_symbols, partition=partition))
+    counts = np.bincount(labels)
+    p = counts[counts > 0] / labels.size
+    return float(-np.sum(p * np.log(p))) + 0.0  # + 0.0 turns -0.0 into 0.0