diff --git a/Project.toml b/Project.toml index 71474f32..16d020c9 100644 --- a/Project.toml +++ b/Project.toml @@ -2,7 +2,7 @@ name = "Associations" uuid = "614afb3a-e278-4863-8805-9959372b9ec2" authors = ["Kristian Agasøster Haaga ", "Tor Einar Møller ", "George Datseris "] repo = "https://github.com/kahaaga/Associations.jl.git" -version = "4.3.0" +version = "4.4.0" [deps] Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697" diff --git a/README.md b/README.md index 25dcb008..668e6b81 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ # Associations -[![CI](https://github.com/juliadynamics/Associations.jl/workflows/CI/badge.svg)](https://github.com/JuliaDynamics/Associations.jl/actions) +[![CI (main)](https://github.com/juliadynamics/Associations.jl/workflows/CI/badge.svg?branch=main)](https://github.com/JuliaDynamics/Associations.jl/actions) [![](https://img.shields.io/badge/docs-latest_tagged-blue.svg)](https://juliadynamics.github.io/Associations.jl/stable/) [![](https://img.shields.io/badge/docs-dev_(main)-blue.svg)](https://juliadynamics.github.io/Associations.jl/dev/) [![codecov](https://codecov.io/gh/JuliaDynamics/Associations.jl/branch/main/graph/badge.svg?token=0b71n6x6AP)](https://codecov.io/gh/JuliaDynamics/Associations.jl) diff --git a/changelog.md b/changelog.md index 49225b64..f3ce3c99 100644 --- a/changelog.md +++ b/changelog.md @@ -2,6 +2,11 @@ From version v4.0 onwards, this package has been renamed to to Associations.jl. +# 4.4 + +- New association measure: `SECMI` (`ShortExpansionConditionalMutualInformation`) +- New independence test: `SECMITest`, which is based on `SECMI`. + # 4.3 - Compatiblity with StateSpaceSets.jl v2.X diff --git a/docs/refs.bib b/docs/refs.bib index 795cd0d0..02987fd5 100644 --- a/docs/refs.bib +++ b/docs/refs.bib @@ -1333,4 +1333,14 @@ @article{Azadkia2021 pages={3070--3102}, year={2021}, publisher={Institute of Mathematical Statistics} +} + +@article{Kubkowski2021, + title={How to gain on power: novel conditional independence tests based on short expansion of conditional mutual information}, + author={Kubkowski, Mariusz and Mielniczuk, Jan and Teisseyre, Pawe{\l}}, + journal={Journal of Machine Learning Research}, + volume={22}, + number={62}, + pages={1--57}, + year={2021} } \ No newline at end of file diff --git a/docs/src/associations.md b/docs/src/associations.md index ad690522..8143384f 100644 --- a/docs/src/associations.md +++ b/docs/src/associations.md @@ -107,6 +107,12 @@ EmbeddingTE PartialMutualInformation ``` +### Short expansion of conditional mutual information + +```@docs +ShortExpansionConditionalMutualInformation +``` + ## [Correlation measures](@id correlation_api) ```@docs diff --git a/docs/src/examples/examples_associations.md b/docs/src/examples/examples_associations.md index e584a0e0..887ecc91 100644 --- a/docs/src/examples/examples_associations.md +++ b/docs/src/examples/examples_associations.md @@ -1093,6 +1093,24 @@ est = MIDecomposition(CMIShannon(base = 2), KSG1(k = 10)) association(est, x, z, y) ``` +## [`ShortExpansionConditionalMutualInformation`](@ref) + +### [[`JointProbabilities`](@ref) with [`CodifyVariables`](@ref) and [`ValueBinning`](@ref)](@id example_ShortExpansionConditionalMutualInformation_JointProbabilities_CodifyVariables_ValueBinning) + +```@example +using Associations +using Test +using Random; rng = Xoshiro(1234) +n = 20 +x = rand(rng, n) +y = randn(rng, n) .+ x .^ 2 +z = randn(rng, n) .* y + +# An estimator for estimating the SECMI measure +est = JointProbabilities(SECMI(base = 2), CodifyVariables(ValueBinning(3))) +association(est, x, z, y) +``` + ### [[`EntropyDecomposition`](@ref) + [`Kraskov`](@ref)](@id example_CMIShannon_EntropyDecomposition_Kraskov) Any [`DifferentialInfoEstimator`](@ref) can also be used to compute conditional diff --git a/docs/src/examples/examples_independence.md b/docs/src/examples/examples_independence.md index ae579245..c29986f9 100644 --- a/docs/src/examples/examples_independence.md +++ b/docs/src/examples/examples_independence.md @@ -469,4 +469,54 @@ connecting `x` and `z`.) independence(test, x, z, y) ``` -The test verifies our expectation. \ No newline at end of file +The test verifies our expectation. +## [[`SECMITest`](@ref)](@id example_SECMITEST) + +## [[`JointProbabilities`](@ref) estimation on numeric data](@id example_SECMITEST_JointProbabilities_CodifyVariables_ValueBinning) + +```@example example_SECMITEst +using Associations +using Test +using Random; rng = Xoshiro(1234) +n = 25 +x = rand(rng, n) +y = randn(rng, n) .+ x .^ 2 +z = randn(rng, n) .* y + +# An estimator for estimating the SECMI measure +est = JointProbabilities(SECMI(base = 2), CodifyVariables(ValueBinning(3))) +test = SECMITest(est; nshuffles = 19) +``` + +When analyzing ``SECMI(x, y | z)``, the expectation is to reject the null hypothesis (independence), since `x` and `y` are connected, regardless of the effect of `z`. + +```@example example_SECMITEst +independence(test, x, y, z) +``` + +We can detect this association, even for `n = 25`! When analyzing ``SECMI(x, z | y)``, we +expect that we can't reject the null (indepdendence), precisely since `x` and `z` are *not* +connected when "conditioning away" `y`. + +```@example example_SECMITEst +independence(test, x, z, y) +``` + +## [[`JointProbabilities`](@ref) estimation on categorical data](@id example_SECMITEST_JointProbabilities_CodifyVariables_UniqueElements) + +Note that this also works for categorical variables. Just use [`UniqueElements`](@ref) to +discretize! + +```@example example_SECMITest_categorical +using Associations +using Test +using Random; rng = Xoshiro(1234) +n = 24 +x = rand(rng, ["vegetables", "candy"], n) +y = [xᵢ == "candy" && rand(rng) > 0.3 ? "yummy" : "yuck" for xᵢ in x] +z = [yᵢ == "yummy" && rand(rng) > 0.6 ? "grown-up" : "child" for yᵢ in y] +d = CodifyVariables(UniqueElements()) +est = JointProbabilities(SECMI(base = 2), d) + +independence(SECMITest(est; nshuffles = 19), x, z, y) +``` \ No newline at end of file diff --git a/docs/src/independence.md b/docs/src/independence.md index 901f7ca1..d9340540 100644 --- a/docs/src/independence.md +++ b/docs/src/independence.md @@ -50,3 +50,10 @@ JDDTestResult CorrTest CorrTestResult ``` + +## [`SECMITest`](@ref) + +```@docs +SECMITest +SECMITestResult +``` \ No newline at end of file diff --git a/src/core.jl b/src/core.jl index db6946c4..85f03896 100644 --- a/src/core.jl +++ b/src/core.jl @@ -32,45 +32,46 @@ if each measure was implemented "in isolation". Concrete subtypes are given as input to [`association`](@ref). Many of these types require an [`AssociationMeasureEstimator`](@ref) to compute. -| Type | [`AssociationMeasure`](@ref) | Pairwise | Conditional | -| -------------------------- | ------------------------------------------- | :------: | :---------: | -| Correlation | [`PearsonCorrelation`](@ref) | ✓ | ✖ | -| Correlation | [`PartialCorrelation`](@ref) | ✓ | ✓ | -| Correlation | [`DistanceCorrelation`](@ref) | ✓ | ✓ | -| Correlation | [`ChatterjeeCorrelation`](@ref) | ✓ | ✖ | -| Correlation | [`AzadkiaChatterjeeCoefficient`](@ref) | ✓ | ✓ | -| Closeness | [`SMeasure`](@ref) | ✓ | ✖ | -| Closeness | [`HMeasure`](@ref) | ✓ | ✖ | -| Closeness | [`MMeasure`](@ref) | ✓ | ✖ | -| Closeness (ranks) | [`LMeasure`](@ref) | ✓ | ✖ | -| Closeness | [`JointDistanceDistribution`](@ref) | ✓ | ✖ | -| Cross-mapping | [`PairwiseAsymmetricInference`](@ref) | ✓ | ✖ | -| Cross-mapping | [`ConvergentCrossMapping`](@ref) | ✓ | ✖ | -| Conditional recurrence | [`MCR`](@ref) | ✓ | ✖ | -| Conditional recurrence | [`RMCD`](@ref) | ✓ | ✓ | -| Shared information | [`MIShannon`](@ref) | ✓ | ✖ | -| Shared information | [`MIRenyiJizba`](@ref) | ✓ | ✖ | -| Shared information | [`MIRenyiSarbu`](@ref) | ✓ | ✖ | -| Shared information | [`MITsallisFuruichi`](@ref) | ✓ | ✖ | -| Shared information | [`PartialCorrelation`](@ref) | ✖ | ✓ | -| Shared information | [`CMIShannon`](@ref) | ✖ | ✓ | -| Shared information | [`CMIRenyiSarbu`](@ref) | ✖ | ✓ | -| Shared information | [`CMIRenyiJizba`](@ref) | ✖ | ✓ | -| Shared information | [`CMIRenyiPoczos`](@ref) | ✖ | ✓ | -| Shared information | [`CMITsallisPapapetrou`](@ref) | ✖ | ✓ | -| Information transfer | [`TEShannon`](@ref) | ✓ | ✓ | -| Information transfer | [`TERenyiJizba`](@ref) | ✓ | ✓ | -| Partial mutual information | [`PartialMutualInformation`](@ref) | ✖ | ✓ | -| Information measure | [`JointEntropyShannon`](@ref) | ✓ | ✖ | -| Information measure | [`JointEntropyRenyi`](@ref) | ✓ | ✖ | -| Information measure | [`JointEntropyTsallis`](@ref) | ✓ | ✖ | -| Information measure | [`ConditionalEntropyShannon`](@ref) | ✓ | ✖ | -| Information measure | [`ConditionalEntropyTsallisAbe`](@ref) | ✓ | ✖ | -| Information measure | [`ConditionalEntropyTsallisFuruichi`](@ref) | ✓ | ✖ | -| Divergence | [`HellingerDistance`](@ref) | ✓ | ✖ | -| Divergence | [`KLDivergence`](@ref) | ✓ | ✖ | -| Divergence | [`RenyiDivergence`](@ref) | ✓ | ✖ | -| Divergence | [`VariationDistance`](@ref) | ✓ | ✖ | +| Type | [`AssociationMeasure`](@ref) | Pairwise | Conditional | +| -------------------------- | ---------------------------------------------------- | :------: | :---------: | +| Correlation | [`PearsonCorrelation`](@ref) | ✓ | ✖ | +| Correlation | [`PartialCorrelation`](@ref) | ✓ | ✓ | +| Correlation | [`DistanceCorrelation`](@ref) | ✓ | ✓ | +| Correlation | [`ChatterjeeCorrelation`](@ref) | ✓ | ✖ | +| Correlation | [`AzadkiaChatterjeeCoefficient`](@ref) | ✓ | ✓ | +| Closeness | [`SMeasure`](@ref) | ✓ | ✖ | +| Closeness | [`HMeasure`](@ref) | ✓ | ✖ | +| Closeness | [`MMeasure`](@ref) | ✓ | ✖ | +| Closeness (ranks) | [`LMeasure`](@ref) | ✓ | ✖ | +| Closeness | [`JointDistanceDistribution`](@ref) | ✓ | ✖ | +| Cross-mapping | [`PairwiseAsymmetricInference`](@ref) | ✓ | ✖ | +| Cross-mapping | [`ConvergentCrossMapping`](@ref) | ✓ | ✖ | +| Conditional recurrence | [`MCR`](@ref) | ✓ | ✖ | +| Conditional recurrence | [`RMCD`](@ref) | ✓ | ✓ | +| Shared information | [`MIShannon`](@ref) | ✓ | ✖ | +| Shared information | [`MIRenyiJizba`](@ref) | ✓ | ✖ | +| Shared information | [`MIRenyiSarbu`](@ref) | ✓ | ✖ | +| Shared information | [`MITsallisFuruichi`](@ref) | ✓ | ✖ | +| Shared information | [`PartialCorrelation`](@ref) | ✖ | ✓ | +| Shared information | [`CMIShannon`](@ref) | ✖ | ✓ | +| Shared information | [`CMIRenyiSarbu`](@ref) | ✖ | ✓ | +| Shared information | [`CMIRenyiJizba`](@ref) | ✖ | ✓ | +| Shared information | [`CMIRenyiPoczos`](@ref) | ✖ | ✓ | +| Shared information | [`CMITsallisPapapetrou`](@ref) | ✖ | ✓ | +| Shared information | [`ShortExpansionConditionalMutualInformation`](@ref) | ✖ | ✓ | +| Information transfer | [`TEShannon`](@ref) | ✓ | ✓ | +| Information transfer | [`TERenyiJizba`](@ref) | ✓ | ✓ | +| Partial mutual information | [`PartialMutualInformation`](@ref) | ✖ | ✓ | +| Information measure | [`JointEntropyShannon`](@ref) | ✓ | ✖ | +| Information measure | [`JointEntropyRenyi`](@ref) | ✓ | ✖ | +| Information measure | [`JointEntropyTsallis`](@ref) | ✓ | ✖ | +| Information measure | [`ConditionalEntropyShannon`](@ref) | ✓ | ✖ | +| Information measure | [`ConditionalEntropyTsallisAbe`](@ref) | ✓ | ✖ | +| Information measure | [`ConditionalEntropyTsallisFuruichi`](@ref) | ✓ | ✖ | +| Divergence | [`HellingerDistance`](@ref) | ✓ | ✖ | +| Divergence | [`KLDivergence`](@ref) | ✓ | ✖ | +| Divergence | [`RenyiDivergence`](@ref) | ✓ | ✖ | +| Divergence | [`VariationDistance`](@ref) | ✓ | ✖ | """ abstract type AssociationMeasure end @@ -88,45 +89,46 @@ Concrete subtypes are given as input to [`association`](@ref). ## Concrete implementations -| AssociationMeasure | Estimators | -| :------------------------------------------ | :--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- | -| [`PearsonCorrelation`](@ref) | Not required | -| [`DistanceCorrelation`](@ref) | Not required | -| [`PartialCorrelation`](@ref) | Not required | -| [`ChatterjeeCorrelation`](@ref) | Not required | -| [`AzadkiaChatterjeeCoefficient`](@ref) | Not required | -| [`SMeasure`](@ref) | Not required | -| [`HMeasure`](@ref) | Not required | -| [`MMeasure`](@ref) | Not required | -| [`LMeasure`](@ref) | Not required | -| [`JointDistanceDistribution`](@ref) | Not required | -| [`PairwiseAsymmetricInference`](@ref) | [`RandomVectors`](@ref), [`RandomSegment`](@ref) | -| [`ConvergentCrossMapping`](@ref) | [`RandomVectors`](@ref), [`RandomSegment`](@ref) | -| [`MCR`](@ref) | Not required | -| [`RMCD`](@ref) | Not required | -| [`MIShannon`](@ref) | [`JointProbabilities`](@ref), [`EntropyDecomposition`](@ref), [`KraskovStögbauerGrassberger1`](@ref), [`KraskovStögbauerGrassberger2`](@ref), [`GaoOhViswanath`](@ref), [`GaoKannanOhViswanath`](@ref), [`GaussianMI`](@ref) | -| [`MIRenyiJizba`](@ref) | [`JointProbabilities`](@ref), [`EntropyDecomposition`](@ref) | -| [`MIRenyiSarbu`](@ref) | [`JointProbabilities`](@ref) | -| [`MITsallisFuruichi`](@ref) | [`JointProbabilities`](@ref), [`EntropyDecomposition`](@ref) | -| [`MITsallisMartin`](@ref) | [`JointProbabilities`](@ref), [`EntropyDecomposition`](@ref) | -| [`CMIShannon`](@ref) | [`JointProbabilities`](@ref), [`EntropyDecomposition`](@ref), [`MIDecomposition`](@ref), [`GaussianCMI`](@ref), [`FPVP`](@ref), [`MesnerShalizi`](@ref), [`Rahimzamani`](@ref) | -| [`CMIRenyiSarbu`](@ref) | [`JointProbabilities`](@ref) | -| [`CMIRenyiJizba`](@ref) | [`JointProbabilities`](@ref), [`EntropyDecomposition`](@ref) | -| [`CMIRenyiPoczos`](@ref) | [`PoczosSchneiderCMI`](@ref) | -| [`CMITsallisPapapetrou`](@ref) | [`JointProbabilities`](@ref) | -| [`TEShannon`](@ref) | [`JointProbabilities`](@ref), [`EntropyDecomposition`](@ref), [`Zhu1`](@ref), [`Lindner`](@ref) | -| [`TERenyiJizba`](@ref) | [`JointProbabilities`](@ref) | -| [`PartialMutualInformation`](@ref) | [`JointProbabilities`](@ref) | -| [`JointEntropyShannon`](@ref) | [`JointProbabilities`](@ref) | -| [`JointEntropyRenyi`](@ref) | [`JointProbabilities`](@ref) | -| [`JointEntropyTsallis`](@ref) | [`JointProbabilities`](@ref) | -| [`ConditionalEntropyShannon`](@ref) | [`JointProbabilities`](@ref) | -| [`ConditionalEntropyTsallisAbe`](@ref) | [`JointProbabilities`](@ref) | -| [`ConditionalEntropyTsallisFuruichi`](@ref) | [`JointProbabilities`](@ref) | -| [`HellingerDistance`](@ref) | [`JointProbabilities`](@ref) | -| [`KLDivergence`](@ref) | [`JointProbabilities`](@ref) | -| [`RenyiDivergence`](@ref) | [`JointProbabilities`](@ref) | -| [`VariationDistance`](@ref) | [`JointProbabilities`](@ref) | +| AssociationMeasure | Estimators | +| :--------------------------------------------------- | :--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- | +| [`PearsonCorrelation`](@ref) | Not required | +| [`DistanceCorrelation`](@ref) | Not required | +| [`PartialCorrelation`](@ref) | Not required | +| [`ChatterjeeCorrelation`](@ref) | Not required | +| [`AzadkiaChatterjeeCoefficient`](@ref) | Not required | +| [`SMeasure`](@ref) | Not required | +| [`HMeasure`](@ref) | Not required | +| [`MMeasure`](@ref) | Not required | +| [`LMeasure`](@ref) | Not required | +| [`JointDistanceDistribution`](@ref) | Not required | +| [`PairwiseAsymmetricInference`](@ref) | [`RandomVectors`](@ref), [`RandomSegment`](@ref) | +| [`ConvergentCrossMapping`](@ref) | [`RandomVectors`](@ref), [`RandomSegment`](@ref) | +| [`MCR`](@ref) | Not required | +| [`RMCD`](@ref) | Not required | +| [`MIShannon`](@ref) | [`JointProbabilities`](@ref), [`EntropyDecomposition`](@ref), [`KraskovStögbauerGrassberger1`](@ref), [`KraskovStögbauerGrassberger2`](@ref), [`GaoOhViswanath`](@ref), [`GaoKannanOhViswanath`](@ref), [`GaussianMI`](@ref) | +| [`MIRenyiJizba`](@ref) | [`JointProbabilities`](@ref), [`EntropyDecomposition`](@ref) | +| [`MIRenyiSarbu`](@ref) | [`JointProbabilities`](@ref) | +| [`MITsallisFuruichi`](@ref) | [`JointProbabilities`](@ref), [`EntropyDecomposition`](@ref) | +| [`MITsallisMartin`](@ref) | [`JointProbabilities`](@ref), [`EntropyDecomposition`](@ref) | +| [`CMIShannon`](@ref) | [`JointProbabilities`](@ref), [`EntropyDecomposition`](@ref), [`MIDecomposition`](@ref), [`GaussianCMI`](@ref), [`FPVP`](@ref), [`MesnerShalizi`](@ref), [`Rahimzamani`](@ref) | +| [`CMIRenyiSarbu`](@ref) | [`JointProbabilities`](@ref) | +| [`CMIRenyiJizba`](@ref) | [`JointProbabilities`](@ref), [`EntropyDecomposition`](@ref) | +| [`CMIRenyiPoczos`](@ref) | [`PoczosSchneiderCMI`](@ref) | +| [`CMITsallisPapapetrou`](@ref) | [`JointProbabilities`](@ref) | +| [`TEShannon`](@ref) | [`JointProbabilities`](@ref), [`EntropyDecomposition`](@ref), [`Zhu1`](@ref), [`Lindner`](@ref) | +| [`TERenyiJizba`](@ref) | [`JointProbabilities`](@ref) | +| [`PartialMutualInformation`](@ref) | [`JointProbabilities`](@ref) | +| [`ShortExpansionConditionalMutualInformation`](@ref) | [`JointProbabilities`](@ref) | +| [`JointEntropyShannon`](@ref) | [`JointProbabilities`](@ref) | +| [`JointEntropyRenyi`](@ref) | [`JointProbabilities`](@ref) | +| [`JointEntropyTsallis`](@ref) | [`JointProbabilities`](@ref) | +| [`ConditionalEntropyShannon`](@ref) | [`JointProbabilities`](@ref) | +| [`ConditionalEntropyTsallisAbe`](@ref) | [`JointProbabilities`](@ref) | +| [`ConditionalEntropyTsallisFuruichi`](@ref) | [`JointProbabilities`](@ref) | +| [`HellingerDistance`](@ref) | [`JointProbabilities`](@ref) | +| [`KLDivergence`](@ref) | [`JointProbabilities`](@ref) | +| [`RenyiDivergence`](@ref) | [`JointProbabilities`](@ref) | +| [`VariationDistance`](@ref) | [`JointProbabilities`](@ref) | """ abstract type AssociationMeasureEstimator end diff --git a/src/independence_tests/independence.jl b/src/independence_tests/independence.jl index 891bc025..0e0dae20 100644 --- a/src/independence_tests/independence.jl +++ b/src/independence_tests/independence.jl @@ -7,6 +7,14 @@ using Statistics: quantile IndependenceTest <: IndependenceTest The supertype for all independence tests. + +## Concrete implementations + +- [`SurrogateAssociationTest`](@ref) +- [`LocalPermutationTest`](@ref) +- [`JointDistanceDistributionTest`](@ref) +- [`CorrTest`](@ref) +- [`SECMITest`](@ref) """ abstract type IndependenceTest{M} end @@ -22,12 +30,7 @@ If `z` is given too, then `test` must provide a conditional association measure. Returns a test `summary`, whose type depends on `test`. -## Compatible independence tests - -- [`SurrogateAssociationTest`](@ref) -- [`LocalPermutationTest`](@ref) -- [`JointDistanceDistributionTest`](@ref) -- [`CorrTest`](@ref) +See [`IndependenceTest`](@ref) for a list of compatible tests. """ function independence(test::IndependenceTest, x...) L = length(x) @@ -85,3 +88,4 @@ end include("parametric/parametric.jl") include("surrogate/SurrogateAssociationTest.jl") include("local_permutation/LocalPermutationTest.jl") +include("secmi/secmi_test.jl") \ No newline at end of file diff --git a/src/independence_tests/secmi/secmi_test.jl b/src/independence_tests/secmi/secmi_test.jl new file mode 100644 index 00000000..2e4af9c2 --- /dev/null +++ b/src/independence_tests/secmi/secmi_test.jl @@ -0,0 +1,124 @@ +export SECMITest +export SECMITestResult + +""" + SECMITest <: IndependenceTest + SECMITest(est; nshuffles = 19, surrogate = RandomShuffle(), rng = Random.default_rng()) + +A test for conditional independence based on the [`ShortExpansionConditionalMutualInformation`](@ref) measure [Kubkowski2021](@cite). + +The first argument `est` must be a [`InformationMeasureEstimator`](@ref) that provides the +[`ShortExpansionConditionalMutualInformation`](@ref) instance. See examples below. + +## Examples + +- [Example 1](@ref example_SECMITEST_JointProbabilities_CodifyVariables_ValueBinning): Independence test for small sample sizes using [`CodifyVariables`](@ref) with + [`ValueBinning`](@ref) discretization. +- [Example 2](@ref example_SECMITEST_JointProbabilities_CodifyVariables_UniqueElements): Independence test for small sample + sizes with categorical data (using [`CodifyVariables`](@ref) with [`UniqueElements`](@ref) discretization). +""" +struct SECMITest{E, S, I, RNG} <: IndependenceTest{E} + # really, this can only be an estimator, but we name it `est_or_measure` for consistency + # with the remaining tests, so we don't have to write custom code. + est_or_measure::E + surrogate::S + nshuffles::I + rng::RNG +end + +function SECMITest(est; nshuffles = 19, surrogate = RandomShuffle(), rng = Random.default_rng()) + return SECMITest(est, surrogate, nshuffles, rng) +end + +""" + SECMITestResult <: IndependenceTestResult + SECMITestResult(secmi₀, secmiₖ, p, μ̂, σ̂, emp_cdf, D𝒩, D𝒳², nshuffles::Int) + +A simple struct that holds the computed parameters of a [`SECMITest`](@ref) when called +with [`independence`](@ref), as described in [Kubkowski2021](@cite). + +## Parameters + +- `p`: The p-value for the test. +- `secmi₀`: The value of the [`ShortExpansionConditionalMutualInformation`](@ref) measure estimated on the original data. +- `secmiₖ`: An ensemble of values for the [`ShortExpansionConditionalMutualInformation`](@ref) measure estimated on triples + `SECMI(X̂, Y, Z)`, where `X̂` indicates a shuffled version of the first variable `X` + and `length(secmiₖ) == nshuffles`. +- `μ̂`: The estimated mean of the `secmiₖ`. +- `σ̂`: The estimated standard deviation of the `secmiₖ`. +- `emp_cdf`: The empirical cumulative distribution function (CDF) of the `secmiₖ`s. +- `D𝒩`: The ``D_{N(\\hat{\\mu}, \\hat{\\sigma})}`` statistic. +- `D𝒳²`: The ``D_{\\chi^2}`` statistic. +""" +struct SECMITestResult{S0, SK, P, MU, S, E, DN, DCHI} <: IndependenceTestResult + n_vars::Int # 3 vars = conditional (always 3) + secmi₀::S0 # the value of the measure, non-permuted + secmiₖ::SK # the values of the measure, permuted `nshuffles` times + pvalue::P # the p-value, computed from the estimatedd parameters below. + μ̂::MU + σ̂::S + emp_cdf::E + D𝒩::DN + D𝒳²::DCHI +end +pvalue(r::SECMITestResult) = r.pvalue + +function Base.show(io::IO, test::SECMITestResult) + print(io, + """\ + `SECMITEST` independence test + $(null_hypothesis_text(test::IndependenceTestResult)) + Estimated parameters: + μ̂ = $(test.μ̂), σ̂ = $(test.σ̂) + D𝒩 = $(test.D𝒩), D𝒳² = $(test.D𝒳²) + + $(pvalue_text_summary(test)) + """ + ) +end + +function independence(test::SECMITest, x, y, z) + (; est_or_measure, surrogate, nshuffles, rng) = test + secmi₀ = association(est_or_measure, x, y, z) + sx = surrogenerator(x, surrogate, rng) + secmiₖ = zeros(test.nshuffles) + for k = 1:nshuffles + secmiₖ[k] = association(est_or_measure, sx(), y, z) + end + μ̂ = 1/nshuffles * sum(secmiₖ) + σ̂ = 1/(nshuffles - 1) * sum((sₖ - μ̂)^2 for sₖ in secmiₖ) + emp_cdf = ecdf(secmiₖ) + F𝒩 = Normal(μ̂, σ̂) + + if μ̂ ≤ 0.0 + p = 1 - cdf(F𝒩, secmi₀) + return SECMITestResult(3, secmi₀, secmiₖ, p, μ̂, σ̂, emp_cdf, nothing, nothing) + else + # Degrees of freedom for Chi squared distribution estimated as the mean of the `secmiₖ` + # (page 18 in Kubkowski et al.). The `Chisq` distribution is only defined for μ̂ > 0, + # so we put μ̂ <= 0.0 in a separate criterion first to avoid errors. + F𝒳² = Chisq(μ̂) + D𝒩, D𝒳² = sup_values(emp_cdf, F𝒩, F𝒳², secmiₖ) + if D𝒩 < D𝒳² + p = 1 - cdf(F𝒩, secmi₀) + else + p = 1 - cdf(F𝒳², secmi₀) + end + return SECMITestResult(3, secmi₀, secmiₖ, p, μ̂, σ̂, emp_cdf, D𝒩, D𝒳²) + + end + + + return SECMITestResult(3, secmi₀, secmiₖ, p, μ̂, σ̂, emp_cdf, D𝒩, D𝒳²) +end + + + +function sup_values(emp_cdf, F𝒩, F𝒳², secmiₖ) + empirical_cdf_values = emp_cdf.(secmiₖ) + normal_cdf_values = cdf.(F𝒩, secmiₖ) + chi2_cdf_values = cdf.(F𝒳², secmiₖ) + D𝒩 = maximum(abs.(empirical_cdf_values .- normal_cdf_values)) + D𝒳² = maximum(abs.(empirical_cdf_values .- chi2_cdf_values)) + return D𝒩, D𝒳² +end diff --git a/src/methods/information/counts_and_probs/counts.jl b/src/methods/information/counts_and_probs/counts.jl index 6c84d944..abd4ecc9 100644 --- a/src/methods/information/counts_and_probs/counts.jl +++ b/src/methods/information/counts_and_probs/counts.jl @@ -204,6 +204,11 @@ function counts(discretization::CodifyVariables{1}, x::Vararg{ArrayOrStateSpaceS return counts(x̂...) end +function counts(d::CodifyVariables{1, UniqueElements}, x::Vararg{ArrayOrStateSpaceSet, N}) where N + o = first(d.outcome_spaces) + return counts(o, x...) +end + as_vec(x::AbstractStateSpaceSet{1}) = [first(xᵢ) for xᵢ in vec(x)] diff --git a/src/methods/information/counts_and_probs/probabilities.jl b/src/methods/information/counts_and_probs/probabilities.jl index 18cf6537..986b4940 100644 --- a/src/methods/information/counts_and_probs/probabilities.jl +++ b/src/methods/information/counts_and_probs/probabilities.jl @@ -125,4 +125,4 @@ end function probabilities(discretization::CodifyVariables, x::Vararg{ArrayOrStateSpaceSet, N}) where N cts = counts(discretization, x...) return probabilities(RelativeAmount(), cts) -end +end \ No newline at end of file diff --git a/src/methods/information/definitions/conditional_mutual_informations/conditional_mutual_informations.jl b/src/methods/information/definitions/conditional_mutual_informations/conditional_mutual_informations.jl index 7b6e7f5f..74560b91 100644 --- a/src/methods/information/definitions/conditional_mutual_informations/conditional_mutual_informations.jl +++ b/src/methods/information/definitions/conditional_mutual_informations/conditional_mutual_informations.jl @@ -41,10 +41,10 @@ function marginal_entropies_cmi4h_differential(est::EntropyDecomposition{<:Condi XYZ = StateSpaceSet(X, Y, Z) modified_est = estimator_with_overridden_parameters(est.definition, est.est) - HXZ = entropy(modified_est, XZ) - HYZ = entropy(modified_est, YZ) - HXYZ = entropy(modified_est, XYZ) - HZ = entropy(modified_est, Z) + HXZ = information(modified_est, XZ) + HYZ = information(modified_est, YZ) + HXYZ = information(modified_est, XYZ) + HZ = information(modified_est, Z) return HXZ, HYZ, HXYZ, HZ end diff --git a/src/methods/information/definitions/information_definitions.jl b/src/methods/information/definitions/information_definitions.jl index bfc84f86..6659e8ff 100644 --- a/src/methods/information/definitions/information_definitions.jl +++ b/src/methods/information/definitions/information_definitions.jl @@ -6,5 +6,6 @@ include("mutual_informations/mutual_informations.jl") include("conditional_mutual_informations/conditional_mutual_informations.jl") include("partial_mutual_information/partial_mutual_information.jl") include("transferentropy/transfer_entropies.jl") +include("short_expansion_conditional_mutual_information/short_expansion_conditional_mutual_information.jl") include("override_parameters.jl") \ No newline at end of file diff --git a/src/methods/information/definitions/short_expansion_conditional_mutual_information/short_expansion_conditional_mutual_information.jl b/src/methods/information/definitions/short_expansion_conditional_mutual_information/short_expansion_conditional_mutual_information.jl new file mode 100644 index 00000000..b8f6ad85 --- /dev/null +++ b/src/methods/information/definitions/short_expansion_conditional_mutual_information/short_expansion_conditional_mutual_information.jl @@ -0,0 +1,84 @@ +using TimeseriesSurrogates +using Random +using StatsBase: ecdf +using Distributions + +export ShortExpansionConditionalMutualInformation, SECMI + +""" + ShortExpansionConditionalMutualInformation <: MultivariateInformationMeasure + ShortExpansionConditionalMutualInformation(; base = 2) + SECMI(; base = 2) # alias + +The short expansion of (Shannon) conditional mutual information (SECMI) measure +from [Kubkowski2021](@citet). + +## Description + +The SECMI measure is defined as + +```math +SECMI(X,Y|Z) = I(X,Y) + \\sum_{k=1}^{m} II(X,Z_k,Y) = (1 - m) I(X,Y) + \\sum_{k=1}^{m} I(X,Y|Z_k). +``` + +This quantity is estimated from data using one of the estimators below from the formula + +```math +\\widehat{SECMI}(X,Y|Z) = \\widehat{I}(X,Y) + \\sum_{k=1}^{m} \\widehat{II}(X,Z_k,Y) = (1 - m) \\widehat{I}(X,Y) + \\sum_{k=1}^{m} \\widehat{I}(X,Y|Z_k). +``` + +## Compatible estimators + +- [`JointProbabilities`](@ref). + +## Estimation + +- [Example 1](@ref example_ShortExpansionConditionalMutualInformation_JointProbabilities_CodifyVariables_ValueBinning): + Estimating [`ShortExpansionConditionalMutualInformation`](@ref) using the [`JointProbabilities`](@ref) estimator using a + [`CodifyVariables`](@ref) with [`ValueBinning`](@ref) discretization. +""" +Base.@kwdef struct ShortExpansionConditionalMutualInformation{B} <: MultivariateInformationMeasure + base::B = 2 +end +const SECMI = ShortExpansionConditionalMutualInformation + +function Base.show(io::IO, m::ShortExpansionConditionalMutualInformation) + msg = "ShortExpansionConditionalMutualInformation(; base = $(m.base))" + print(io, msg) +end +min_inputs_vars(::ShortExpansionConditionalMutualInformation) = 3 +max_inputs_vars(::ShortExpansionConditionalMutualInformation) = Inf + +# Assumes 1st dimension of `probs` corresponds to X, 2nd dimension of `probs` +# corresponds to Y, and dimensions `3:ndims(probs)` correspond to marginals Zₖ, +function association(definition::SECMI, probs::Probabilities{T, N}) where {T, N} + @assert N >= 3 + (; base) = definition + + def_mi = MIShannon(; base = base) + def_cmi = CMIShannon(; base = base) + + m = ndims(probs) - 2 + pXY = marginal(probs, dims = 1:2) + mi_XY = association(def_mi, pXY) + cmis = 0.0 + for k = 1:m + # association(definition::CMIShannon, pxyz::Probabilities{T, 3}) + # is the call signature, so we simply replace the last variable + # with the marginal Zₖ for each Î(X, Y | Zₖ) in the sum + cmis += association(def_cmi, marginal(probs, dims = (1, 2, 2 + k))) + end + return (1 - m) * mi_XY + cmis +end + +function association(est::JointProbabilities{<:SECMI}, x, y, z...) + probs = probabilities(est.discretization, x, y, z...) + return association(est.definition, probs) +end + +# SECMI operates column-wise on the conditional variable, so if a statespace set +# is given, then we operate on columns +function association(est::JointProbabilities{<:SECMI}, x, y, z::AbstractStateSpaceSet) + probs = probabilities(est.discretization, x, y, columns(z)...) + return association(est.definition, probs) +end diff --git a/test/causal_graphs/oce.jl b/test/causal_graphs/oce.jl index c8a13983..6854dfa2 100644 --- a/test/causal_graphs/oce.jl +++ b/test/causal_graphs/oce.jl @@ -29,6 +29,7 @@ using DynamicalSystemsBase GaussianCMI(), AzadkiaChatterjeeCoefficient(), PartialCorrelation(), + JointProbabilities(SECMI(), CodifyVariables(OrdinalPatterns())) ] for mi_est in mi_ests utest = SurrogateAssociationTest(mi_est; rng, nshuffles = 2) @@ -37,9 +38,21 @@ using DynamicalSystemsBase alg = OCE(; utest, ctest, τmax = 1) parents = infer_graph(alg, X; verbose = false) @test parents isa Vector{<:OCESelectedParents} + + if (cmi_est isa JointProbabilities{<:SECMI}) + ctest = SECMITest(cmi_est; nshuffles = 2) + alg = OCE(; utest, ctest, τmax = 1) + parents = infer_graph(alg, X; verbose = false) + @test parents isa Vector{<:OCESelectedParents} + end end end end +# ---------------------------------------------------------------- +# Some more that don't fint easily inside the loop. +# ---------------------------------------------------------------- + + # ---------------------------------------------------------------- # A few examples with more data and more iterations # ---------------------------------------------------------------- diff --git a/test/independence/SECMITest.jl b/test/independence/SECMITest.jl new file mode 100644 index 00000000..9874b32a --- /dev/null +++ b/test/independence/SECMITest.jl @@ -0,0 +1,31 @@ +n = 40 +using Associations +using Test +using Random; rng = Xoshiro(1234) +x = rand(rng, n) +y = randn(rng, n) .+ x .^ 2 +z = randn(rng, n) .* y + +# An estimator for estimating the SECMI measure +est = JointProbabilities(SECMI(base = 2), CodifyVariables(ValueBinning(3))) +test = SECMITest(est; nshuffles = 19, rng = rng) + +# Do a test and check that we can reject null or not as expected +α = 0.05 +@test pvalue(independence(test, x, y, z)) < α +@test pvalue(independence(test, x, z, y)) > α + +out = repr(independence(SECMITest(est; nshuffles = 2), x, y, z)) +@test occursin("D𝒳²", out) + +# Categorical +rng = Xoshiro(1234) +n = 12 +x = rand(rng, ["vegetables", "candy"], n) +y = [xᵢ == "candy" && rand(rng) > 0.3 ? "yummy" : "yuck" for xᵢ in x] +z = [yᵢ == "yummy" && rand(rng) > 0.6 ? "grown-up" : "child" for yᵢ in y] +d = CodifyVariables(UniqueElements()) +est = JointProbabilities(SECMI(base = 2), d) +test = SECMITest(est; nshuffles = 19, rng = rng) +@test pvalue(independence(test, x, y, z)) < α +@test pvalue(independence(test, x, z, y)) > α \ No newline at end of file diff --git a/test/independence/SurrogateAssociationTest/chatterjee_correlation.jl b/test/independence/SurrogateAssociationTest/chatterjee_correlation.jl index 23b61580..1fac5481 100644 --- a/test/independence/SurrogateAssociationTest/chatterjee_correlation.jl +++ b/test/independence/SurrogateAssociationTest/chatterjee_correlation.jl @@ -7,12 +7,12 @@ rng = Xoshiro(1234) test = SurrogateAssociationTest(ChatterjeeCorrelation(), nshuffles = 19, rng = rng) # We expect that we *cannot* reject the null hypothesis for independent variables -x = rand(rng, 1:10, 100) -y = rand(rng, 1:10, 100) +x = rand(rng, 1:10, 200) +y = rand(rng, 1:10, 200) α = 0.05 @test pvalue(independence(test, x, y)) > α # We expect that we *can* reject the null hypothesis for for dependent variables. -w = rand(rng, 1:10, 100) -z = rand(rng, 1:10, 100) .* sin.(w) .+ cos.(w) +w = rand(rng, 1:10, 200) +z = rand(rng, 1:10, 200) .* sin.(w) .+ cos.(w) @test pvalue(independence(test, z, w)) < α diff --git a/test/independence/independence.jl b/test/independence/independence.jl index 2bafe60a..623ed9c4 100644 --- a/test/independence/independence.jl +++ b/test/independence/independence.jl @@ -1,5 +1,5 @@ include("LocalPermutationTest/local_permutation_test.jl") include("SurrogateAssociationTest/SurrogateAssociationTest.jl") include("JointDistanceDistributionTest.jl") -#include("PATest.jl") include("CorrTest.jl") +include("SECMITest.jl") \ No newline at end of file diff --git a/test/methods/information/information.jl b/test/methods/information/information.jl index 94670347..d31b80c5 100644 --- a/test/methods/information/information.jl +++ b/test/methods/information/information.jl @@ -6,6 +6,7 @@ include("conditional_entropies/conditional_entropies.jl") include("mutual_informations/mutual_informations.jl") include("conditional_mutual_informations/conditional_mutual_informations.jl") include("transfer_entropies/transfer_entropies.jl") +include("secmi/secmi.jl") # Estimators of the information measures. include("estimators/estimators.jl") diff --git a/test/methods/information/secmi/secmi.jl b/test/methods/information/secmi/secmi.jl new file mode 100644 index 00000000..c660dc1f --- /dev/null +++ b/test/methods/information/secmi/secmi.jl @@ -0,0 +1,21 @@ +using Associations +using Test +using Random; rng = Xoshiro(8284) + +# Categorical +n = 30 +x = rand(rng, ["vegetables", "candy"], n) +y = [xᵢ == "candy" && rand(rng) > 0.3 ? "yummy" : "yuck" for xᵢ in x] +z = [yᵢ == "yummy" && rand(rng) > 0.6 ? "grown-up" : "child" for yᵢ in y] +w = [yᵢ == "yummy" && rand(rng) > 0.4 ? "tv" : "nature" for yᵢ in y] + +d = CodifyVariables(UniqueElements()) +est = JointProbabilities(SECMI(base = 2), d) + +# Test both signatures +@test association(est, x, y, z) ≥ 0.0 +@test association(est, x, y, StateSpaceSet(z, w)) ≥ 0.0 + +# When the conditioning set has cardinality 1, we should get the regular CMI +est_cmi = JointProbabilities(CMIShannon(base = 2), d) +@test association(est_cmi, x, y, z) ≈ association(est, x, y, z)