Dispersion and reverse dispersion probability estimators (#96)

* Dispersion and reverse dispersion probability estimators

* Fix tests

* ReverseDispersion should be a complexity measure, plus addressing some comments

* Export reverse_dispersion

* Reverse dispersion

* Remove reference to `ReverseDispersion`

* Update docstring

* Add analytical tests

* Improve tests and add doc examples

* Fix tests

* Better docs and doctests

* Remove file

* Fix tests

* Update src/complexity_measures/reverse_dispersion_entropy.jl

Co-authored-by: George Datseris <[email protected]>

* Addressing review comments (#99)

* Address review comments

* Throw error

* Remove non-used field

* Fix tests

* Merge docs and move individual methods to their respective files

* Tsallis reduces to Shannon entropy for q -> 1.

* Normalized entropy API, including utility functions

* Analytical test cases for Tsallis

* Include `k` in `maxentropy_tsallis`

Co-authored-by: George Datseris <[email protected]>

Author: kahaaga

docs/src/complexity_measures.md

@@ -1,7 +1,14 @@
-# Complexity measures
+# [Complexity measures](@id complexity_measures)
 
 ## Sample entropy
 
 ## Approximate entropy
 
+## Reverse dispersion entropy
+
+```@docs
+reverse_dispersion
+distance_to_whitenoise
+```
+
 ## Disequilibrium

docs/src/entropies.md

@@ -13,10 +13,21 @@ entropy_tsallis
 ```
 
 ## Shannon entropy (convenience)
+
 ```@docs
 entropy_shannon
 ```
 
+## Normalization
+
+The generic [`entropy_normalized`](@ref) normalizes any entropy value to the entropy of a
+uniform distribution. We also provide [maximum entropy](@ref maximum_entropy) functions
+that are useful for manual normalization.
+
+```@docs
+entropy_normalized
+```
+
 ## Indirect entropies
 Here we list functions which compute Shannon entropies via alternate means, without explicitly computing some probability distributions and then using the Shannon formulat.
 
@@ -33,4 +44,4 @@ entropy_permutation
 entropy_spatial_permutation
 entropy_wavelet
 entropy_dispersion
-```
+```

docs/src/examples.md

@@ -145,3 +145,60 @@ for a in (ax, ay, az); axislegend(a); end
 for a in (ax, ay); hidexdecorations!(a; grid=false); end
 fig
 ```
+
+## [Dispersion and reverse dispersion entropy](@id dispersion_examples)
+
+Here we reproduce parts of figure 3 in Li et al. (2019), computing reverse and regular dispersion entropy for a time series consisting of normally distributed noise with a single spike in the middle of the signal. We compute the entropies over a range subsets of the data, using a sliding window consisting of 70 data points, stepping the window 10 time steps at a time.
+
+Note: the results here are not exactly the same as in the original paper, because Li et 
+al. (2019) base their examples on randomly generated numbers and do not provide code that 
+specify random number seeds.
+
+```@example
+using Entropies, DynamicalSystems, Random, CairoMakie, Distributions
+
+n = 1000
+ts = 1:n
+x = [i == n ÷ 2 ? 50.0 : 0.0 for i in ts]
+rng = Random.default_rng()
+s = rand(rng, Normal(0, 1), n)
+y = x .+ s
+
+ws = 70
+windows = [t:t+ws for t in 1:10:n-ws]
+rdes = zeros(length(windows))
+des = zeros(length(windows))
+pes = zeros(length(windows))
+
+m, c = 2, 6
+est_de = Dispersion(symbolization = GaussianSymbolization(c), m = m, τ = 1)
+
+for (i, window) in enumerate(windows)
+    rdes[i] = reverse_dispersion(y[window], est_de; normalize = true)
+    des[i] = entropy_renyi_norm(y[window], est_de)
+end
+
+fig = Figure()
+
+a1 = Axis(fig[1,1]; xlabel = "Time step", ylabel = "Value")
+lines!(a1, ts, y)
+display(fig)
+
+a2 = Axis(fig[2, 1]; xlabel = "Time step", ylabel = "Value")
+p_rde = scatterlines!([first(w) for w in windows], rdes,
+    label = "Reverse dispersion entropy",
+    color = :black,
+    markercolor = :black, marker = '●')
+p_de = scatterlines!([first(w) for w in windows], des,
+    label = "Dispersion entropy",
+    color = :red,
+    markercolor = :red, marker = 'x', markersize = 20)
+
+axislegend(position = :rc)
+ylims!(0, max(maximum(pes), 1))
+fig
+```
+
+[^Rostaghi2016]: Rostaghi, M., & Azami, H. (2016). Dispersion entropy: A measure for time-series analysis. IEEE Signal Processing Letters, 23(5), 610-614.
+[^Li2019]: Li, Y., Gao, X., & Wang, L. (2019). Reverse dispersion entropy: a new
+    complexity measure for sensor signal. Sensors, 19(23), 5203.

docs/src/index.md

@@ -42,7 +42,7 @@ Thus, any of the implemented [probabilities estimators](@ref estimators) can be
 
 
 ### Complexity measures
-Other complexity measures, which strictly speaking don't compute entropies, and may or may not explicitly compute probability distributions, appear in the [Complexity measures](@ref) section.
+Other complexity measures, which strictly speaking don't compute entropies, and may or may not explicitly compute probability distributions, appear in the [Complexity measures](@ref complexity_measures) section.
 
 
 ## Input data

docs/src/probabilities.md

@@ -24,6 +24,12 @@ SymbolicPermutation
 SpatialSymbolicPermutation
 ```
 
+## Dispersion (symbolic)
+
+```@docs
+Dispersion
+```
+
 ## Visitation frequency (binning)
 
 ```@docs

docs/src/utils.md

@@ -26,3 +26,17 @@ OrdinalPattern
 ```@docs
 Entropies.encode_motif
 ```
+
+### Normalization
+
+```@docs
+alphabet_length
+```
+
+#### [Maximum entropy](@id maximum_entropy)
+
+```@docs
+maxentropy_tsallis
+maxentropy_renyi
+maxentropy_shannon
+```

src/Entropies.jl

@@ -16,6 +16,8 @@ include("symbolization/symbolize.jl")
 include("probabilities.jl")
 include("probabilities_estimators/probabilities_estimators.jl")
 include("entropies/entropies.jl")
+include("complexity_measures/complexity_measures.jl")
+
 include("deprecations.jl")
 
 

src/complexity_measures/complexity_measures.jl

@@ -0,0 +1 @@
+include("reverse_dispersion_entropy.jl")

src/complexity_measures/reverse_dispersion_entropy.jl

@@ -0,0 +1,76 @@
+export reverse_dispersion
+export distance_to_whitenoise
+
+# Note: this is not an entropy estimator, so we don't use the entropy_xxx_norm interface
+# for normalization, even though we rely on `alphabet_length`.
+"""
+    distance_to_whitenoise(p::Probabilities, estimator::Dispersion; 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([`reverse_dispersion`](@ref);
+Li et al., 2019[^Li2019]).
+
+[^Li2019]: Li, Y., Gao, X., & Wang, L. (2019). Reverse dispersion entropy: a new
+    complexity measure for sensor signal. Sensors, 19(23), 5203.
+"""
+function distance_to_whitenoise(p::Probabilities, est::Dispersion; normalize = false)
+    # We can safely skip non-occurring symbols, because they don't contribute
+    # to the sum in eq. 3 in Li et al. (2019)
+    Hrde = sum(abs2, p) - (1 / alphabet_length(est))
+
+    if normalize
+        return Hrde / (1 - (1 / alphabet_length(est)))
+    else
+        return Hrde
+    end
+end
+
+# Note again: this is a *complexity measure*, not an entropy estimator, so we don't use
+# the entropy_xxx_norm interface for normalization, even though we rely on `alphabet_length`.
+"""
+    reverse_dispersion(x::AbstractVector{T}, est::Dispersion = Dispersion();
+        normalize = true) where T <: Real
+
+Compute the reverse dispersion entropy complexity measure (Li et al., 2019)[^Li2019].
+
+## Description
+
+Li et al. (2021)[^Li2019] defines the reverse dispersion entropy as
+
+```math
+H_{rde} = \\sum_{i = 1}^{c^m} \\left(p_i - \\dfrac{1}{{c^m}} \\right)^2 =
+\\left( \\sum_{i=1}^{c^m} p_i^2 \\right) - \\dfrac{1}{c^{m}}
+```
+where the probabilities ``p_i`` are obtained precisely as for the [`Dispersion`](@ref)
+probability estimator. Relative frequencies of dispersion patterns are computed using the
+given `symbolization` scheme , which defaults to symbolization using the normal cumulative
+distribution function (NCDF), as implemented by [`GaussianSymbolization`](@ref), using
+embedding dimension `m` and embedding delay `τ`.
+Recommended parameter values[^Li2018] are `m ∈ [2, 3]`, `τ = 1` for the embedding, and
+`c ∈ [3, 4, …, 8]` categories for the Gaussian mapping.
+
+If `normalize == true`, then the reverse dispersion entropy is normalized to `[0, 1]`.
+
+The minimum value of ``H_{rde}`` is zero and occurs precisely when the dispersion
+pattern distribution is flat, which occurs when all ``p_i``s are equal to ``1/c^m``.
+Because ``H_{rde} \\geq 0``, ``H_{rde}`` can therefore be said to be a measure of how far
+the dispersion pattern probability distribution is from white noise.
+
+[^Li2019]: Li, Y., Gao, X., & Wang, L. (2019). Reverse dispersion entropy: a new
+    complexity measure for sensor signal. Sensors, 19(23), 5203.
+"""
+function reverse_dispersion(x::AbstractVector{T}, est::Dispersion = Dispersion();
+        normalize = true) where T <: Real
+
+    p = probabilities(x, est)
+
+    # The following step combines distance information with the probabilities, so
+    # from here on, it is not possible to use `renyi_entropy` or similar methods, because
+    # we're not dealing with probabilities anymore.
+    Hrde = distance_to_whitenoise(p, est, normalize = normalize)
+end

src/entropies/convenience_definitions.jl

@@ -53,6 +53,20 @@ function entropy_wavelet(x; wavelet = Wavelets.WT.Daubechies{12}(), base = MathC
     entropy_renyi(x, est; base, q = 1)
 end
 
-function entropy_dispersion(args...)
+"""
+    entropy_dispersion(x; m = 2, τ = 1, s = GaussianSymbolization(3),
+        base = MathConstants.e)
+
+Compute the dispersion entropy. This function is just a convenience call to:
+```julia
+est = Dispersion(m = m, τ = τ, s = s)
+entropy_shannon(x, est; base)
+```
+See [`Dispersion`](@ref) for more info.
+"""
+function entropy_dispersion(x; wavelet = Wavelets.WT.Daubechies{12}(),
+        base = MathConstants.e)
 
-end
+    est = Dispersion(m = m, τ = τ, s = s)
+    entropy_renyi(x, est; base, q = 1)
+end