Probabilities

ProbabilitiesEstimator

The supertype for all probabilities estimators.

In ComplexityMeasures.jl, probability distributions are estimated from data by defining a set of possible outcomes $\Omega = \{\omega_1, \omega_2, \ldots, \omega_L \}$, and assigning to each outcome $\omega_i$ a probability $p(\omega_i)$, such that $\sum_{i=1}^N p(\omega_i) = 1$. It is the role of a ProbabilitiesEstimator to

1. Define $\Omega$, the "outcome space", which is the set of all possible outcomes over which probabilities are estimated. The cardinality of this set can be obtained using total_outcomes.
2. Define how probabilities $p_i = p(\omega_i)$ are assigned to outcomes $\omega_i$.

In practice, probability estimation is done by calling probabilities with some input data and one of the following probabilities estimators. The result is a Probabilities p (Vector-like), where each element p[i] is the probability of the outcome ω[i]. Use probabilities_and_outcomes if you need both the probabilities and the outcomes, and use outcome_space to obtain $\Omega$ alone. The element type of $\Omega$ varies between estimators, but it is guranteed to be hashable. This allows for conveniently tracking the probability of a specific event across experimental realizations, by using the outcome as a dictionary key and the probability as the value for that key (or, alternatively, the key remains the outcome and one has a vector of probabilities, one for each experimental realization).

Some estimators can deduce $\Omega$ without knowledge of the input, such as SymbolicPermutation. For others, knowledge of input is necessary for concretely specifying $\Omega$, such as ValueHistogram with RectangularBinning. This only matters for the functions outcome_space and total_outcomes.

All currently implemented probability estimators are listed in a nice table in the probabilities estimators section of the online documentation.

probabilities(est::ProbabilitiesEstimator, x::Array_or_Dataset) → p::Probabilities

Compute a probability distribution over the set of possible outcomes defined by the probabilities estimator est, given input data x, which is typically an Array or a Dataset; see Input data for ComplexityMeasures.jl. Configuration options are always given as arguments to the chosen estimator.

To obtain the outcomes corresponding to these probabilities, use outcomes.

Due to performance optimizations, whether the returned probablities contain 0s as entries or not depends on the estimator. E.g., in ValueHistogram 0s are skipped, while in SymbolicPermutation 0 are not, because we get them for free.

probabilities(x::Array_or_Dataset) → p::Probabilities

Estimate probabilities by directly counting the elements of x, assuming that Ω = sort(unique(x)), i.e. that the outcome space is the unique elements of x. This is mostly useful when x contains categorical data.

See also: Probabilities, ProbabilitiesEstimator.

probabilities!(s, args...)

Similar to probabilities(args...), but allows pre-allocation of temporarily used containers s.

Only works for certain estimators. See for example SymbolicPermutation.

Probabilities <: AbstractVector
Probabilities(x) → p

Probabilities is a simple wrapper around x::AbstractVector{<:Real} that ensures its values sum to 1, so that p can be interpreted as probability mass function.

probabilities_and_outcomes(est, x)

Return probs, outs, where probs = probabilities(x, est) and outs[i] is the outcome with probability probs[i]. The element type of outs depends on the estimator. outs is a subset of the outcome_space of est.

See also outcomes, total_outcomes.

outcomes(est::ProbabilitiesEstimator, x)

Return all (unique) outcomes contained in x according to the given estimator. Equivalent with probabilities_and_outcomes(x, est)[2], but for some estimators it may be explicitly extended for better performance.

outcome_space(est::ProbabilitiesEstimator, x) → Ω

Return a container containing all possible outcomes of est for input x.

For some estimators the concrete outcome space is known without knowledge of input x, in which case the function dispatches to outcome_space(est). In general it is recommended to use the 2-argument version irrespectively of estimator.

total_outcomes(est::ProbabilitiesEstimator, x)

Return the length (cardinality) of the outcome space $\Omega$ of est.

For some estimators the concrete outcome space is known without knowledge of input x, in which case the function dispatches to total_outcomes(est). In general it is recommended to use the 2-argument version irrespectively of estimator.

missing_outcomes(est::ProbabilitiesEstimator, x) → n_missing::Int

Estimate a probability distribution for x using the given estimator, then count the number of missing (i.e. zero-probability) outcomes.

See also: MissingDispersionPatterns.

CountOccurrences()

A probabilities/entropy estimator 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.

ValueHistogram(b::AbstractBinning) <: ProbabilitiesEstimator

A probability estimator 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:

• RectangularBinning
• FixedRectangularBinning

The ValueHistogram 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.

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

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

Outcomes

The outcome space for ValueHistogram is the unique bins constructed from b. Each bin is identified by its left (lowest-value) corner. The bins are in data units, not integer (cartesian indices units), and are returned as SVectors. For FixedRectangularBinning this is well-defined from the binning, but for RectangularBinning input x is needed for a well-defined outcome_space.

RectangularBinning(ϵ) <: AbstractBinning

Rectangular box partition of state space using the scheme ϵ, deducing the coordinates of the grid axis minima from the input data. Generally it is preferred to use FixedRectangularBinning instead, as it has a well defined outcome space without knowledge of input data.

Binning instructions are deduced from the type of ϵ as follows:

1. ϵ::Int divides each coordinate axis into ϵ equal-length intervals that cover all data.
2. ϵ::Float64 divides each coordinate axis into intervals of fixed size ϵ, starting from the axis minima until the data is completely covered by boxes.
3. ϵ::Vector{Int} divides the i-th coordinate axis into ϵ[i] equal-length intervals that cover all data.
4. ϵ::Vector{Float64} divides the i-th coordinate axis into intervals of fixed size ϵ[i], starting from the axis minima until the data is completely covered by boxes.

FixedRectangularBinning <: AbstractBinning
FixedRectangularBinning(ϵmin::NTuple, ϵmax::NTuple, N::Int)

Rectangular box partition of state space where the extent of the grid is explicitly specified by ϵmin and emax, and along each dimension, the grid is subdivided into N subintervals. Points falling outside the partition do not contribute to probabilities. This binning type leads to a well-defined outcome space without knowledge of input data, see ValueHistogram.

ϵmin/emax must be NTuple{D, <:Real} for input of D-dimensional data.

FixedRectangularBinning(ϵmin::Real, ϵmax::Real, N::Int, D::Int = 1)

This is a convenience method where each dimension of the binning has the same extent and the input data are D dimensional, which defaults to 1 (timeseries).

SymbolicPermutation <: ProbabilitiesEstimator
SymbolicPermutation(; m = 3, τ = 1, lt::Function = ComplexityMeasures.isless_rand)

A probabilities estimator based on ordinal permutation patterns.

When passed to probabilities the output depends on the input data type:

• Univariate data. If applied to a univariate timeseries (Vector), then the timeseries is first embedded using embedding delay τ and dimension m, resulting in embedding vectors $\{ \bf{x}_i \}_{i=1}^{N-(m-1)\tau}$. Then, for each $\bf{x}_i$, we find its permutation pattern $\pi_{i}$. Probabilities are then estimated as the frequencies of the encoded permutation symbols by using CountOccurrences. When giving the resulting probabilities to entropy, the original permutation entropy is computed ^{[BandtPompe2002]}.
• Multivariate data. If applied to a an D-dimensional Dataset, then no embedding is constructed, and we each vector $\bf{x}_i$ of the dataset directly to its permutation pattern $\pi_{i}$, $\pi_{i}$ by comparing the relative magnitudes of the elements of $\bf{x}_i$. Like above, probabilities are estimated as the frequencies of the permutation symbols. In this case, m is ignored, but m must still match the dimension of the dataset for optimization. The resulting probabilities can be used to compute multivariate permutation entropy^[He2016], although here we don't perform any further subdivision of the permutation patterns (as in Figure 3 of^[He2016]).

Internally, SymbolicPermutation uses the OrdinalPatternEncoding to represent ordinal patterns as integers for efficient computations.

Outcome space

The outcome space Ω for SymbolicPermutation is the set of length-m ordinal patterns (i.e. permutations) that can be formed by the integers 1, 2, …, m, ordered lexicographically. There are factorial(m) such patterns.

For example, the outcome [3, 1, 2] corresponds to the ordinal pattern of having first the largest value, then the lowest value, and then the value in between.

In-place symbolization

SymbolicPermutation also implements the in-place probabilities! for Dataset input (or embedded vector input). The length of the pre-allocated symbol vector must match the length of the dataset. For example

using DelayEmbeddings, ComplexityMeasures
m, N = 2, 100
est = SymbolicPermutation(; m, τ)
x = Dataset(rand(N, m) # timeseries example
πs_ts = zeros(Int, N) # length must match length of `x`
p = probabilities!(πs_ts, est, x)

See SymbolicWeightedPermutation and SymbolicAmplitudeAwarePermutation for estimators that not only consider ordinal (sorting) patterns, but also incorporate information about within-state-vector amplitudes. For a version of this estimator that can be used on high-dimensional arrays, see SpatialSymbolicPermutation.

SymbolicWeightedPermutation <: ProbabilitiesEstimator
SymbolicWeightedPermutation(; τ = 1, m = 3, lt::Function = ComplexityMeasures.isless_rand)

A variant of SymbolicPermutation that also incorporates amplitude information, based on the weighted permutation entropy^{[Fadlallah2013]}. The outcome space and keywords are the same as in SymbolicPermutation.

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.

SymbolicAmplitudeAwarePermutation <: ProbabilitiesEstimator
SymbolicAmplitudeAwarePermutation(; τ = 1, m = 3, A = 0.5, lt = ComplexityMeasures.isless_rand)

A variant of SymbolicPermutation that also incorporates amplitude information, based on the amplitude-aware permutation entropy^[Azami2016]. The outcome space and keywords are the same as in SymbolicPermutation.

Description

Similarly to SymbolicWeightedPermutation, 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.

Dispersion(; c = 5, m = 2, τ = 1, check_unique = true)

A probability estimator based on dispersion patterns, originally used by Rostaghi & Azami, 2016^{[Rostaghi2016]} to compute the "dispersion entropy", which characterizes the complexity and irregularity of a time series.

Recommended parameter values^[Li2018] are m ∈ [2, 3], τ = 1 for the embedding, and c ∈ [3, 4, …, 8] categories for the Gaussian symbol mapping.

Description

Assume we have a univariate time series $X = \{x_i\}_{i=1}^N$. First, this time series is encoded into a symbol timeseries $S$ using the Gaussian encoding GaussianCDFEncoding with empirical mean μ and empirical standard deviation σ (both determined from $X$), and c as given to Dispersion.

Then, $S$ is embedded into an $m$-dimensional time series, using an embedding lag of $\tau$, which yields a total of $N - (m - 1)\tau$ delay vectors $z_i$, or "dispersion patterns". Since each element of $z_i$ can take on c different values, and each delay vector has m entries, there are c^m possible dispersion patterns. This number is used for normalization when computing dispersion entropy.

The returned probabilities are simply the frequencies of the unique dispersion patterns present in $S$ (i.e., the CountOccurences of $S$).

Outcome space

The outcome space for Dispersion is the unique delay vectors whose elements are the the symbols (integers) encoded by the Gaussian CDF, i.e., the unique elements of $S$.

Data requirements and parameters

The input must have more than one unique element for the Gaussian mapping to be well-defined. Li et al. (2018) recommends that x has at least 1000 data points.

If check_unique == true (default), then it is checked that the input has more than one unique value. If check_unique == false and the input only has one unique element, then a InexactError is thrown when trying to compute probabilities.

For a version of this estimator that can be used on high-dimensional arrays, see SpatialDispersion.

TransferOperator <: ProbabilitiesEstimator
TransferOperator(b::RectangularBinning)

A probability estimator based on binning data into rectangular boxes dictated by the given binning scheme b, then approximating the transfer (Perron-Frobenius) operator over the bins, then taking the invariant measure associated with that transfer operator as the bin probabilities. Assumes that the input data are sequential (time-ordered).

This implementation follows the grid estimator approach in Diego et al. (2019)^[Diego2019].

Outcome space

The outcome space for TransferOperator is the set of unique bins constructed from b. Bins are identified by their left (lowest-value) corners, are given in data units, and are returned as SVectors.

Bin ordering

Bins returned by probabilities_and_outcomes are ordered according to first appearance (i.e. the first time the input (multivariate) timeseries visits the bin). Thus, if

b = RectangularBinning(4)
est = TransferOperator(b)
probs, outcomes = probabilities_and_outcomes(x, est) # x is some timeseries

then probs[i] is the invariant measure (probability) of the bin outcomes[i], which is the i-th bin visited by the timeseries with nonzero measure.

Description

The transfer operator $P^{N}$is computed as an N-by-N matrix of transition probabilities between the states defined by the partition elements, where N is the number of boxes in the partition that is visited by the orbit/points.

If $\{x_t^{(D)} \}_{n=1}^L$ are the $L$ different $D$-dimensional points over which the transfer operator is approximated, $\{ C_{k=1}^N \}$ are the $N$ different partition elements (as dictated by ϵ) that gets visited by the points, and $\phi(x_t) = x_{t+1}$, then

\[P_{ij} = \dfrac {\#\{ x_n | \phi(x_n) \in C_j \cap x_n \in C_i \}} {\#\{ x_m | x_m \in C_i \}},\]

where $\#$ denotes the cardinal. The element $P_{ij}$ thus indicates how many points that are initially in box $C_i$ end up in box $C_j$ when the points in $C_i$ are projected one step forward in time. Thus, the row $P_{ik}^N$ where $k \in \{1, 2, \ldots, N \}$ gives the probability of jumping from the state defined by box $C_i$ to any of the other $N$ states. It follows that $\sum_{k=1}^{N} P_{ik} = 1$ for all $i$. Thus, $P^N$ is a row/right stochastic matrix.

Invariant measure estimation from transfer operator

The left invariant distribution $\mathbf{\rho}^N$ is a row vector, where $\mathbf{\rho}^N P^{N} = \mathbf{\rho}^N$. Hence, $\mathbf{\rho}^N$ is a row eigenvector of the transfer matrix $P^{N}$ associated with eigenvalue 1. The distribution $\mathbf{\rho}^N$ approximates the invariant density of the system subject to binning, and can be taken as a probability distribution over the partition elements.

In practice, the invariant measure $\mathbf{\rho}^N$ is computed using invariantmeasure, which also approximates the transfer matrix. The invariant distribution is initialized as a length-N random distribution which is then applied to $P^{N}$. The resulting length-N distribution is then applied to $P^{N}$ again. This process repeats until the difference between the distributions over consecutive iterations is below some threshold.

See also: RectangularBinning, invariantmeasure.

InvariantMeasure(to, ρ)

Minimal return struct for invariantmeasure that contains the estimated invariant measure ρ, as well as the transfer operator to from which it is computed (including bin information).

See also: invariantmeasure.

invariantmeasure(x::AbstractDataset, binning::RectangularBinning) → 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, Plots, ComplexityMeasures
D = 4
ds = Systems.lorenz96(D; F = 32.0)
N, dt = 20000, 0.1
orbit = trajectory(ds, N*dt; dt = dt, Ttr = 10.0)

# 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]. Analogous to binhist.

See also: InvariantMeasure.

transfermatrix(iv::InvariantMeasure) → (M::AbstractArray{<:Real, 2}, bins::Vector{<:SVector})

Return the transfer matrix/operator and corresponding bins. Here, bins[i] corresponds to the i-th row/column of the transfer matrix. Thus, the entry M[i, j] is the probability of jumping from the state defined by bins[i] to the state defined by bins[j].

See also: TransferOperator.

NaiveKernel(ϵ::Real; method = KDTree, w = 0, metric = Euclidean()) <: ProbabilitiesEstimator

Estimate probabilities/entropy using a "naive" kernel density estimation approach (KDE), as discussed in Prichard and Theiler (1995) ^{[PrichardTheiler1995]}.

Probabilities $P(\mathbf{x}, \epsilon)$ are assigned to every point $\mathbf{x}$ by counting how many other points occupy the space spanned by a hypersphere of radius ϵ around $\mathbf{x}$, according to:

\[P_i( X, \epsilon) \approx \dfrac{1}{N} \sum_{s} B(||X_i - X_j|| < \epsilon),\]

where $B$ gives 1 if the argument is true. Probabilities are then normalized.

Keyword arguments

• method = KDTree: the search structure supported by Neighborhood.jl. Specifically, use KDTree to use a tree-based neighbor search, or BruteForce for the direct distances between all points. KDTrees heavily outperform direct distances when the dimensionality of the data is much smaller than the data length.
• w = 0: the Theiler window, which excludes indices $s$ that are within $|i - s| ≤ w$ from the given point $x_i$.
• metric = Euclidean(): the distance metric.

Outcome space

The outcome space Ω for NaiveKernel are the indices of the input data, eachindex(x). Hence, input x is needed for a well-defined outcome_space. The reason to not return the data points themselves is because duplicate data points may not get assigned same probabilities (due to having different neighbors).

WaveletOverlap([wavelet]) <: ProbabilitiesEstimator

Apply the maximal overlap discrete wavelet transform (MODWT) to a signal, then compute probabilities as the (normalized) energies at different wavelet scales. These probabilities are used to compute the wavelet entropy, according to Rosso et al. (2001)^[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.

PowerSpectrum() <: ProbabilitiesEstimator

Calculate the power spectrum of a timeseries (amplitude square of its Fourier transform), and return the spectrum normalized to sum = 1 as probabilities. The Shannon entropy of these probabilities is typically referred in the literature as spectral entropy, e.g. ^[Llanos2016],^[Tian2017].

The closer the spectrum is to flat, i.e., white noise, the higher the entropy. However, you can't compare entropies of timeseries with different length, because the binning in spectral space depends on the length of the input.

Outcome space

The outcome space Ω for PowerSpectrum is the set of frequencies in Fourier space. They should be multiplied with the sampling rate of the signal, which is assumed to be 1. Input x is needed for a well-defined outcome_space.

Diversity(; m::Int, τ::Int, nbins::Int)

A ProbabilitiesEstimator based on the cosine similarity. It can be used with entropy to compute the diversity entropy of an input timeseries^[Wang2020].

The implementation here allows for τ != 1, which was not considered in the original paper.

Description

Diversity probabilities are computed as follows.

1. From the input time series x, using embedding lag τ and embedding dimension m, construct the embedding $Y = \{\bf x_i \} = \{(x_{i}, x_{i+\tau}, x_{i+2\tau}, \ldots, x_{i+m\tau - 1}\}_{i = 1}^{N-mτ}$.
2. Compute $D = \{d(\bf x_t, \bf x_{t+1}) \}_{t=1}^{N-mτ-1}$, where $d(\cdot, \cdot)$ is the cosine similarity between two m-dimensional vectors in the embedding.
3. Divide the interval [-1, 1] into nbins equally sized subintervals.
4. Construct a histogram of cosine similarities $d \in D$ over those subintervals.
5. Sum-normalize the histogram to obtain probabilities.

Outcome space

The outcome space for Diversity is the bins of the [-1, 1] interval, and the return configuration is the same as in ValueHistogram (left bin edge).

SpatialSymbolicPermutation <: ProbabilitiesEstimator
SpatialSymbolicPermutation(stencil, x; periodic = true)

A symbolic, permutation-based probabilities estimator for spatiotemporal systems that generalises SymbolicPermutation to high-dimensional arrays. The order m of the permutation pattern is extracted from the stencil, see below.

SpatialSymbolicPermutation is based on the 2D and 3D spatiotemporal permutation entropy estimators by by Ribeiro et al. (2012)^{[Ribeiro2012]} and Schlemmer et al. (2018)^{[Schlemmer2018]}), respectively, but is here implemented as a pure probabilities probabilities estimator that is generalized for D-dimensional input array x, with arbitrary regions (stencils) to get patterns form and (possibly) periodic boundary conditions.

See below for ways to specify the stencil. If periodic = true, then the stencil wraps around at the ends of the array. If false, then collected regions with indices which exceed the array bounds are skipped.

In combination with entropy and entropy_normalized, this probabilities estimator can be used to compute generalized spatiotemporal permutation EntropyDefinition of any type.

Outcome space

The outcome space Ω for SpatialSymbolicPermutation is the set of length-m ordinal patterns (i.e. permutations) that can be formed by the integers 1, 2, …, m, ordered lexicographically. There are factorial(m) such patterns. Here m refers to the number of points included in stencil.

Stencils

The stencil defines what local area to use to group hypervoxels. Each grouping of hypervoxels is mapped to an order-m permutation pattern, which is then mapped to an integer as in SymbolicPermutation. The stencil is moved around the input array, in a sense "scanning" the input array, to collect all possible groupings allowed by the boundary condition (periodic or not).

Stencils are passed in one of the following three ways:

1. As vectors of CartesianIndex which encode the offset of indices to include in the stencil, with respect to the current array index when scanning over the array. For example stencil = CartesianIndex.([(0,0), (0,1), (1,1), (1,0)]). Don't forget to include the zero offset index if you want to include the hypervoxel itself, which is almost always the case. Here the stencil creates a 2x2 square extending to the bottom and right of the pixel (directions here correspond to the way Julia prints matrices by default). When passing a stencil as a vector of CartesianIndex, m = length(stencil).
2. As a D-dimensional array (where D matches the dimensionality of the input data) containing 0s and 1s, where if stencil[index] == 1, the corresponding pixel is included, and if stencil[index] == 0, it is not included. To generate the same estimator as in 1., use stencil = [1 1; 1 1]. When passing a stencil as a D-dimensional array, m = sum(stencil)
3. As a Tuple containing two Tuples, both of length D, for D-dimensional data. The first tuple specifies the extent of the stencil, where extent[i] dictates the number of hypervoxels to be included along the ith axis and lag[i] the separation of hypervoxels along the same axis. This method can only generate (hyper)rectangular stencils. To create the same estimator as in the previous examples, use here stencil = ((2, 2), (1, 1)). When passing a stencil using extent and lag, m = prod(extent).

SpatialDispersion <: ProbabilitiesEstimator
SpatialDispersion(stencil, x::AbstractArray;
    periodic = true,
    c = 5,
    skip_encoding = false,
    L = nothing,
)

A dispersion-based probabilities estimator that generalises Dispersion for input data that are high-dimensional arrays.

SpatialDispersion is based on Azami et al. (2019)^[Azami2019]'s 2D square dispersion (Shannon) entropy estimator, but is here implemented as a pure probabilities probabilities estimator that is generalized for N-dimensional input data x, with arbitrary neighborhood regions (stencils) and (optionally) periodic boundary conditions.

In combination with entropy and entropy_normalized, this probabilities estimator can be used to compute (normalized) generalized spatiotemporal dispersion EntropyDefinition of any type.

Arguments

• stencil. Defines what local area (hyperrectangle), or which points within this area, to include around each hypervoxel (i.e. pixel in 2D). The examples below demonstrate different ways of specifying stencils. For details, see SpatialSymbolicPermutation. See SpatialSymbolicPermutation for more information about stencils.
• x::AbstractArray. The input data. Must be provided because we need to know its size for optimization and bound checking.

Keyword arguments

• periodic::Bool. If periodic == true, then the stencil should wrap around at the end of the array. If periodic = false, then pixels whose stencil exceeds the array bounds are skipped.
• c::Int. Determines how many discrete categories to use for the Gaussian encoding.
• skip_encoding. If skip_encoding == true, encoding is ignored, and dispersion patterns are computed directly from x, under the assumption that L is the alphabet length for x (useful for categorical or integer data). Thus, if skip_encoding == true, then L must also be specified. This is useful for categorical or integer-valued data.
• L. If L == nothing (default), then the number of total outcomes is inferred from stencil and encoding. If L is set to an integer, then the data is considered pre-encoded and the number of total outcomes is set to L.

Outcome space

The outcome space for SpatialDispersion is the unique delay vectors whose elements are the the symbols (integers) encoded by the Gaussian CDF. Hence, the outcome space is all m-dimensional delay vectors whose elements are all possible values in 1:c. There are c^m such vectors.

Description

Estimating probabilities/entropies from higher-dimensional data is conceptually simple.

1. Discretize each value (hypervoxel) in x relative to all other values xᵢ ∈ x using the provided encoding scheme.
2. Use stencil to extract relevant (discretized) points around each hypervoxel.
3. Construct a symbol these points.
4. Take the sum-normalized histogram of the symbol as a probability distribution.
5. Optionally, compute entropy or entropy_normalized from this probability distribution.

Usage

Here's how to compute spatial dispersion entropy using the three different ways of specifying stencils.

x = rand(50, 50) # first "time slice" of a spatial system evolution

# Cartesian stencil
stencil_cartesian = CartesianIndex.([(0,0), (1,0), (1,1), (0,1)])
est = SpatialDispersion(stencil_cartesian, x)
entropy_normalized(est, x)

# Extent/lag stencil
extent = (2, 2); lag = (1, 1); stencil_ext_lag = (extent, lag)
est = SpatialDispersion(stencil_ext_lag, x)
entropy_normalized(est, x)

# Matrix stencil
stencil_matrix = [1 1; 1 1]
est = SpatialDispersion(stencil_matrix, x)
entropy_normalized(est, x)

To apply this to timeseries of spatial data, simply loop over the call (broadcast), e.g.:

imgs = [rand(50, 50) for i = 1:100]; # one image per second over 100 seconds
stencil = ((2, 2), (1, 1)) # a 2x2 stencil (i.e. dispersion patterns of length 4)
est = SpatialDispersion(stencil, first(imgs))
h_vs_t = entropy_normalized.(Ref(est), imgs)

Computing generalized spatiotemporal dispersion entropy is trivial, e.g. with Renyi:

x = reshape(repeat(1:5, 500) .+ 0.1*rand(500*5), 50, 50)
est = SpatialDispersion(stencil, x)
entropy(Renyi(q = 2), est, x)

See also: SpatialSymbolicPermutation, GaussianCDFEncoding, symbolize.