Docstrings · ComplexityMeasures.jl

ComplexityMeasures.ComplexityMeasures — Module

ComplexityMeasures.jl

ComplexityMeasures.jl is a Julia-based software for calculating 1000s of various kinds of probabilities, entropies, and other so-called complexity measures from a single-variable input datasets. For relational measures across many input datasets see its extension CausalityTools.jl. If you are a user of other programming languages (Python, R, MATLAB, ...), you can still use ComplexityMeasures.jl due to Julia's interoperability. For example, for Python use juliacall.

A careful comparison with alternative widely used software shows that ComplexityMeasures.jl outclasses the alternatives in several objective aspects of comparison, such as computational performance, overall amount of measures, reliability, and extendability. See the associated publication for more details.

The key features that it provides can be summarized as:

A rigorous framework for extracting probabilities from data, based on the mathematical formulation of probability spaces.
Several (12+) outcome spaces, i.e., ways to discretize data into probabilities.
Several estimators for estimating probabilities given an outcome space, which correct theoretically known estimation biases.
Several definitions of information measures, such as various flavours of entropies (Shannon, Tsallis, Curado...), extropies, and other complexity measures, that are used in the context of nonlinear dynamics, nonlinear timeseries analysis, and complex systems.
Several discrete and continuous (differential) estimators for entropies, which correct theoretically known estimation biases.
An extendable interface and well thought out API accompanied by dedicated developer documentation. This makes it trivial to define new outcome spaces, or new estimators for probabilities, information measures, or complexity measures and integrate them with everything else in the software without boilerplate code.

ComplexityMeasures.jl can be used as a standalone package, or as part of other projects in the JuliaDynamics organization, such as DynamicalSystems.jl or CausalityTools.jl.

To install it, run import Pkg; Pkg.add("ComplexityMeasures").

All further information is provided in the documentation, which you can either find online or build locally by running the docs/make.jl file.

Previously, this package was called Entropies.jl.

ComplexityMeasures.our_abstract_types — Constant

const our_abstract_types

The types from out package which we want to pretty-print.

ComplexityMeasures.AbstractBinning — Type

AbstractBinning

Supertype encompassing RectangularBinning and FixedRectangularBinning.

ComplexityMeasures.AddConstant — Type

AddConstant <: ProbabilitiesEstimator
AddConstant(; c = 1.0)

A generic add-constant probabilities estimator for counting-based OutcomeSpaces, where several literature estimators can be obtained tuning c. Currently $c$ can only be a scalar.

c = 1.0 is the Laplace estimator, or the "add-one" estimator.

Description

Probabilities for the $k$-th outcome $\omega_{k}$ are estimated as

\[p(\omega_k) = \dfrac{(n_k + c)}{n + mc},\]

where $m$ is the cardinality of the outcome space, and $n$ is the number of (encoded) input data points, and $n_k$ is the number of times the outcome $\omega_{k}$ is observed in the (encoded) input data points.

If the AddConstant estimator used with probabilities_and_outcomes, then $m$ is set to the number of observed outcomes. If used with allprobabilities_and_outcomes, then $m$ is set to the number of possible outcomes.

Unobserved outcomes are assigned nonzero probability!

Looking at the formula above, if $n_k = 0$, then unobserved outcomes are assigned a non-zero probability of $\dfrac{c}{n + mc}$. This means that if the estimator is used with allprobabilities_and_outcomes, then all outcomes, even those that are not observed, are assigned non-zero probabilities. This might affect your results if using e.g. missing_outcomes.

ComplexityMeasures.AlizadehArghami — Type

AlizadehArghami <: DifferentialInfoEstimator
AlizadehArghami(definition = Shannon(); m::Int = 1)

The AlizadehArghami estimator computes the Shannon differential information of a timeseries using the method from Alizadeh2010, with logarithms to the base specified in definition.

The AlizadehArghami estimator belongs to a class of differential entropy estimators based on order statistics. It only works for timeseries input.

Description

Assume we have samples $\bar{X} = \{x_1, x_2, \ldots, x_N \}$ from a continuous random variable $X \in \mathbb{R}$ with support $\mathcal{X}$ and density function$f : \mathbb{R} \to \mathbb{R}$. AlizadehArghami estimates the Shannon differential entropy

\[H(X) = \int_{\mathcal{X}} f(x) \log f(x) dx = \mathbb{E}[-\log(f(X))].\]

However, instead of estimating the above integral directly, it makes use of the equivalent integral, where $F$ is the distribution function for $X$:

\[H(X) = \int_0^1 \log \left(\dfrac{d}{dp}F^{-1}(p) \right) dp.\]

This integral is approximated by first computing the order statistics of $\bar{X}$ (the input timeseries), i.e. $x_{(1)} \leq x_{(2)} \leq \cdots \leq x_{(n)}$. The AlizadehArghami Shannon differential entropy estimate is then the the Vasicek estimate $\hat{H}_{V}(\bar{X}, m, n)$, plus a correction factor

\[\hat{H}_{A}(\bar{X}, m, n) = \hat{H}_{V}(\bar{X}, m, n) + \dfrac{2}{n}\left(m \log(2) \right).\]

ComplexityMeasures.AmplitudeAwareOrdinalPatterns — Type

AmplitudeAwareOrdinalPatterns <: OutcomeSpace
AmplitudeAwareOrdinalPatterns{m}(τ = 1, A = 0.5, lt = ComplexityMeasures.isless_rand)

A variant of OrdinalPatterns that also incorporates amplitude information, based on the amplitude-aware permutation entropy Azami2016. The outcome space and arguments are the same as in OrdinalPatterns.

Description

Similarly to WeightedOrdinalPatterns, a weight $w_i$ is attached to each ordinal pattern extracted from each state (or delay) vector $\mathbf{x}_i = (x_1^i, x_2^i, \ldots, x_m^i)$ as

\[w_i = \dfrac{A}{m} \sum_{k=1}^m |x_k^i | + \dfrac{1-A}{d-1} \sum_{k=2}^d |x_{k}^i - x_{k-1}^i|,\]

with $0 \leq A \leq 1$. When $A=0$ , only internal differences between the elements of $\mathbf{x}_i$ are weighted. Only mean amplitude of the state vector elements are weighted when $A=1$. With, $0<A<1$, a combined weighting is used.

ComplexityMeasures.ApproximateEntropy — Type

ApproximateEntropy <: ComplexityEstimator
ApproximateEntropy([x]; r = 0.2std(x), kwargs...)

An estimator for the approximate entropy Pincus1991 complexity measure, used with complexity.

The keyword argument r is mandatory if an input timeseries x is not provided.

Keyword arguments

r::Real: The radius used when querying for nearest neighbors around points. Its value should be determined from the input data, for example as some proportion of the standard deviation of the data.
m::Int = 2: The embedding dimension.
τ::Int = 1: The embedding lag.
base::Real = MathConstants.e: The base to use for the logarithm. Pincus (1991) uses the natural logarithm.

Description

Approximate entropy (ApEn) is defined as

\[ApEn(m ,r) = \lim_{N \to \infty} \left[ \phi(x, m, r) - \phi(x, m + 1, r) \right].\]

Approximate entropy is estimated for a timeseries x, by first embedding x using embedding dimension m and embedding lag τ, then searching for similar vectors within tolerance radius r, using the estimator described below, with logarithms to the given base (natural logarithm is used in Pincus, 1991).

Specifically, for a finite-length timeseries x, an estimator for $ApEn(m ,r)$ is

\[ApEn(m, r, N) = \phi(x, m, r, N) - \phi(x, m + 1, r, N),\]

where N = length(x) and

\[\phi(x, k, r, N) = \dfrac{1}{N-(k-1)\tau} \sum_{i=1}^{N - (k-1)\tau} \log{\left( \sum_{j = 1}^{N-(k-1)\tau} \dfrac{\theta(d({\bf x}_i^m, {\bf x}_j^m) \leq r)}{N-(k-1)\tau} \right)}.\]

Here, $\theta(\cdot)$ returns 1 if the argument is true and 0 otherwise, $d({\bf x}_i, {\bf x}_j)$ returns the Chebyshev distance between vectors ${\bf x}_i$ and ${\bf x}_j$, and the k-dimensional embedding vectors are constructed from the input timeseries $x(t)$ as

\[{\bf x}_i^k = (x(i), x(i+τ), x(i+2τ), \ldots, x(i+(k-1)\tau)).\]

Flexible embedding lag

In the original paper, they fix τ = 1. In our implementation, the normalization constant is modified to account for embeddings with τ != 1.

ComplexityMeasures.BayesianRegularization — Type

BayesianRegularization <: ProbabilitiesEstimator
BayesianRegularization(; a = 1.0)

The BayesianRegularization estimator is used with probabilities and related functions to estimate probabilities an m-element counting-based OutcomeSpace using Bayesian regularization of cell counts Hausser2009. See ProbabilitiesEstimator for usage.

Outcome space requirements

This estimator only works with counting-compatible outcome spaces.

Description

The BayesianRegularization estimator estimates the probability of the $k$-th outcome $\omega_{k}$ is

\[\omega_{k}^{\text{BayesianRegularization}} = \dfrac{n_k + a_k}{n + A},\]

where $n$ is the number of samples in the input data, $n_k$ is the observed counts for the outcome $\omega_{k}$, and $A = \sum_{i=1}^k a_k$.

Picking a

There are many common choices of priors, some of which are listed in Hausser2009. They include

a == 0, which is equivalent to the RelativeAmount estimator.
a == 0.5 (Jeffrey's prior)
a == 1 (Bayes-Laplace uniform prior)

a can also be chosen as a vector of real numbers. Then, if used with allprobabilities_and_outcomes, it is required that length(a) == total_outcomes(o, x), where x is the input data and o is the OutcomeSpace. If used with probabilities, then length(a) must match the number of observed outcomes (you can check this using probabilities_and_outcomes). The choice of a can severely impact the estimation errors of the probabilities, and the errors depend both on the choice of a and on the sampling scenario Hausser2009.

Assumptions

The BayesianRegularization estimator assumes a fixed and known m. Thus, using it with probabilities_and_outcomes and allprobabilities_and_outcomes will yield different results, depending on whether all outcomes are observed in the input data or not. For probabilities_and_outcomes, m is the number of observed outcomes. For allprobabilities_and_outcomes, m = total_outcomes(o, x), where o is the OutcomeSpace and x is the input data.

Note

If used with allprobabilities_and_outcomes, then outcomes which have not been observed may be assigned non-zero probabilities. This might affect your results if using e.g. missing_outcomes.

Examples

using ComplexityMeasures
x = cumsum(randn(100))
ps_bayes = probabilities(BayesianRegularization(a = 0.5), OrdinalPatterns{3}(), x)

See also: RelativeAmount, Shrinkage.

ComplexityMeasures.BubbleEntropy — Type

BubbleEntropy <: ComplexityEstimator
BubbleEntropy(; m = 3, τ = 1, definition = Renyi(q = 2))

The BubbleEntropy complexity estimator Manis2017 is just a difference between two entropies, each computed with the BubbleSortSwaps outcome space, for embedding dimensions m + 1 and m, respectively.

Manis2017 use the Renyi entropy of order q = 2 as the information measure definition, but here you can use any InformationMeasure. Manis2017 formulates the "bubble entropy" as the normalized measure below, while here you can also compute the unnormalized measure.

Definition

For input data x, the "bubble entropy" is computed by first embedding the input data using embedding dimension m and embedding delay τ (call the embedded pts y), and then computing the difference between the two entropies:

\[BubbleEn_T(τ) = H_T(y, m + 1) - H_T(y, m)\]

where $H_T(y, m)$ and $H_T(y, m + 1)$ are entropies of type $T$ (e.g. Renyi) computed with the input data x embedded to dimension $m$ and $m+1$, respectively. Use complexity to compute this non-normalized version. Use complexity_normalized to compute the normalized difference of entropies:

\[BubbleEn_H(τ)^{norm} = \dfrac{H_T(x, m + 1) - H_T(x, m)}{max(H_T(x, m + 1)) - max(H_T(x, m))},\]

where the maximum of the entropies for dimensions m and m + 1 are computed using information_maximum.

Example

using ComplexityMeasures
x = rand(1000)
est = BubbleEntropy(m = 5, τ = 3)
complexity(est, x)

ComplexityMeasures.BubbleSortSwaps — Type

BubbleSortSwaps <: CountBasedOutcomeSpace
BubbleSortSwaps(; m = 3, τ = 1)

The BubbleSortSwaps outcome space is based on Manis2017's paper on "bubble entropy".

Description

BubbleSortSwaps does the following:

Embeds the input data using embedding dimension m and embedding lag τ
For each state vector in the embedding, counting how many swaps are necessary for the bubble sort algorithm to sort state vectors.

For counts_and_outcomes, we then define a distribution over the number of necessary swaps. This distribution can then be used to estimate probabilities using probabilities_and_outcomes, which again can be used to estimate any InformationMeasure. An example of how to compute the "Shannon bubble entropy" is given below.

Outcome space

The outcome_space for BubbleSortSwaps are the integers 0:N, where N = (m * (m - 1)) / 2 + 1 (the worst-case number of swaps). Hence, the number of total_outcomes is N + 1.

Implements

codify. Returns the number of swaps required for each embedded state vector.

Examples

With the BubbleSortSwaps outcome space, we can easily compute a "bubble entropy" inspired by Manis2017. Note: this is not actually a new entropy - it is just a new way of discretizing the input data. To reproduce the bubble entropy complexity measure from Manis2017, see BubbleEntropy.

Examples

using ComplexityMeasures
x = rand(100000)
o = BubbleSortSwaps(; m = 5) # 5-dimensional embedding vectors
information(Shannon(; base = 2), o, x)

# We can also compute any other "bubble quantity", for example the 
# "Tsallis bubble extropy", with arbitrary probabilities estimators:
information(TsallisExtropy(), BayesianRegularization(), o, x)

ComplexityMeasures.BubbleSortSwapsEncoding — Type

BubbleSortSwapsEncoding <: Encoding
BubbleSortSwapsEncoding{m}()

BubbleSortSwapsEncoding is used with encode to encode a length-m input vector x into an integer in the range ω ∈ 0:((m*(m-1)) ÷ 2), by counting the number of swaps required for the bubble sort algorithm to sort x in ascending order.

decode is not implemented for this encoding.

Example

using ComplexityMeasures
x = [1, 5, 3, 1, 2]
e = BubbleSortSwapsEncoding{5}() # constructor type argument must match length of vector 
encode(e, x)

ComplexityMeasures.ChaoShen — Type

ChaoShen <: DiscreteInfoEstimatorShannon
ChaoShen(definition::Shannon = Shannon())

The ChaoShen estimator is used with information to compute the discrete Shannon entropy according to Chao2003.

Description

This estimator is a modification of the HorvitzThompson estimator that multiplies each plugin probability estimate by an estimate of sample coverage. If $f_1$ is the number of singletons (outcomes that occur only once) in a sample of length $N$, then the sample coverage is $C = 1 - \dfrac{f_1}{N}$. The Chao-Shen estimator of Shannon entropy is then

\[H_S^{CS} = -\sum_{i=1}^M \left( \dfrac{C p_i \log(C p_i)}{1 - (1 - C p_i)^N} \right),\]

where $N$ is the sample size and $M$ is the number of outcomes. If $f_1 = N$, then $f_1$ is set to $f_1 = N - 1$ to ensure positive entropy Arora2022.

ComplexityMeasures.CombinationEncoding — Type

CombinationEncoding <: Encoding
CombinationEncoding(encodings)

A CombinationEncoding takes multiple Encodings and creates a combined encoding that can be used to encode inputs that are compatible with the given encodings.

Encoding/decoding

When used with encode, each Encoding in encodings returns integers in the set 1, 2, …, n_e, where n_e is the total number of outcomes for a particular encoding. For k different encodings, we can thus construct the cartesian coordinate (c₁, c₂, …, cₖ) (cᵢ ∈ 1, 2, …, n_i), which can uniquely be identified by an integer. We can thus identify each unique combined encoding with a single integer.

When used with decode, the integer symbol is converted to its corresponding cartesian coordinate, which is used to retrieve the decoded symbols for each of the encodings, and a tuple of the decoded symbols are returned.

The total number of outcomes is prod(total_outcomes(e) for e in encodings).

Examples

using ComplexityMeasures

# We want to encode the vector `x`.
x = [0.9, 0.2, 0.3]

# To do so, we will use a combination of first-difference encoding, amplitude encoding,
# and ordinal pattern encoding.

encodings = (
    RelativeFirstDifferenceEncoding(0, 1; n = 2),
    RelativeMeanEncoding(0, 1; n = 5),
    OrdinalPatternEncoding(3) # x is a three-element vector
    )
c = CombinationEncoding(encodings)

# Encode `x` as integer
ω = encode(c, x)

# Decode symbol (into a vector of decodings, one for each encodings `e ∈ encodings`).
# In this particular case, the first two element will be left-bin edges, and
# the last element will be the decoded ordinal pattern (indices that would sort `x`).
d = decode(c, ω)

ComplexityMeasures.ComplexityEstimator — Type

ComplexityEstimator

Supertype for estimators for various complexity measures that are not entropies in the strict mathematical sense.

See complexity for all available estimators.

ComplexityMeasures.CompositeDownsampling — Type

CompositeDownsampling <: MultiScaleAlgorithm
CompositeDownsampling(; f::Function = Statistics.mean, scales = 1:8)

Composite multi-scale algorithm for multiscale entropy analysis Wu2013, used with multiscale to compute, for example, composite multiscale entropy (CMSE).

Description

Given a scalar-valued input time series x, the composite multiscale algorithm, like RegularDownsampling, downsamples and coarse-grains x by splitting it into non-overlapping windows of length s, and then constructing downsampled time series by applying the function f to each of the resulting length-s windows.

However, Wu2013 realized that for each scale s, there are actually s different ways of selecting windows, depending on where indexing starts/ends. These s different downsampled time series D_t(s, f) at each scale s are constructed as follows:

\[\{ D_{k}(s) \} = \{ D_{t, k}(s) \}_{t = 1}^{L}, = \{ f \left( \bf x_{t, k} \right) \} = \left\{ {f\left( (x_i)_{i = (t - 1)s + k}^{ts + k - 1} \right)} \right\}_{t = 1}^{L},\]

where L = floor((N - s + 1) / s) and 1 ≤ k ≤ s, such that $D_{i, k}(s)$ is the i-th element of the k-th downsampled time series at scale s.

Finally, compute $\dfrac{1}{s} \sum_{k = 1}^s g(D_{k}(s))$, where g is some summary function, for example information or complexity.

Keyword Arguments

scales. The downsampling levels. If scales is set to an integer, then this integer is taken as maximum number of scales (i.e. levels of downsampling), and downsampling is done over levels 1:scales. Otherwise, downsampling is done over the provided scales (which may be a range, or some specific scales (e.g. scales = [1, 5, 6]). The maximum scale level is length(x) ÷ 2, but to avoid applying the method to time series that are extremely short, consider limiting the maximum scale (e.g. scales = length(x) ÷ 5).

Relation to RegularDownsampling

The downsampled time series $D_{t, 1}(s)$ constructed using the composite multiscale method is equivalent to the downsampled time series $D_{t}(s)$ constructed using the RegularDownsampling method, for which k == 1 is fixed, such that only a single time series is returned.

Outcome space	Principle	Input data	Counting-compatible
`UniqueElements`	Count of unique elements	`Any`	✔
`ValueBinning`	Binning (histogram)	`Vector`, `StateSpaceSet`	✔
`OrdinalPatterns`	Ordinal patterns	`Vector`, `StateSpaceSet`	✔
`SpatialOrdinalPatterns`	Ordinal patterns in space	`Array`	✔
`Dispersion`	Dispersion patterns	`Vector`	✔
`SpatialDispersion`	Dispersion patterns in space	`Array`	✔
`CosineSimilarityBinning`	Cosine similarity	`Vector`	✔
`BubbleSortSwaps`	Swap counts when sorting	`Vector`	✔
`SequentialPairDistances`	Sequential state vector distances	`Vector`, `StateSpaceSet`	✔
`TransferOperator`	Binning (transfer operator)	`Vector`, `StateSpaceSet`	✖
`NaiveKernel`	Kernel density estimation	`StateSpaceSet`	✖
`WeightedOrdinalPatterns`	Ordinal patterns	`Vector`, `StateSpaceSet`	✖
`AmplitudeAwareOrdinalPatterns`	Ordinal patterns	`Vector`, `StateSpaceSet`	✖
`WaveletOverlap`	Wavelet transform	`Vector`	✖
`PowerSpectrum`	Fourier transform	`Vector`	✖

See also: RectangularBinning.

ComplexityMeasures.Tsallis — Type

Tsallis <: InformationMeasure
Tsallis(q; k = 1.0, base = 2)
Tsallis(; q = 1.0, k = 1.0, base = 2)

The Tsallis generalized order-q entropy Tsallis1988, used with information to compute an entropy.

base only applies in the limiting case q == 1, in which the Tsallis entropy reduces to Shannon entropy.

Description

The Tsallis entropy is a generalization of the Boltzmann-Gibbs entropy, with k standing for the Boltzmann constant. It is defined as

\[S_q(p) = \frac{k}{q - 1}\left(1 - \sum_{i} p[i]^q\right)\]

The maximum value of the Tsallis entropy is $k(L^{1 - q} - 1)/(1 - q)$, with $L$ the total_outcomes.

ComplexityMeasures.TsallisExtropy — Type

TsallisExtropy <: InformationMeasure
TsallisExtropy(; base = 2)

The Tsallis extropy Xue2023.

Description

TsallisExtropy is used with information to compute

\[J_T(P) = k \dfrac{N - 1 - \sum_{i=1}^N ( 1 - p[i])^q}{q - 1}\]

for a probability distribution $P = \{p_1, p_2, \ldots, p_N\}$, with the $\log$ at the given base. Alternatively, TsallisExtropy can be used with information_normalized, which ensures that the computed extropy is on the interval $[0, 1]$ by normalizing to to the maximal Tsallis extropy, given by

\[J_T(P) = \dfrac{(N - 1)N^{q - 1} - (N - 1)^q}{(q - 1)N^{q - 1}}\]

ComplexityMeasures.UniqueElements — Type

UniqueElements()

An OutcomeSpace based on straight-forward counting of distinct elements in a univariate time series or multivariate dataset. This is the same as giving no estimator to probabilities.

Outcome space

The outcome space is the unique sorted values of the input. Hence, input x is needed for a well-defined outcome_space.

Implements

codify. Used for encoding inputs where ordering matters (e.g. time series).

ComplexityMeasures.UniqueElementsEncoding — Type

UniqueElementsEncoding <: Encoding
UniqueElementsEncoding(x)

UniqueElementsEncoding is a generic encoding that encodes each xᵢ ∈ unique(x) to one of the positive integers. The xᵢ are encoded according to the order of their first appearance in the input data.

The constructor requires the input data x, since the number of possible symbols is length(unique(x)).

Example

using ComplexityMeasures
x = ['a', 2, 5, 2, 5, 'a']
e = UniqueElementsEncoding(x)
encode.(Ref(e), x) == [1, 2, 3, 2, 3, 1] # true

ComplexityMeasures.ValueBinning — Type

ValueBinning(b::AbstractBinning) <: OutcomeSpace

An OutcomeSpace based on binning the values of the data as dictated by the binning scheme b and formally computing their histogram, i.e., the frequencies of points in the bins. An alias to this is VisitationFrequency. Available binnings are subtypes of AbstractBinning.

The ValueBinning estimator has a linearithmic time complexity (n log(n) for n = length(x)) and a linear space complexity (l for l = dimension(x)). This allows computation of probabilities (histograms) of high-dimensional datasets and with small box sizes ε without memory overflow and with maximum performance. For performance reasons, the probabilities returned never contain 0s and are arbitrarily ordered.

ValueBinning(ϵ::Union{Real,Vector})

A convenience method that accepts same input as RectangularBinning and initializes this binning directly.

Outcomes

The outcome space for ValueBinning is the unique bins constructed from b. Each bin is identified by its left (lowest-value) corner, because bins are always left-closed-right-open intervals [a, b). The bins are in data units, not integer (cartesian indices units), and are returned as SVectors, i.e., same type as input data.

For convenience, outcome_space returns the outcomes in the same array format as the underlying binning (e.g., Matrix for 2D input).

For FixedRectangularBinning the outcome_space is well-defined from the binning, but for RectangularBinning input x is needed as well.

Implements

codify. Used for encoding inputs where ordering matters (e.g. time series).

ComplexityMeasures.ValueHistogram — Type

ValueHistogram

An alias for ValueBinning.

ComplexityMeasures.Vasicek — Type

Vasicek <: DifferentialInfoEstimator
Vasicek(definition = Shannon(); m::Int = 1)

The Vasicek estimator computes the Shannon differential information of a timeseries using the method from Vasicek1976, with logarithms to the base specified in definition.

The Vasicek estimator belongs to a class of differential entropy estimators based on order statistics, of which Vasicek1976 was the first. It only works for timeseries input.

Description

Assume we have samples $\bar{X} = \{x_1, x_2, \ldots, x_N \}$ from a continuous random variable $X \in \mathbb{R}$ with support $\mathcal{X}$ and density function$f : \mathbb{R} \to \mathbb{R}$. Vasicek estimates the Shannon differential entropy

\[H(X) = \int_{\mathcal{X}} f(x) \log f(x) dx = \mathbb{E}[-\log(f(X))].\]

However, instead of estimating the above integral directly, it makes use of the equivalent integral, where $F$ is the distribution function for $X$,

\[H(X) = \int_0^1 \log \left(\dfrac{d}{dp}F^{-1}(p) \right) dp\]

This integral is approximated by first computing the order statistics of $\bar{X}$ (the input timeseries), i.e. $x_{(1)} \leq x_{(2)} \leq \cdots \leq x_{(n)}$. The Vasicek Shannon differential entropy estimate is then

\[\hat{H}_V(\bar{X}, m) = \dfrac{1}{n} \sum_{i = 1}^n \log \left[ \dfrac{n}{2m} (\bar{X}_{(i+m)} - \bar{X}_{(i-m)}) \right]\]

Usage

In practice, choice of m influences how fast the entropy converges to the true value. For small value of m, convergence is slow, so we recommend to scale m according to the time series length n and use m >= n/100 (this is just a heuristic based on the tests written for this package).

ComplexityMeasures.VisitationFrequency — Type

VisitationFrequency

An alias for ValueBinning.

ComplexityMeasures.WaveletOverlap — Type

WaveletOverlap([wavelet]) <: OutcomeSpace

An OutcomeSpace based on the maximal overlap discrete wavelet transform (MODWT).

When used with probabilities, the MODWT is applied to a signal, then probabilities are computed as the (normalized) energies at different wavelet scales. These probabilities are used to compute the wavelet entropy according to Rosso2001. Input timeseries x is needed for a well-defined outcome space.

By default the wavelet Wavelets.WT.Daubechies{12}() is used. Otherwise, you may choose a wavelet from the Wavelets package (it must subtype OrthoWaveletClass).

Outcome space

The outcome space for WaveletOverlap are the integers 1, 2, …, N enumerating the wavelet scales. To obtain a better understanding of what these mean, we prepared a notebook you can view online. As such, this estimator only works for timeseries input and input x is needed for a well-defined outcome_space.

ComplexityMeasures.WeightedOrdinalPatterns — Type

WeightedOrdinalPatterns <: OutcomeSpace
WeightedOrdinalPatterns{m}(τ = 1, lt::Function = ComplexityMeasures.isless_rand)

A variant of OrdinalPatterns that also incorporates amplitude information, based on the weighted permutation entropy Fadlallah2013. The outcome space and arguments are the same as in OrdinalPatterns.

Description

For each ordinal pattern extracted from each state (or delay) vector, a weight is attached to it which is the variance of the vector. Probabilities are then estimated by summing the weights corresponding to the same pattern, instead of just counting the occurrence of the same pattern.

An implementation note

Note: in equation 7, section III, of the original paper, the authors write

\[w_j = \dfrac{1}{m}\sum_{k=1}^m (x_{j-(k-1)\tau} - \mathbf{\hat{x}}_j^{m, \tau})^2.\]

*But given the formula they give for the arithmetic mean, this is not the variance of the delay vector $\mathbf{x}_i$, because the indices are mixed: $x_{j+(k-1)\tau}$ in the weights formula, vs. $x_{j+(k+1)\tau}$ in the arithmetic mean formula. Here, delay embedding and computation of the patterns and their weights are completely separated processes, ensuring that we compute the arithmetic mean correctly for each vector of the input dataset (which may be a delay-embedded timeseries).

ComplexityMeasures.Zhu — Type

Zhu <: DifferentialInfoEstimator
Zhu(; definition = Shannon(), k = 1, w = 0)

The Zhu estimator Zhu2015 is an extension to KozachenkoLeonenko, and computes the Shannon differential information of a multi-dimensional StateSpaceSet, with logarithms to the base specified in definition.

Description

\[H(X) = \int_{\mathcal{X}} f(x) \log f(x) dx = \mathbb{E}[-\log(f(X))]\]

by approximating densities within hyperrectangles surrounding each point xᵢ ∈ x using using k nearest neighbor searches. w is the Theiler window, which determines if temporal neighbors are excluded during neighbor searches (defaults to 0, meaning that only the point itself is excluded when searching for neighbours).

ComplexityMeasures.ZhuSingh — Type

ZhuSingh <: DifferentialInfoEstimator
ZhuSingh(definition = Shannon(); k = 1, w = 0)

The ZhuSingh estimator Zhu2015 computes the Shannon differential information of a multi-dimensional StateSpaceSet, with logarithms to the base specified in definition.

Description

Assume we have samples $\{\bf{x}_1, \bf{x}_2, \ldots, \bf{x}_N \}$ from a continuous random variable $X \in \mathbb{R}^d$ with support $\mathcal{X}$ and density function$f : \mathbb{R}^d \to \mathbb{R}$. ZhuSingh estimates the Shannon differential entropy

\[H(X) = \int_{\mathcal{X}} f(x) \log f(x) dx = \mathbb{E}[-\log(f(X))].\]

Like Zhu, this estimator approximates probabilities within hyperrectangles surrounding each point xᵢ ∈ x using using k nearest neighbor searches. However, it also considers the number of neighbors falling on the borders of these hyperrectangles. This estimator is an extension to the entropy estimator in Singh2003.

ComplexityMeasures.AAPE — Function

AAPE(x, A::Real = 0.5, m::Int = length(x))

Encode relative amplitude information of the elements of a.

A = 1 emphasizes only average values.
A = 0 emphasizes changes in amplitude values.
A = 0.5 equally emphasizes average values and changes in the amplitude values.

ComplexityMeasures.allcounts_and_outcomes — Method

allcounts_and_outcomes(o::OutcomeSpace, x::Array_or_SSSet) → (cts::Counts{<:Integer, 1}, Ω)

Like counts_and_outcomes, but ensures that all outcomes Ωᵢ ∈ Ω, where Ω = outcome_space(o, x)), are included.

Outcomes that do not occur in the data x get a 0 count.

ComplexityMeasures.allprobabilities_and_outcomes — Method

allprobabilities_and_outcomes(est::ProbabilitiesEstimator, x::Array_or_SSSet) → (p::Probabilities, outs)
allprobabilities_and_outcomes(o::OutcomeSpace, x::Array_or_SSSet) → (p::Probabilities, outs)

The same as probabilities_and_outcomes, but ensures that outcomes with 0 probability are explicitly added in the returned vector. This means that p[i] is the probability of ospace[i], with ospace =outcome_space(est, x).

This function is useful in cases where one wants to compare the probability mass functions of two different input data x, y under the same estimator. E.g., to compute the KL-divergence of the two PMFs assumes that the obey the same indexing. This is not true for probabilities even with the same est, due to the skipping of 0 entries, but it is true for allprobabilities_and_outcomes.

ComplexityMeasures.apply_multiscale — Function

apply_multiscale(alg::MultiScaleAlgorithm, f::Function, args...)

Define multiscale dispatch for the function f (either information, complexity or their normalized variants) to downsampled timeseries resulting from coarse-graining last(args) (the input data) using coarse-graining algorithm alg with arguments args[1:end-1] (the estimation parameters).

ComplexityMeasures.ball_volume — Method

Volume of a unit ball in R^d.

ComplexityMeasures.base10_to_factorial — Function

base10_to_factorial(s::Int,
    ndigits::Int = ndigits_in_factorial_base(s)) → f::SVector{ndigits, Int}

Convert a base-10 integer to its factorial number system representation. f is a vector where f[k] is the multiplier of factorial(k - 1).

For example, the base-10 integer 567, in the factorial number system, is $4\cdot 5! + 3\cdot 4! + 2\cdot 3! + 1\cdot 2! + 1\cdot 1! + 0\cdot 0!$. For this example, base10_to_factorial would return the SVector [4, 3, 2, 1, 1, 0].

ndigits fixes the number of digits in f (this just prepends a zero to f for each extraneous radix/base). This is useful when using factorial number for decoding Lehmer codes into permutations

ComplexityMeasures.cartesian_bin_index — Method

cartesian_bin_index(e::RectangularBinEncoding, point::SVector)

Return the cartesian index of the given point within the binning encapsulated in e. Internal function called by encode.

ComplexityMeasures.center_neighborhood! — Method

center_neighborhood!(C, c, xᵢ, neighbors)

Center the point xᵢ, as well as each of its neighboring points nⱼ ∈ neighbors, to the (precomputed) centroid c of the points {xᵢ, n₁, n₂, …, nₖ}, and store the centered vectors in the pre-allocated vector of vectors C.

ComplexityMeasures.codify — Function

codify(o::OutcomeSpace, x::Vector) → s::Vector{Int}
codify(o::OutcomeSpace, x::AbstractStateSpaceSet{D}) → s::NTuple{D, Vector{Int}

Codify x according to the outcome space o. If x is a Vector, then a Vector{<:Integer} is returned. If x is a StateSpaceSet{D}, then symbolization is done column-wise and an NTuple{D, Vector{<:Integer}} is returned, where D = dimension(x).

Description

The reason this function exists is that we don't always want to encode the entire input x at once. Sometimes, it is desirable to first apply some transformation to x first, then apply Encodings in a point-wise manner in the transformed space. (the OutcomeSpace dictates this transformation). This is useful for encoding timeseries data.

The length of the returned s depends on the OutcomeSpace. Some outcome spaces preserve the input data length (e.g. UniqueElements), while some outcome spaces (e.g. OrdinalPatterns) do e.g. delay embeddings before encoding, so that length(s) < length(x).

ComplexityMeasures.complexity — Method

complexity(c::ComplexityEstimator, x) → m::Real

Estimate a complexity measure according to c for input data x, where c is an instance of any subtype of ComplexityEstimator:

ApproximateEntropy.
LempelZiv76.
MissingDispersionPatterns.
ReverseDispersion.
SampleEntropy.
BubbleEntropy.
StatisticalComplexity.

ComplexityMeasures.complexity_normalized — Method

complexity_normalized(c::ComplexityEstimator, x) → m::Real ∈ [a, b]

The same as complexity, but the result is normalized to the interval [a, b], where [a, b] depends on c.

ComplexityMeasures.compute_ϕ — Method

compute_ϕ(x::AbstractVector{T}; r = 0.2 * Statistics.std(x), k::Int = 2,
    τ::Int = 1, base = MathConstants.e) where T <: Real

Construct the embedding

\[u = \{{\bf u}_n \}_{n = 1}^{N - k + 1} = \{[x(i), x(i + 1), \ldots, x(i + k - 1)]\}_{n = 1}^{N - k + 1}\]

and use a tree-and-nearest-neighbor search approach to compute

\[\phi^k(r) = \dfrac{1}{N - kτ + 1} \sum_{i}^{N - kτ + 1} \log_{b}{(C_i^k(r))},\]

taking logarithms to base $b$, and where

\[C_i^k(r) = \textrm{number of } j \textrm{ such that } d({\bf u}_i, {\bf u}_j) < r,\]

where $d$ is the maximum (Chebyshev) distance, r is the tolerance, and N is the length of the original scalar-valued time series x.

ComplexityMeasures.convert_logunit — Method

convert_logunit(h_a::Real, base_from, base_to) → h_b

Convert a number h_a computed with logarithms to base base_from to an entropy h_b computed with logarithms to base base_to. This can be used to convert the "unit" of an entropy.

ComplexityMeasures.counts — Function

counts(o::OutcomeSpace, x) → cts::Counts

Compute the same counts as in the counts_and_outcomes function, with two differences:

Do not explicitly return the outcomes.
If the outcomes are not estimated for free while estimating the counts, a special integer type is used to enumerate the outcomes, to avoid the computational cost of estimating the outcomes.

ComplexityMeasures.counts_and_outcomes — Method

counts_and_outcomes(o::OutcomeSpace, x) → (cts::Counts, Ω)

Discretize/encode x (which must be sortable) into a finite set of outcomes Ω specified by the provided OutcomeSpace o, and then count how often each outcome Ωᵢ ∈ Ω (i.e. each "discretized value", or "encoded symbol") appears.

Return a tuple where the first element is a Counts instance, which is vector-like and contains the counts, and where the second element Ω are the outcomes corresponding to the counts, such that cts[i] is the count for the outcome Ω[i].

The outcomes are actually included in cts, and you can use the outcomes function on the cts to get them. counts_and_outcomes returns both for backwards compatibility.

counts_and_outcomes(x) → cts::Counts

If no OutcomeSpace is specified, then UniqueElements is used as the outcome space.

Description

For OutcomeSpaces that uses encode to discretize, it is possible to count how often each outcome $\omega_i \in \Omega$, where $\Omega$ is the set of possible outcomes, is observed in the discretized/encoded input data. Thus, we can assign to each outcome $\omega_i$ a count $f(\omega_i)$, such that $\sum_{i=1}^N f(\omega_i) = N$, where $N$ is the number of observations in the (encoded) input data. counts returns the counts $f(\omega_i)_{obs}$ and outcomes only for the observed outcomes $\omega_i^{obs}$ (those outcomes that actually appear in the input data). If you need the counts for unobserved outcomes as well, use allcounts_and_outcomes.

ComplexityMeasures.decode — Function

decode(c::Encoding, i::Integer) -> ω

Decode an encoded element i into the outcome ω ∈ Ω it corresponds to. Ω is the outcome_space that uses encoding c.

ComplexityMeasures.distance_to_whitenoise — Method

distance_to_whitenoise(estimator::ReverseDispersion, p::Probabilities;
    normalize = false)

Compute the distance of the probability distribution p from a uniform distribution, given the parameters of estimator (which must be known beforehand).

If normalize == true, then normalize the value to the interval [0, 1] by using the parameters of estimator.

Used to compute reverse dispersion entropy(ReverseDispersion; Li et al., 2019Li2019).

ComplexityMeasures.downsample — Method

downsample(algorithm::MultiScaleAlgorithm, s::Int, x)

Downsample and coarse-grain x to scale s according to the given MultiScaleAlgorithm. The return type depends on algorithm.

ComplexityMeasures.encode — Function

encode(c::Encoding, χ) -> i::Int

Encode an element χ ∈ x of input data x (those given to e.g., counts) into the positive integers using encoding c. The special value of i = -1 is used as a return value for inappropriate elements χ that cannot be encoded according to c.

ComplexityMeasures.entropy — Method

entropy(args...)

entropy is nothing more than a call to information that will simply throw an error if used with an information measure that is not an entropy.

ComplexityMeasures.entropy_approx — Method

entropy_approx(x; m = 2, τ = 1, r = 0.2 * Statistics.std(x), base = MathConstants.e)

Convenience syntax for computing the approximate entropy (Pincus, 1991) for timeseries x.

This is just a wrapper for complexity(ApproximateEntropy(; m, τ, r, base), x) (see also ApproximateEntropy).

ComplexityMeasures.entropy_complexity — Method

entropy_complexity(c::StatisticalComplexity, x) → (h, compl)

Return a information measure h and the corresponding StatisticalComplexity value compl.

Useful when wanting to plot data on the "entropy-complexity plane". See also entropy_complexity_curves.

ComplexityMeasures.entropy_complexity_curves — Method

entropy_complexity_curves(c::StatisticalComplexity;
    num_max=1, num_min=1000) -> (min_entropy_complexity, max_entropy_complexity)

Calculate the maximum complexity-entropy curve for the statistical complexity according to Rosso2007 for num_max * total_outcomes(c.o) different values of the normalized information measure of choice (in case of the maximum complexity curves) and num_min different values of the normalized information measure of choice (in case of the minimum complexity curve).

This function can also be used to compute the maximum "complexity-extropy curve" if c.hest is e.g. ShannonExtropy, which is the equivalent of the complexity-entropy curves, but using extropy instead of entropy.

Description

The way the statistical complexity is designed, there is a minimum and maximum possible complexity for data with a given value of an information measure. The calculation time of the maximum complexity curve grows as O(total_outcomes(c.o)^2), and thus takes very long for high numbers of outcomes. This function is inspired by S. Sippels implementation in statcomp Sippel2016.

This function will work with any ProbabilitiesEstimator where total_outcomes is known a priori.

ComplexityMeasures.entropy_dispersion — Method

entropy_dispersion(x; base = 2, kwargs...)

Compute the dispersion entropy. This function is just a convenience call to:

est = Dispersion(kwargs...)
information(Shannon(base), est, x)

See Dispersion for more info.

ComplexityMeasures.entropy_distribution — Method

entropy_distribution(x; τ = 1, m = 3, n = 3, base = 2)

Compute the distribution entropy Li2015 of x using embedding dimension m with delay/lag τ, using the Chebyshev distance metric, and using an n-element equally-spaced binning over the distribution of distances to estimate probabilities.

This function is just a convenience call to:

x = rand(1000000)
o = SequentialPairDistances(x, n, m, τ, metric = Chebyshev())
h = information(Shannon(base = 2), o, x)

See SequentialPairDistances for more info.

ComplexityMeasures.entropy_permutation — Method

entropy_permutation(x; τ = 1, m = 3, base = 2)

Compute the permutation entropy of x of order m with delay/lag τ. This function is just a convenience call to:

est = OrdinalPatterns(; m, τ)
information(Shannon(base), est, x)

See OrdinalPatterns for more info. Similarly, one can use WeightedOrdinalPatterns or AmplitudeAwareOrdinalPatterns for the weighted/amplitude-aware versions.

ComplexityMeasures.entropy_sample — Method

entropy_sample(x; r = 0.2std(x), m = 2, τ = 1, normalize = true)

Convenience syntax for estimating the (normalized) sample entropy (Richman & Moorman, 2000) of timeseries x.

This is just a wrapper for complexity(SampleEntropy(; r, m, τ, base), x).

ComplexityMeasures.entropy_wavelet — Method

entropy_wavelet(x; wavelet = Wavelets.WT.Daubechies{12}(), base = 2)

Compute the wavelet entropy. This function is just a convenience call to:

est = WaveletOverlap(wavelet)
information(Shannon(base), est, x)

See WaveletOverlap for more info.

ComplexityMeasures.fasthist! — Method

fasthist!(x) → c::Vector{Int}

Count the occurrences c of the unique data values in x, so that c[i] is the number of times the value sort!(unique(x))[i] occurs. Hence, this method is useful mostly when x contains integer or categorical data.

Prior to counting, x is sorted, so this function also mutates x. Therefore, it is called with copy in higher level API when necessary. This function works for any x for which sort!(x) works.

ComplexityMeasures.fasthist! — Method

fasthist!(x, weights) → c::Vector{Real}

Similar to fasthist!(x), but here the weights are summed up for each unique entry of x. x is sorted just like in fasthist!(x).

ComplexityMeasures.fasthist — Method

fasthist(c::RectangularBinEncoding, x::Vector_or_SSSet)

Intermediate method that runs fasthist! in the encoded space and returns the encoded space histogram (counts) and corresponding bins. Also skips any instances of out-of-bound points for the histogram.

ComplexityMeasures.hidefields — Method

hidefields(::Type{T})

Returns an iterable of symbols incidating fields to hide for instances of type T.

ComplexityMeasures.information — Method

information(est::DifferentialInfoEstimator, x) → h::Real

Estimate a differential information measure using the provided DifferentialInfoEstimator and input data x.

Description

The overwhelming majority of differential estimators estimate the Shannon entropy. If the same estimator can estimate different information measures (e.g. it can estimate both Shannon and Tsallis), then the information measure is provided as an argument to the estimator itself.

See the table of differential information measure estimators in the docs for all differential information measure estimators.

Currently, unlike for the discrete information measures, this method doesn't involve explicitly first computing a probability density function and then passing this density to an information measure definition. But in the future, we want to establish a density API similar to the probabilities API.

Examples

To compute the differential version of a measure, give it as the first argument to a DifferentialInfoEstimator and pass it to information.

x = randn(1000)
h_sh = information(Kraskov(Shannon()), x)
h_vc = information(Vasicek(Shannon()), x)

A normal distribution has a base-e Shannon differential entropy of 0.5*log(2π) + 0.5 nats.

est = Kraskov(k = 5, base = ℯ) # Base `ℯ` for nats.
h = information(est, randn(2_000_000))
abs(h - 0.5*log(2π) - 0.5) # ≈ 0.0001

ComplexityMeasures.information — Method

information([die::DiscreteInfoEstimator,] [est::ProbabilitiesEstimator,] o::OutcomeSpace, x) → h::Real
information(o::OutcomeSpace, x) → h::Real

Estimate a discrete information measure from input data x using the provided DiscreteInfoEstimator and ProbabilitiesEstimator over the given OutcomeSpace.

As an alternative, you can provide an InformationMeasure for the first argument (die) which will default to PlugIn estimation) for the information estimation. You may also skip the first argument (die), in which case Shannon() will be used. You may also skip the second argument (est), which will default to the RelativeAmount probabilities estimator. Note that some information measure estimators (e.g., GeneralizedSchuermann) operate directly on counts and hence ignore est.

information([e::DiscreteInfoEstimator,] p::Probabilities) → h::Real
information([e::DiscreteInfoEstimator,] c::Counts) → h::Real

Like above, but estimate the information measure from the pre-computed Probabilities p or Counts. Counts are converted into probabilities using RelativeAmount, unless the estimator e uses counts directly.

See also: information_maximum, information_normalized for a normalized version.

Examples (naive estimation)

The simplest way to estimate a discrete measure is to provide the InformationMeasure directly in combination with an OutcomeSpace. This will use the "naive" PlugIn estimator for the measure, and the "naive" RelativeAmount estimator for the probabilities.

x = randn(100) # some input data
o = ValueBinning(RectangularBinning(5)) # a 5-bin histogram outcome space
h_s = information(Shannon(), o, x)

Here are some more examples:

x = [rand(Bool) for _ in 1:10000] # coin toss
ps = probabilities(x) # gives about [0.5, 0.5] by definition
h = information(ps) # gives 1, about 1 bit by definition (Shannon entropy by default)
h = information(Shannon(), ps) # syntactically equivalent to the above
h = information(Shannon(), UniqueElements(), x) # syntactically equivalent to above
h = information(Renyi(2.0), ps) # also gives 1, order `q` doesn't matter for coin toss
h = information(OrdinalPatterns(;m=3), x) # gives about 2, again by definition

Examples (bias-corrected estimation)

It is known that both PlugIn estimation for information measures and RelativeAmount estimation for probabilities are biased. The scientific literature abounds with estimators that correct for this bias, both on the measure-estimation level and on the probability-estimation level. We thus provide the option to use any DiscreteInfoEstimator in combination with any ProbabilitiesEstimator for improved estimates. Note that custom probabilites estimators will only work with counting-compatible OutcomeSpace.

x = randn(100)
o = ValueBinning(RectangularBinning(5))

# Estimate Shannon entropy estimation using various dedicated estimators
h_s = information(MillerMadow(Shannon()), RelativeAmount(), o, x)
h_s = information(HorvitzThompson(Shannon()), Shrinkage(), o, x)
h_s = information(Schuermann(Shannon()), Shrinkage(), o, x)

# Estimate information measures using the generic `Jackknife` estimator
h_r = information(Jackknife(Renyi()), Shrinkage(), o, x)
j_t = information(Jackknife(TsallisExtropy()), BayesianRegularization(), o, x)
j_r = information(Jackknife(RenyiExtropy()), RelativeAmount(), o, x)

ComplexityMeasures.information_maximum — Method

information_maximum(e::InformationMeasure, o::OutcomeSpace [, x])

Return the maximum value of the given information measure can have, given input data x and the given outcome space (the OutcomeSpace may also be specified by a ProbabilitiesEstimator).

Like in outcome_space, for some outcome spaces, the possible outcomes are known without knowledge of input x, in which case the function dispatches to information_maximum(e, o).

information_maximum(e::InformationMeasure, L::Int)

The same as above, but computed directly from the number of total outcomes L.

ComplexityMeasures.information_normalized — Method

information_normalized([e::DiscreteInfoEstimator,] [est::ProbabilitiesEstimator,] o::OutcomeSpace, x) → h::Real

Estimate the normalized version of the given discrete information measure, This is just the value of information divided its maximum possible value given o.

The same convenience syntaxes as in information can be used here.

Notice that there is no method information_normalized(e::DiscreteInfoEstimator, probs::Probabilities), because there is no way to know the number of possible outcomes (i.e., the total_outcomes) from probs.

Normalized values

For the PlugIn estimator, it is guaranteed that h̃ ∈ [0, 1]. For any other estimator, we can't guarantee this, since the estimator might over-correct. You should know what you're doing if using anything but PlugIn to estimate normalized values.

ComplexityMeasures.invariantmeasure — Method

invariantmeasure(x::AbstractStateSpaceSet, binning::RectangularBinning;
    rng = Random.default_rng()) → iv::InvariantMeasure

Estimate an invariant measure over the points in x based on binning the data into rectangular boxes dictated by the binning, then approximate the transfer (Perron-Frobenius) operator over the bins. From the approximation to the transfer operator, compute an invariant distribution over the bins. Assumes that the input data are sequential.

Details on the estimation procedure is found the TransferOperator docstring.

Example

using DynamicalSystems
henon_rule(x, p, n) = SVector{2}(1.0 - p[1]*x[1]^2 + x[2], p[2]*x[1])
henon = DeterministicIteratedMap(henon_rule, zeros(2), [1.4, 0.3])
orbit, t = trajectory(ds, 20_000; Ttr = 10)

# Estimate the invariant measure over some coarse graining of the orbit.
iv = invariantmeasure(orbit, RectangularBinning(15))

# Get the probabilities and bins
invariantmeasure(iv)

Probabilities and bin information

invariantmeasure(iv::InvariantMeasure) → (ρ::Probabilities, bins::Vector{<:SVector})

From a pre-computed invariant measure, return the probabilities and associated bins. The element ρ[i] is the probability of visitation to the box bins[i].

Transfer operator approach vs. naive histogram approach

Why bother with the transfer operator instead of using regular histograms to obtain probabilities?

In fact, the naive histogram approach and the transfer operator approach are equivalent in the limit of long enough time series (as $n \to \intfy$), which is guaranteed by the ergodic theorem. There is a crucial difference, however:

The naive histogram approach only gives the long-term probabilities that orbits visit a certain region of the state space. The transfer operator encodes that information too, but comes with the added benefit of knowing the transition probabilities between states (see transfermatrix).