Inference

dit.inference estimates distributions and information quantities from samples. Plug-in counts live alongside bias-corrected entropy estimators, Markov-order tests and selection, surrogate and bootstrap generators, and k-nearest-neighbor / Kraskov–Stögbauer–Grassberger estimators for differential entropy and total correlation.

The kNN / KSG implementations optionally use scikit-learn (installed with the dit[optional] extra) when it is available.

Every randomized function takes a prng argument: None (draw from dit.math.prng), an integer seed, a numpy.random.Generator, or a numpy.random.RandomState.

From samples

distribution_from_data() builds a joint over words of length L from a sequence of symbols. dist_from_timeseries() treats each column of a multivariate series as a variable and appends a history_length past together with the present. Passing history_length='auto' chooses the history length with select_markov_order() (see below).

In [1]: from dit.inference import distribution_from_data, entropy_0, entropy_1

In [2]: data = [0, 0, 0, 1, 1, 1]

In [3]: d = distribution_from_data(data, L=1, base='linear')

In [4]: d.outcomes
Out[4]: (0, 1)

In [5]: entropy_0(data)
Out[5]: 1.0

The distributions returned here are plug-in (maximum-likelihood) estimates, so any measure computed from them inherits the plug-in bias. When the number of distinct words is not small compared to the number of samples, prefer the count-level estimators below. A Distribution indexes outcomes densely by each variable’s alphabet, so for large alphabets and word lengths the count-level estimators are also the only practical option.

Independent trials

Wrap separate recordings in Trials to pool them. Counts are summed over trials, but no word ever spans the boundary between two trials, which concatenating them would do. Every sample-based function in this module accepts Trials: block entropies, Markov-order tests (each trial gets its own Whittle surrogate set), surrogates, and the stationary bootstrap, which resamples whole trials.

from dit.inference import Trials, select_markov_order

select_markov_order(Trials([run_1, run_2, run_3]), max_order=4)

Undersampling

Estimators raise an UndersamplingWarning when there are more than a fifth as many distinct words as windows, or when the word space len(alphabet) ** L exceeds the number of windows. Markov-order test results report their effective sample size (MarkovOrderTest.n_windows).

Entropy estimators

entropy_from_counts() estimates an entropy from outcome counts, and block_entropy() and conditional_entropy_rate() apply it to the overlapping words of a sequence. The estimator argument selects:

'plugin' — the maximum-likelihood estimate (as entropy_0()).
'miller_madow' — plug-in plus \((K - 1) / 2N\) nats for \(K\) observed outcomes [Mil55].
'digamma' — \(\psi(N) - \sum_i \frac{n_i}{N} \psi(n_i)\) [Gra88] (as entropy_1()).
'grassberger' — Grassberger’s (2003) estimator [Gra03]. entropy_2() is the same estimator with \(\psi(N)\) in place of \(\ln N\).
'chao_shen' — the coverage-adjusted estimator, which accounts for unseen outcomes [CS03].
'nsb' — the Nemenman–Shafee–Bialek estimator, a mixture of Dirichlet priors chosen to make the prior on the entropy nearly uniform [NSB02]. It needs the number of possible outcomes (alphabet_size); block_entropy() supplies len(alphabet) ** L.

conditional_entropy_rate() estimates \(h_L = \operatorname{H}[X_L \mid X_{0:L}]\), which decreases to the entropy rate and equals it once \(L\) reaches the Markov order [CT06].

Two further entropy-rate estimates check it:

lz_entropy_rate() uses match lengths [KASW98]. It needs no history length, but converges slowly and usually from below.
entropy_rate_posterior() draws transition matrices of an order-\(k\) chain from their Dirichlet posterior and reports the entropy rate of each draw [SCHubler07], giving credible intervals.

Ordinal patterns

ordinal_patterns() maps each delay window \((x_{t-(m-1)\tau}, \ldots, x_t)\) of a real-valued series to the permutation that sorts it [BP02], giving \(m!\) symbols that are invariant to monotone transforms of the series. relative_rank() encodes the present \(x_t\) as its rank among the \(m\) preceding values, so that a present symbol shares no values with the preceding pattern; this avoids the leakage that biases symbolic transfer entropy [Kug12]. Ties are ranked by order of occurrence, broken by negligible noise, or kept as their own patterns [BQMS12].

permutation_entropy() is the entropy of the pattern distribution, and weighted_permutation_entropy() weights each window by its variance, which restores some amplitude information [FCKPrincipe13] and serves as a model-free predictability measure [GJB14].

Consecutive windows with \(\tau = 1\) share \(m - 1\) values, so each pattern can be followed by only \(m\) of the \(m!\) patterns. Those forbidden transitions are a property of the encoding, not of the process, so history-based estimates on ordinal symbols should use non-overlapping patterns.

Markov order

Choosing a history length by hand conflates estimator bias with dynamics. markov_order_test() tests the null hypothesis that a sequence is \(n\)-th order Markov against order \(n + 1\).

With method='exact' (the default), the null distribution comes from whittle_surrogates(): sequences drawn uniformly from all those sharing the observed \((n + 1)\)-gram counts and initial word. Those counts are sufficient for an \(n\)-th order chain, so the test is exact at any sample size [PH14]. Whittle’s formula [Bil61, Whi55] (whittle_count()) gives the size of that set. Pethel & Hahs show the asymptotic chi-squared test [AG57] (method='asymptotic') is badly anti-conservative at short lengths and higher orders. Surrogates are sampled in linear time as random Eulerian trails of the de Bruijn multigraph [KMUW96].

select_markov_order() tests \(n = 0, 1, \ldots\) in turn and returns the first order that is not rejected ('exact' or 'chi2'), or minimizes AIC or BIC over orders fitted on common targets [Kat81, Ton75]. Processes with infinite Markov order, such as strictly sofic processes, have no true order; all methods return longer histories as the sample grows.

In [6]: import numpy as np

In [7]: from dit.inference import markov_order_test, select_markov_order

In [8]: rng = np.random.default_rng(0)

In [9]: x = np.cumsum(rng.random(500) < 0.2) % 2  # a sticky first-order chain

In [10]: markov_order_test(x, order=0, prng=0).pvalue < 0.05
Out[10]: True

In [11]: select_markov_order(x, max_order=4, prng=0)
Out[11]: 1

Conditional mutual information and resampling

conditional_mutual_information() estimates \(I[X : Y \mid Z]\) from paired samples with any of the entropy estimators above; lagged copies of a series give time-delayed and transfer-entropy-style quantities.

stationary_bootstrap() [PR94] resamples a series in blocks of geometrically distributed length (whole trials for Trials). Each block junction creates words that were never observed, so statistics of lagged windows are biased toward independence unless the mean block length is much longer than the window.

Transfer entropy, its surrogate tests and confidence intervals, false discovery rate control, and network inference live in the companion package infoflow, which builds on these estimators and null generators.

API

class Trials(iterable=(), /)[source]

Independent realizations (trials) of the same process.

Wrap a list of sequences in Trials to tell the estimators in dit.inference that they are separate recordings: counts are pooled across trials, but no word ever spans the boundary between two trials. Trials may have different lengths.

Examples

>>> from dit.inference import Trials, block_entropy
>>> block_entropy(Trials([[0, 1], [1, 0]]), 2)  # (1, 1) never straddles the trials
1.0

class UndersamplingWarning[source]: Warned when counts are too sparse for plug-in (and most bias-corrected) information estimates to be trusted.

distribution_from_data(d, L, trim=True, base=None)[source]

Returns a distribution over words of length L from d.

The returned distribution is the naive estimate of the distribution, which assigns probabilities equal to the number of times a particular word appeared in the data divided by the total number of times a word could have appeared in the data.

Roughly, it corresponds to the stationary distribution of a maximum likelihood estimate of the transition matrix of an (L-1)th order Markov chain.

Parameters:

d (list or Trials) – A list of symbols to be converted into a distribution, or independent Trials, whose word counts are pooled without any word spanning two trials.
L (integer) – The length of the words for the distribution.
trim (bool) – If true, then words with zero probability are trimmed from the distribution.
base (int or string) – The desired base of the returned distribution. If None, then the value of dit.ditParams[‘base’] is used.

Notes

With trim=True only the observed words are built, rather than all len(alphabet) ** L of them. A Distribution still indexes its outcomes densely by each variable’s alphabet, so for large alphabets and word lengths use the count-level estimators (block_entropy(), conditional_entropy_rate()), which never build a distribution.

binned(ts, bins=2, style='maxent')[source]

Discretize a real-valued list.

Parameters:

ts (ndarray) – The real-valued array to bin
bins (int) – The number of bins to map the data into.
style (str, {'maxent', 'uniform'}) – The method of discretizing the data. Defaults to ‘maxent’.

Returns:

symb – The discretized time-series.

Return type:

ndarray

Raises:

ValueError – Raised if style is not a recognized method.

ordinal_patterns(ts, order=3, delay=1, ties='first', prng=None)[source]

Map each delay window of a real-valued series to its ordinal pattern.

The window ending at time t is \((x_{t-(m-1)\tau}, \ldots, x_{t-\tau}, x_t)\) with \(m\) = order and \(\tau\) = delay; its pattern is the permutation that sorts it [BP02]. Patterns are invariant to monotone transforms of the series.

Parameters:

ts (array_like or Trials) – A real-valued series, or independent Trials.
order (int) – The window length \(m \geq 2\).
delay (int) – The spacing \(\tau \geq 1\) between window elements.
ties ({'first', 'noise', 'distinct'}) – How equal values are ranked. 'first' ranks the earlier occurrence lower (the usual convention); 'noise' adds negligible random noise; 'distinct' keeps ties as their own patterns [BQMS12].
prng (None, int, Generator, RandomState) – Source of randomness for ties='noise'.

Returns:

patterns – Integer pattern codes, one per window (len(ts) - (order - 1) * delay); entry i is the window ending at t = i + (order - 1) * delay. With 'first' or 'noise' codes are the lexicographic index of the permutation, in 0 .. order! - 1. With 'distinct' they are base-order codes of the dense rank vector.

Return type:

np.ndarray or Trials

relative_rank(ts, order=3, delay=1, ties='first', prng=None)[source]

The rank of each value among the order values preceding it.

For time t this is the number of values among \(x_{t-\tau}, \ldots, x_{t-m\tau}\) that are smaller than \(x_t\) (or, with ties='first', not larger), in 0 .. order. Paired with ordinal_patterns() of the preceding window, it encodes the present of a series without sharing values with its past, which removes the leakage that biases symbolic transfer entropy [Kug12].

Parameters:

ts (array_like or Trials) – A real-valued series, or independent trials.
order (int) – The number \(m\) of preceding values.
delay (int) – The spacing \(\tau\) between them.
ties ({'first', 'noise', 'distinct'}) – 'first' counts earlier equal values as smaller; 'distinct' counts only strictly smaller values; 'noise' adds negligible noise first.
prng (None, int, Generator, RandomState) – Source of randomness for ties='noise'.

Returns:

ranks – Entry i is the rank of \(x_t\) with t = i + order * delay.

Return type:

np.ndarray or Trials

permutation_entropy(ts, order=3, delay=1, normalize=False, estimator='plugin', ties='first', prng=None)[source]

The permutation entropy of a real-valued series [BP02].

This is the entropy, in bits, of the distribution of ordinal_patterns() of order \(m\) and delay \(\tau\). For piecewise-monotone maps, the permutation entropy per symbol converges to the Kolmogorov–Sinai entropy as \(m \to \infty\).

Parameters:

ts (array_like or Trials) – The series, or independent trials (patterns are pooled).
order (int) – The pattern order \(m\) and delay \(\tau\).
delay (int) – The pattern order \(m\) and delay \(\tau\).
normalize (bool) – Divide by \(\log_2 m!\), giving a value in \([0, 1]\).
estimator (str) – See entropy_from_counts(); 'nsb' uses \(m!\) outcomes.
ties ({'first', 'noise', 'distinct'}) – See ordinal_patterns().
prng (None, int, Generator, RandomState) – Source of randomness for ties='noise'.

Returns:

h

Return type:

float

weighted_permutation_entropy(ts, order=3, delay=1, normalize=False, ties='first', prng=None)[source]

The weighted permutation entropy of a real-valued series [FCKPrincipe13].

Each window’s pattern is weighted by the variance of the window, so patterns of large fluctuations count more than patterns of noise-level wiggles. This restores some of the amplitude information that permutation_entropy() discards. Garland, James & Bradley use it as a model-free measure of predictability [GJB14].

Parameters:

ts (array_like or Trials) – The series, or independent trials (weights are pooled).
order (int) – The pattern order \(m\) and delay \(\tau\).
delay (int) – The pattern order \(m\) and delay \(\tau\).
normalize (bool) – Divide by \(\log_2 m!\).
ties ({'first', 'noise', 'distinct'}) – See ordinal_patterns().
prng (None, int, Generator, RandomState) – Source of randomness for ties='noise'.

Returns:

h

Return type:

float

dist_from_timeseries(observations, history_length=1, base='linear', max_history=None, prng=None)[source]

Infer a distribution from time series observations. For each variable, infer a history_length past and a single observation present.

Parameters:

observations (list of tuples, ndarray, Trials) – A sequence of observations in time order, or independent Trials whose windows are pooled without any window spanning two trials.
history_length (int or str) – The history length to utilize. If a string, the history length is the Markov order of the joint series estimated by select_markov_order(): 'auto' (the exact surrogate test) or one of its methods 'exact', 'chi2', 'aic', 'bic'.
base (float, str) – The base to use for the distribution. Defaults to ‘linear’.
max_history (int, None) – The largest history length considered when history_length is a string. Defaults to the largest L with M**(L + 1) <= N / 5, for M joint symbols and N observations.
prng (None, int, Generator, RandomState) – Source of randomness for 'auto' / 'exact'.

Returns:

ts – A distribution with the first half of the indices as the pasts of the various time series, and the second half their present values.

Return type:

Distribution

entropy_0(data, length=1)[source]

Estimate the entropy of length length subsequences in data.

Parameters:

data (iterable) – An iterable of samples.
length (int) – The length to group samples into.

Returns:

h0 – An estimate of the entropy.

Return type:

float

Notes

This returns the naive estimate of the entropy.

entropy_1(data, length=1)[source]

Estimate the entropy of length length subsequences in data.

Parameters:

data (iterable) – An iterable of samples.
length (int) – The length to group samples into.

Returns:

h1 – An estimate of the entropy.

Return type:

float

Notes

This is the digamma correction \(\hat{H} = \psi(N) - \sum_i \frac{n_i}{N} \psi(n_i)\), the leading term of Grassberger’s estimator [Gra88] (not the Miller–Madow correction; see entropy_from_counts()).

If M is the alphabet size and N is the number of samples, then the bias of this estimator is:: B ~ M/N

entropy_2(data, length=1)[source]

Estimate the entropy of length length subsequences in data.

Parameters:

data (iterable) – An iterable of samples.
length (int) – The length to group samples into.

Returns:

h2 – An estimate of the entropy.

Return type:

float

Notes

This is Grassberger’s (2003) estimator [Gra03], \(\hat{H} = \psi(N) - \sum_i \frac{n_i}{N} G(n_i)\) with \(G(n) = \psi(n) + \tfrac{1}{2}(-1)^n [\psi(\tfrac{n+1}{2}) - \psi(\tfrac{n}{2})]\), using \(\psi(N)\) in place of Grassberger’s \(\ln N\). It is not the Nemenman–Shafee–Bialek (NSB) estimator.

If M is the alphabet size and N is the number of samples, then the bias of this estimator is:: B ~ (M+1)/(2N)

entropy_from_counts(counts, estimator='plugin', alphabet_size=None)[source]

Estimate an entropy, in bits, from a vector of outcome counts.

Parameters:

counts (array_like) – The number of times each outcome was observed. Zeros are ignored.
estimator (str) –
One of:
- 'plugin' — the maximum-likelihood (naive) estimate.
- 'miller_madow' — plug-in plus \((K - 1) / 2N\) nats for \(K\) observed outcomes [Mil55].
- 'digamma' — \(\psi(N) - \sum_i \frac{n_i}{N} \psi(n_i)\) [Gra88] (as entropy_1()).
- 'grassberger' — Grassberger’s (2003) estimator [Gra03].
- 'chao_shen' — coverage-adjusted Horvitz–Thompson estimator, which accounts for unseen outcomes [CS03].
- 'nsb' — the Nemenman–Shafee–Bialek Bayesian estimator, which averages Dirichlet priors so that the prior on the entropy is nearly uniform [NSB02]. Needs alphabet_size.
alphabet_size (int, None) – The number of possible outcomes \(K\), used by 'nsb'. Defaults to len(counts), so pass a count vector that includes the zero counts or set this explicitly.

Returns:

h – The estimated entropy in bits.

Return type:

float

Warns:

UndersamplingWarning – If more than a fifth as many outcomes were observed as samples.

block_entropy(data, length, estimator='plugin')[source]

Estimate \(H[X_{0:L}]\), the entropy of overlapping length-length words.

Parameters:

data (iterable or Trials) – The sequence, or independent Trials (words never span two trials). Rows of a 2D array are joint symbols.
length (int) – The word length \(L\).
estimator (str) – See entropy_from_counts().

Returns:

h – The estimated block entropy in bits.

Return type:

float

conditional_entropy_rate(data, L, estimator='plugin')[source]

Estimate \(h_L = H[X_L \mid X_{0:L}] = H[X_{0:L+1}] - H[X_{0:L}]\).

As L grows, \(h_L\) decreases to the entropy rate; for a Markov chain of order \(R\) it equals the entropy rate for all \(L \geq R\) [CT06]. Use select_markov_order() to choose L. Bias-corrected estimator choices matter most when the number of observed words is not small compared to the sample size.

Parameters:

data (iterable or Trials) – The sequence, or independent Trials (words never span two trials). Rows of a 2D array are joint symbols.
L (int) – The history length.
estimator (str) – See entropy_from_counts().

Returns:

h – The estimated conditional entropy in bits.

Return type:

float

lz_entropy_rate(data)[source]

Estimate the entropy rate in bits per symbol from match lengths.

For each position \(i\), \(\Lambda_i\) is the length of the shortest string starting at \(i\) that does not appear starting at any earlier position. The increasing-window estimator of Kontoyiannis et al. [KASW98] is

\[\hat{h} = \left[\frac{1}{n} \sum_{i} \frac{\Lambda_i}{\log_2 (i + 1)}\right]^{-1}.\]

It is consistent for stationary ergodic processes and needs no history length. That makes it a check on conditional_entropy_rate(), whose answer depends on the chosen \(L\). Like the Lempel–Ziv compression rate [CT06], it converges slowly: \(\Lambda_i\) exceeds \(\log_2 i / h\) by a roughly constant amount, so the estimate usually approaches the entropy rate from below (e.g. about 0.9 bits for a fair coin at \(n = 10^4\)). Treat it as a cross-check rather than a precise estimate.

Parameters:: data (iterable or Trials) – The sequence. Rows of a 2D array are joint symbols. For Trials, the per-trial sums are pooled; matches are sought only within each trial.
Returns:: h – The estimated entropy rate in bits per symbol.
Return type:: float

Notes

Positions whose match runs to the end of the data are excluded, since their \(\Lambda_i\) is censored. The cost is roughly \(O(n^2 \log n)\) character comparisons, done by str.find.

entropy_rate_posterior(data, order, prior=0.5, n_samples=1000, prng=None)[source]

Posterior samples of the entropy rate under an order-order Markov model.

Each history’s next-symbol distribution gets an independent symmetric Dirichlet(prior) prior, so its posterior is Dirichlet(counts + prior). For each posterior draw of the transition probabilities, the entropy rate \(h_\mu = \sum_u \pi(u) H[X_k \mid X_{0:k} = u]\) is computed from the draw’s own stationary distribution \(\pi\) [SCHubler07].

Parameters:

data (iterable or Trials) – The observed sequence, or independent Trials. Rows of a 2D array are joint symbols.
order (int) – The Markov order \(k\); choose it with select_markov_order().
prior (float) – The Dirichlet concentration per symbol. 0.5 is the Jeffreys prior, 1 the uniform prior.
n_samples (int) – The number of posterior draws.
prng (None, int, Generator, RandomState) – Source of randomness.

Returns:

posterior

Return type:

EntropyRatePosterior

Notes

Every history, observed or not, is included, and a positive prior makes the sampled chain irreducible, so each draw has a unique stationary distribution. Histories never observed contribute prior-dominated rows weighted by their (small) stationary probability. The posterior reflects uncertainty given the order; it does not account for choosing the order from the same data.

class EntropyRatePosterior(samples: ndarray, order: int)[source]

Posterior samples of the entropy rate, from entropy_rate_posterior().

samples

Entropy rates (bits per symbol) of chains drawn from the posterior.

Type:: np.ndarray

order

The Markov order of the model.

Type:: int

markov_order_test(data, order=1, statistic='entropy_rate', method='exact', n_surrogates=1000, prng=None)[source]

Test the null hypothesis that data is Markov of order order against the alternative that it is of order order + 1.

Parameters:

data (iterable or Trials) – The observed sequence, or independent Trials. Rows of a 2D array are joint symbols.
order (int) – The Markov order of the null hypothesis.
statistic ({'entropy_rate', 'chi2'}) – 'entropy_rate' uses the plug-in conditional entropy \(H[X_{n+1} \mid X_{0:n+1}]\) with n = order; small values are evidence against the null. 'chi2' uses Pearson’s statistic for the (order + 2)-gram counts; large values are evidence against the null.
method ({'exact', 'asymptotic'}) – 'exact' compares against whittle_surrogates(), which is valid at any sample size [PH14]. 'asymptotic' uses the chi-squared limit [AG57]; for 'entropy_rate' this is the G-test, since \(G = 2 N \ln 2 \, (h_n - h_{n+1})\).
n_surrogates (int) – The number of surrogates for the exact test.
prng (None, int, Generator, RandomState) – Source of randomness.

Returns:

result

Return type:

MarkovOrderTest

Notes

Pethel & Hahs show that the asymptotic test is badly anti-conservative for short sequences and higher orders (e.g. a size of 0.22 instead of 0.05 for third-order chains on four symbols with 400 samples), while the exact test keeps its nominal size [PH14]. The exact p-value is (1 + #{surrogates at least as extreme}) / (1 + n_surrogates). An UndersamplingWarning is raised when the (order + 2)-gram counts are sparse; the exact test stays valid, but its power is low.

For Trials, counts are pooled over trials and each surrogate draws every trial independently from its own surrogate set, which is the exact conditional null given all trials’ counts.

class MarkovOrderTest(order: int, statistic: str, value: float, pvalue: float, dof: int, null: ndarray | None = None, n_windows: int | None = None)[source]

The result of markov_order_test().

order

The Markov order of the null hypothesis.

Type:: int

statistic

The name of the test statistic.

Type:: str

value

The observed value of the statistic.

Type:: float

pvalue

The p-value of the observed statistic.

Type:: float

dof

Degrees of freedom of the asymptotic chi-squared distribution.

Type:: int

null

The statistic evaluated on each surrogate (exact tests only).

Type:: np.ndarray or None

n_windows

The number of (order + 2)-windows (the effective sample size).

Type:: int or None

select_markov_order(data, max_order, method='exact', alpha=0.05, n_surrogates=1000, prng=None)[source]

Estimate the Markov order of data.

Parameters:

data (iterable or Trials) – The observed sequence, or independent Trials. Rows of a 2D array are joint symbols.
max_order (int) – The largest order considered.
method ({'exact', 'chi2', 'aic', 'bic'}) – 'exact' and 'chi2' test order n against n + 1 for n = 0, 1, ... and return the first order that is not rejected at level alpha ('exact' uses markov_order_test() with Whittle surrogates [PH14]; 'chi2' its asymptotic G-test). 'aic' and 'bic' minimize the information criterion over orders fitted on a common set of targets [Kat81, Ton75].
alpha (float) – The significance level for the testing methods.
n_surrogates (int) – The number of surrogates per test for 'exact'.
prng (None, int, Generator, RandomState) – Source of randomness.

Returns:

order – The estimated order.

Return type:

int

Notes

An UndersamplingWarning is raised when max_order + 1-words are sparse. The selected order then reflects what the data can resolve, and may be smaller than the process’s true order.

A process with infinite Markov order (e.g. a strictly sofic process such as the even process) has no true order; every method returns larger orders as the sample grows.

whittle_count(data, order=1)[source]

The natural log of the number of sequences sharing the observed (order + 1)-gram counts, first order-word, and last order-word.

Parameters:

data (iterable or Trials) – The observed sequence. Rows of a 2D array are joint symbols. For Trials, the log of the product of the per-trial counts.
order (int) – The Markov order of the null hypothesis.

Returns:

log_count – The natural log of the size of the surrogate set.

Return type:

float

Notes

This is Whittle’s formula [Bil61, Whi55], as used by [PH14]:

\[N_{uv}(F) = \frac{\prod_i F_{i\cdot}!}{\prod_{ij} F_{ij}!} C_{vu}\]

where \(F\) counts transitions between observed order-words and \(C_{vu}\) is the \((v, u)\) cofactor of \(\delta_{ij} - F_{ij} / F_{i\cdot}\). For order=0 it reduces to the multinomial coefficient of the symbol counts.

whittle_surrogates(data, order=1, n=1, prng=None)[source]

Draw surrogates that exactly preserve the (order + 1)-gram counts.

Each surrogate is drawn uniformly from the set of sequences with the same first order symbols and the same count of every (order + 1)-gram as data (and hence the same last order-word). Under the null hypothesis that data is an order-th order Markov chain, the observed sequence is itself a uniform draw from this set, so statistics computed on the surrogates give an exact null distribution at any sample size [PH14].

Parameters:

data (iterable or Trials) – The observed sequence. Rows of a 2D array are joint symbols. Each of several Trials is drawn independently from its own surrogate set.
order (int) – The Markov order to preserve. order=0 is a random permutation.
n (int) – The number of surrogates.
prng (None, int, Generator, RandomState) – Source of randomness.

Returns:

surrogates – Shape (n,) + np.shape(data), or a list of n Trials.

Return type:

np.ndarray or list of Trials

Notes

[PH14] sample by weighting each step with Whittle’s formula (whittle_count()), which costs a determinant per symbol. This implementation samples the same uniform distribution in linear time by drawing a random Eulerian trail of the order-word de Bruijn multigraph from a random last-exit arborescence [KMUW96].

shift_surrogates(data, n=1, prng=None)[source]

Circularly shift data by uniformly random nonzero offsets.

Shifting a source series relative to a target preserves each series’ own dynamics while destroying their alignment, which gives a null for directed measures such as transfer entropy.

Parameters:

data (iterable or Trials) – The sequence to shift (along the first axis); each of several Trials is shifted independently.
n (int) – The number of surrogates.
prng (None, int, Generator, RandomState) – Source of randomness.

Returns:

surrogates – Shape (n,) + np.shape(data), or a list of n Trials.

Return type:

np.ndarray or list of Trials

block_surrogates(data, block_length, n=1, prng=None)[source]

Randomly permute contiguous blocks of data.

Dependencies shorter than block_length are (mostly) kept, while longer ones are destroyed.

Parameters:

data (iterable or Trials) – The sequence to shuffle (along the first axis); blocks never move between Trials.
block_length (int) – The length of each block; the final block may be shorter.
n (int) – The number of surrogates.
prng (None, int, Generator, RandomState) – Source of randomness.

Returns:

surrogates – Shape (n,) + np.shape(data), or a list of n Trials.

Return type:

np.ndarray or list of Trials

conditional_mutual_information(x, y, z=None, estimator='plugin')[source]

Estimate \(I[X : Y \mid Z]\) in bits from paired samples.

Parameters:

x (array_like) – Samples of each variable; rows of 2D arrays are joint outcomes.
y (array_like) – Samples of each variable; rows of 2D arrays are joint outcomes.
z (array_like, None) – Samples of the conditioning variable. If None, estimate \(I[X : Y]\).
estimator (str) – The entropy estimator applied to each term; see entropy_from_counts().

Returns:

cmi

Return type:

float

stationary_bootstrap(data, n=1, mean_block_length=None, prng=None)[source]

Resample a series with the stationary bootstrap [PR94].

Each resample concatenates blocks that start at uniformly random positions and have geometrically distributed lengths, wrapping around the end of the series. Rows of multivariate data are resampled together. For Trials, whole trials are resampled with replacement instead (a cluster bootstrap), so no resample joins two trials.

Parameters:

data (array_like or Trials) – The series (along the first axis), or independent trials.
n (int) – The number of resamples.
mean_block_length (float, None) – The mean block length; defaults to N ** (1 / 3).
prng (None, int, Generator, RandomState) – Source of randomness.

Returns:

resamples – Shape (n,) + np.shape(data), or a list of n Trials.

Return type:

np.ndarray or list of Trials

differential_entropy_knn(data, rvs=None, k=4, noise=1e-10, prng=None)[source]

Compute the differential entropy of data using a k-nearest neighbors density estimator.

Parameters:

data (np.ndarray) – The data.
rvs (list) – The columns of data to use as the random variable. If None, use all.
k (int) – The number of nearest neighbors to use.
noise (float) – The standard deviation of the normally-distributed noise to add to data.
prng (None, int, Generator, RandomState) – Source of randomness for the noise.

Returns:

h – The estimated entropy.

Return type:

float

Notes

The entropy is returned in units of bits.

total_correlation_ksg(data, rvs, crvs=None, k=4, noise=1e-10, prng=None)

Compute the total correlation from observations. The total correlation is computed between the columns specified in rvs, given the columns specified in crvs. This utilizes the KSG kNN density estimator, and works on discrete, continuous, and mixed data.

Parameters:

data (np.array) – Real valued time series data.
rvs (iterable of iterables) – The columns for which the total correlation is to be computed.
crvs (iterable) – The columns upon which the total correlation should be conditioned.
k (int) – The number of nearest neighbors to use in estimating the local kernel density.
noise (float) – The standard deviation of the normally-distributed noise to add to the data.
prng (None, int, Generator, RandomState) – Source of randomness for the noise.

Returns:

tc – The total correlation of rvs given crvs.

Return type:

float

Notes

The total correlation is computed in bits, not nats as most KSG estimators do.

This implementation uses scikit-learn.