draft r interface for dmi

[SpatLyu]

This PR introduces the R interface for Delayed Mutual Information (DMI) analysis, providing a complete pipeline from R-level function calls to the underlying C++ implementation powered by the KSG (Kraskov–Stögbauer–Grassberger) estimator. The implementation supports parallel computation across multiple lag steps.

---

Background & Theory

Mutual Information

Mutual Information (MI) quantifies the statistical dependence between two random variables $X$ and $Y$:

$$ I(X; Y) = H(X) + H(Y) - H(X, Y) $$

where $H(\cdot)$ denotes the Shannon entropy. Unlike Pearson correlation, MI captures both linear and nonlinear dependencies, making it a more general measure of association.

KSG Estimator

Estimating MI for continuous variables directly from the definition is intractable because the probability density function is typically unknown. The KSG estimator (Kraskov et al., 2004) addresses this by leveraging $k$-nearest neighbor distances in the joint and marginal spaces:

$$ \hat{I}(X; Y) = \psi(k) - \frac{1}{N} \sum_{i=1}^{N} \left[ \psi(n_x^{(i)}) + \psi(n_y^{(i)}) \right] + \psi(N) $$

where:

• $k$ is the number of nearest neighbors
• $\psi(\cdot)$ is the digamma function
• $n_x^{(i)}$ and $n_y^{(i)}$ are the number of points within the Chebyshev distance to the $k$-th neighbor of point $i$ in the marginal spaces $X$ and $Y$, respectively
• $N$ is the total number of samples

This estimator is asymptotically unbiased and does not require explicit density estimation or binning.

Delayed Mutual Information

Delayed Mutual Information extends MI to a time series context. Given a time series ${x_t}$, DMI measures the information shared between the present state $x_t$ and its past state $x_{t-\tau}$ at a given lag $\tau$:

$$ \text{DMI}(\tau) = I(x_t; x_{t-\tau}) $$

Computing DMI across a range of lags $\tau \in {\tau_1, \tau_2, \dots, \tau_m}$ reveals the temporal dependency structure of the series, which is widely used in:

• Determining optimal embedding delays for state-space reconstruction
• Identifying characteristic timescales in dynamical systems
• Nonlinear time series analysis

---

Algorithm

Step 1: Construct the Lag Matrix

Given:

• An input time series $\mathbf{v} = [v_0, v_1, \dots, v_{N-1}]$
• A set of prediction indices $\text{pred} = [p_1, p_2, \dots, p_n]$ (0-based)
• A vector of lag steps $\boldsymbol{\tau} = [\tau_1, \tau_2, \dots, \tau_m]$

A matrix $\mathbf{M}$ of size $(m + 1) \times n$ is constructed:

Row	Content	Description
0	$v_{p_j}$	Reference (current) values
$i$	$v_{p_j - \tau_{i-1}}$	Values lagged by $\tau_{i-1}$

Entries where $p_j - \tau_{i-1} < 0$ are filled with NaN and excluded from estimation.

Step 2: Input Validation & Preprocessing (R → C++)

• Convert R 1-based indices to C++ 0-based indices
• Remove indices where $p_j < \max(\boldsymbol{\tau})$ to ensure all lagged values are valid
• Sort and deduplicate the prediction indices
• Filter out indices where the target value is NaN

Step 3: Pairwise MI Estimation

For each lag $\tau_i$ ($i = 1, \dots, m$):

1. Extract row 0 (reference) and row $i$ (lagged) from $\mathbf{M}$
2. Remove columns where either value is NaN
3. Compute $\hat{I}(x_t; x_{t-\tau_i})$ using the KSG estimator

When threads > 1, MI computations across different lags are executed in parallel using RcppThread::parallelFor. Each lag's computation is independent, making this embarrassingly parallel.

---

R Interface

Function Signature

# now only works for time series data
dmi(data, target, tau = 1:10, pred = NULL, k = 3,
    base = exp(1), normalize = FALSE, threads = length(tau), ...)

Return Value

A named numeric vector of length length(tau), where each element is the estimated MI at the corresponding lag. Names follow the format "tau:<value>".

---

References

• Fraser, A.M., Swinney, H.L., 1986. Independent coordinates for strange attractors from mutual information. Physical Review A 33, 1134–1140.
• Kraskov, A., Stögbauer, H., Grassberger, P., 2004. Estimating mutual information. Physical Review E 69, 066138.