From c10ed326cde70b5196fea16296e620fabd15a4db Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20Agas=C3=B8ster=20Haaga?= Date: Tue, 30 Nov 2021 16:46:13 +0100 Subject: [PATCH 01/29] Draft of Leanings module --- src/CausalityTools.jl | 1 + src/Leanings/Leanings.jl | 5 +++ src/Leanings/penchants.jl | 92 +++++++++++++++++++++++++++++++++++++++ 3 files changed, 98 insertions(+) create mode 100644 src/Leanings/Leanings.jl create mode 100644 src/Leanings/penchants.jl diff --git a/src/CausalityTools.jl b/src/CausalityTools.jl index 9a5c11d4d..9dbbb8f36 100644 --- a/src/CausalityTools.jl +++ b/src/CausalityTools.jl @@ -15,6 +15,7 @@ module CausalityTools include("CrossMappings/CrossMappings.jl") include("SMeasure/smeasure.jl") include("PredictiveAsymmetry/PredictiveAsymmetry.jl") + include("Leanings/Leanings.jl") include("example_systems/ExampleSystems.jl") using Requires diff --git a/src/Leanings/Leanings.jl b/src/Leanings/Leanings.jl new file mode 100644 index 000000000..20858c23c --- /dev/null +++ b/src/Leanings/Leanings.jl @@ -0,0 +1,5 @@ +using Reexport + +@reexport module Leanings + include("penchants.jl") +end \ No newline at end of file diff --git a/src/Leanings/penchants.jl b/src/Leanings/penchants.jl new file mode 100644 index 000000000..684ed095c --- /dev/null +++ b/src/Leanings/penchants.jl @@ -0,0 +1,92 @@ + +x = [0, 0, 1, 0, 0, 1, 0, 0, 1, 0] +y = [0, 0, 0, 1, 0, 0, 1, 0, 0, 1] +export lean + +""" + lean_ex(C, E, l) + +Computes the leaning of the `l`-assignment ``\\{C, E\\} = \\{X, Y\\}, \\{X_{t-l}, Y_t\\}``. +""" +function lean(X, Y, l = 1) + @assert length(X) == length(Y) + L = length(X) + + # Leave reasonably many for probability computations. `l` can't be too large. + @assert l < L ÷ 2 + + #P(y_t = 1 | x_{t-1} = 1) = n_ec/n_c + + # Number of time the state y_t = 1 AND x_{t-1} = 1 appears in {X, Y} + n_ecs = Vector{Int}(undef, 0) + n_cs = Vector{Int}(undef, 0) + n_es = Vector{Int}(undef, 0) + + states = unique([X; Y]) + + # Define an iterator over all penchants + penchants = Iterators.product(states, states) + @show states + @show penchants |> collect + + + + + for penchant in penchants + @show penchant + sx, sy = penchant[1], penchant[2] + + # Number of times the assumed cause has appeared + n_c = 0 + + # Number of times the assumed effect has appeared + n_e = 0 + + # Number of time the state y_t = a AND x_{t-1} = b appears in {X, Y}, + # where (a, b) is one of the length(states)^2 possible states. + n_ec = 0 + for t = l+1:L + effect = Y[t] + cause = X[t-l] + + if effect == sx && cause == sy + n_ec += 1 + end + + + if effect == sx + n_e += 1 + end + + if cause == sy + n_c += 1 + end + end + push!(n_ecs, n_ec) + push!(n_cs, n_c) + push!(n_es, n_e) + end + + @show Ps_cond = n_ecs ./ n_cs + @show Ps_e = n_es ./ (L - l) + @show Ps_c = n_cs ./ (L - l) + + # for t = l+1:L + # yt = Y[t] + # xtl = X[t-l] + # @show Y[t], X[t-l] + # if yt == 1 && xtl == 1 + # n_ec += 1 + # end + + # if X[t-l] == 1 + # n_c += 1 + # end + # end + + # @show n_ec + # @show n_c + + # p = n_ec/n_c +end + From 1711df73939495b628f5434bf986fd6f8042d9c6 Mon Sep 17 00:00:00 2001 From: Kristian Haaga Date: Wed, 1 Dec 2021 05:48:01 +0100 Subject: [PATCH 02/29] Leanings, with tests and docs --- docs/make.jl | 3 + docs/src/leanings.md | 9 ++ src/Leanings/penchants.jl | 166 +++++++++++++++++++++++----------- test/methods/test_leanings.jl | 9 ++ test/runtests.jl | 1 + 5 files changed, 135 insertions(+), 53 deletions(-) create mode 100644 docs/src/leanings.md create mode 100644 test/methods/test_leanings.jl diff --git a/docs/make.jl b/docs/make.jl index 1c618fa45..a28af6007 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -48,6 +48,9 @@ PAGES = [ "generalized_entropy.md", "info_estimators.md", ], + "Exploratory" => [ + "leanings.md" + ], "example_systems.md", "Utilities" => [ diff --git a/docs/src/leanings.md b/docs/src/leanings.md new file mode 100644 index 000000000..99587e144 --- /dev/null +++ b/docs/src/leanings.md @@ -0,0 +1,9 @@ +# Leanings + +Leanings are probability-based exploratory causal inference measures that rely on the estimation of probabilities directly from the data. Here, we provide +implementations of the *penchants and leanings* algoriths from McCracken and Weigel (2016). + +```@docs +mean_observed_penchant +mean_leaning +``` \ No newline at end of file diff --git a/src/Leanings/penchants.jl b/src/Leanings/penchants.jl index 684ed095c..0cd6db50b 100644 --- a/src/Leanings/penchants.jl +++ b/src/Leanings/penchants.jl @@ -1,39 +1,41 @@ x = [0, 0, 1, 0, 0, 1, 0, 0, 1, 0] y = [0, 0, 0, 1, 0, 0, 1, 0, 0, 1] -export lean +export mean_observed_penchant, mean_leaning -""" - lean_ex(C, E, l) +using Statistics -Computes the leaning of the `l`-assignment ``\\{C, E\\} = \\{X, Y\\}, \\{X_{t-l}, Y_t\\}``. -""" -function lean(X, Y, l = 1) - @assert length(X) == length(Y) - L = length(X) - - # Leave reasonably many for probability computations. `l` can't be too large. - @assert l < L ÷ 2 +# We just assume that the time series have been pre-binned or symbolized +# However, `get_states`` might turn out much more complicated if we decide +# to implement more complicated cause-effect assignments +get_states(x, y) = unique([x; y]) - #P(y_t = 1 | x_{t-1} = 1) = n_ec/n_c +# Getting the penchants might also become more complex for other +# complicated cause-effect assignments, so separate it out in a function. +function get_penchants(x, y) + states = get_states(x, y) - # Number of time the state y_t = 1 AND x_{t-1} = 1 appears in {X, Y} - n_ecs = Vector{Int}(undef, 0) - n_cs = Vector{Int}(undef, 0) - n_es = Vector{Int}(undef, 0) + # All possible 2-permutations (with repetitions) of the unique + # states in `x ∪ y`. + return Iterators.product(states, states) +end - states = unique([X; Y]) - - # Define an iterator over all penchants - penchants = Iterators.product(states, states) - @show states - @show penchants |> collect +function penchant_counts(x, y, l = 1) + @assert length(x) == length(y) + n = length(x) + # Leave reasonably many for probability computations. `l` can't be too large. + @assert l < n ÷ 2 + + # Number of time the state y_t = 1 AND x_{t-1} = 1 appears in {x, y} + ns_ec = Vector{Int}(undef, 0) + ns_c = Vector{Int}(undef, 0) + ns_e = Vector{Int}(undef, 0) - + # For each possible penchant, count occurrences. + penchants = get_penchants(x, y) - for penchant in penchants - @show penchant + @inbounds for penchant in penchants sx, sy = penchant[1], penchant[2] # Number of times the assumed cause has appeared @@ -42,12 +44,14 @@ function lean(X, Y, l = 1) # Number of times the assumed effect has appeared n_e = 0 - # Number of time the state y_t = a AND x_{t-1} = b appears in {X, Y}, - # where (a, b) is one of the length(states)^2 possible states. + # Number of times the cause AND effect has appeared simulaneously n_ec = 0 - for t = l+1:L - effect = Y[t] - cause = X[t-l] + + for t = l+1:n + # "we required that the cause must precede the effect". This is the + # standard l-assignment + effect = y[t] + cause = x[t-l] if effect == sx && cause == sy n_ec += 1 @@ -62,31 +66,87 @@ function lean(X, Y, l = 1) n_c += 1 end end - push!(n_ecs, n_ec) - push!(n_cs, n_c) - push!(n_es, n_e) + push!(ns_ec, n_ec) + push!(ns_c, n_c) + push!(ns_e, n_e) + end + + return ns_c, ns_e, ns_ec +end + + +function mean_observed_ρ(ρs, κs, ns_ec, L, weighted = false) + ρ = 0.0 + ct = 0 + @inbounds for i = 1:length(κs) + if κs[i] != 0.0 + if weighted + ρ += ρs[i] * ns_ec[i] + else + ρ += ρs[i] + end + ct += 1 + end end - @show Ps_cond = n_ecs ./ n_cs - @show Ps_e = n_es ./ (L - l) - @show Ps_c = n_cs ./ (L - l) + if weighted + ρ / L + else + return ρ / ct + end +end + +""" + mean_observed_penchant(x, y, l, weighted = false) → ρ̄ + +Computes the *mean observed penchant* of the `l`-assignment ``\\{C, E\\} = \\{x, y\\}, \\{x_{t-l}, y_t\\}``. + +If `weighted == true`, then compute the *weighted mean observed penchant*. + +[^McCrackenWeigel2016]: McCracken, J. M., & Weigel, R. S. (2016). Nonparametric causal inference for bivariate time series. Physical Review E, 93(2), 022207. +""" +function mean_observed_penchant(x, y, l = 1; weighted = false) + n = length(x) + ns_c, ns_e, ns_ec = penchant_counts(x, y, l) + + L = n - l # "library length", resulting from shifting the time series l steps + Ps_e = ns_e ./ L + Ps_c = ns_c ./ L + κs = ns_ec ./ ns_c + ρs = κs.* (1 .+ Ps_c ./ (1 .- Ps_c)) .- (Ps_e) ./ (1 .- Ps_c) - # for t = l+1:L - # yt = Y[t] - # xtl = X[t-l] - # @show Y[t], X[t-l] - # if yt == 1 && xtl == 1 - # n_ec += 1 - # end - - # if X[t-l] == 1 - # n_c += 1 - # end - # end - - # @show n_ec - # @show n_c - - # p = n_ec/n_c + # Mean observed penchant. We can determine which penchants are observed + # by looking for zeros in `κs`. If `κs[i] == 0`, then the `i`-th penchant + # is unobserved. Why? When the penchant is unobserved, there are no + # instances of the cause and effect occurring together, so + # `P(effect | cause) = P(effect, cause) / P(cause) = 0 / P(cause) = 0` + # As explained in section D in the paper, including unobserved penchants + # would be meaningless for causal inference. + mean_observed_ρ(ρs, κs, ns_ec, L, weighted) end +""" + + mean_leaning(x, y, l = 1, weighted = true) → ρ̄ + +Compute the *mean observed leaning* `ρ̄` (McCracken & Weigel, 2016)[^McCrackenWeigel2016] between scalar time series `x` and `y`, +based on the standard `l`-assignment of the `l`-assignment ``\\{C, E\\} = \\{x, y\\}, \\{x_{t-l}, y_t\\}``. + +If `ρ̄ > 0`, then the probability that `x` drives `y` is higher than the probability +that `y` drives `x`. Vice versa, if `ρ̄ > 0`, then the probability that `y` drives `x` is higher than the probability +that `x` drives `y`. + +If `weighted == true`, then each penchant is weighted by the number of times it appears in +the data, which yields the *weighted mean observed leaning*. This the default behavior, since +"...[the weighted mean observed leaning] accounts for the frequency of observed cause-effect pairs +within the data, which is assumed to be a predictor of causal influence" (McCracken & Weigel, 2016). + +This implementation *does not* take care of the tolerance domains discussed in the paper. +Pre-process your time series, using appropriate binning or symbolization schemes, +so that `x` and `y` both consists of discrete symbols. + +[^McCrackenWeigel2016]: McCracken, J. M., & Weigel, R. S. (2016). Nonparametric causal inference for bivariate time series. Physical Review E, 93(2), 022207. +""" +function mean_leaning(x, y, l = 1, weighted = true) + lean(x, y, l, weighted = weighted) - lean(y, x, l, weighted = weighted) +end \ No newline at end of file diff --git a/test/methods/test_leanings.jl b/test/methods/test_leanings.jl new file mode 100644 index 000000000..a2aa56c9d --- /dev/null +++ b/test/methods/test_leanings.jl @@ -0,0 +1,9 @@ + +x = [0, 0, 1, 0, 0, 1, 0, 0, 1, 0] +y = [0, 0, 0, 1, 0, 0, 1, 0, 0, 1] + +@test mean_observed_penchant(x, y, 1, weighted = false) ≈ 1.0 +@test mean_observed_penchant(y, x, 1, weighted = false) ≈ 1/7 + +@test mean_observed_penchant(x, y, 1, weighted = true) ≈ 1.0 +@test mean_observed_penchant(y, x, 1, weighted = true) ≈ 3/63 \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 58099aebe..1c3c61e7a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -8,6 +8,7 @@ include("methods/test_smeasure.jl") include("methods/test_joint_distance_distribution.jl") include("methods/test_predictive_asymmetry.jl") include("methods/test_crossmapping.jl") +include("methods/test_leanings.jl") include("systems/discrete/test_discrete_systems.jl") include("systems/continuous/test_continuous_systems.jl") From 677a6760b6629a20ff452fd0e348f88de3991228 Mon Sep 17 00:00:00 2001 From: Kristian Haaga Date: Wed, 1 Dec 2021 05:49:25 +0100 Subject: [PATCH 03/29] Correct function name --- src/Leanings/penchants.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/Leanings/penchants.jl b/src/Leanings/penchants.jl index 0cd6db50b..3a5186256 100644 --- a/src/Leanings/penchants.jl +++ b/src/Leanings/penchants.jl @@ -148,5 +148,6 @@ so that `x` and `y` both consists of discrete symbols. [^McCrackenWeigel2016]: McCracken, J. M., & Weigel, R. S. (2016). Nonparametric causal inference for bivariate time series. Physical Review E, 93(2), 022207. """ function mean_leaning(x, y, l = 1, weighted = true) - lean(x, y, l, weighted = weighted) - lean(y, x, l, weighted = weighted) + return mean_observed_penchant(x, y, l, weighted = weighted) - + mean_observed_penchant(y, x, l, weighted = weighted) end \ No newline at end of file From 208e4a5edc74403b1aa98030a1505b1d173d1ab2 Mon Sep 17 00:00:00 2001 From: Kristian Haaga Date: Wed, 1 Dec 2021 05:51:08 +0100 Subject: [PATCH 04/29] Use testset --- test/methods/test_leanings.jl | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/test/methods/test_leanings.jl b/test/methods/test_leanings.jl index a2aa56c9d..1f5fea0e1 100644 --- a/test/methods/test_leanings.jl +++ b/test/methods/test_leanings.jl @@ -2,8 +2,10 @@ x = [0, 0, 1, 0, 0, 1, 0, 0, 1, 0] y = [0, 0, 0, 1, 0, 0, 1, 0, 0, 1] -@test mean_observed_penchant(x, y, 1, weighted = false) ≈ 1.0 -@test mean_observed_penchant(y, x, 1, weighted = false) ≈ 1/7 +@testset "Leanings" begin + @test mean_observed_penchant(x, y, 1, weighted = false) ≈ 1.0 + @test mean_observed_penchant(y, x, 1, weighted = false) ≈ 1/7 -@test mean_observed_penchant(x, y, 1, weighted = true) ≈ 1.0 -@test mean_observed_penchant(y, x, 1, weighted = true) ≈ 3/63 \ No newline at end of file + @test mean_observed_penchant(x, y, 1, weighted = true) ≈ 1.0 + @test mean_observed_penchant(y, x, 1, weighted = true) ≈ 3/63 +end \ No newline at end of file From a1c620bdf80e7cd3ed79e5a9028c3e732ce5ff4b Mon Sep 17 00:00:00 2001 From: Kristian Haaga Date: Wed, 1 Dec 2021 06:14:05 +0100 Subject: [PATCH 05/29] Remove example arrays --- src/Leanings/penchants.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/Leanings/penchants.jl b/src/Leanings/penchants.jl index 3a5186256..34f457cb4 100644 --- a/src/Leanings/penchants.jl +++ b/src/Leanings/penchants.jl @@ -1,6 +1,4 @@ -x = [0, 0, 1, 0, 0, 1, 0, 0, 1, 0] -y = [0, 0, 0, 1, 0, 0, 1, 0, 0, 1] export mean_observed_penchant, mean_leaning using Statistics From 1b0cd8de27d9e986a294c3f13b12d864b1a3290d Mon Sep 17 00:00:00 2001 From: Kristian Haaga Date: Wed, 1 Dec 2021 06:14:13 +0100 Subject: [PATCH 06/29] Add return statement --- src/Leanings/penchants.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Leanings/penchants.jl b/src/Leanings/penchants.jl index 34f457cb4..89e56940f 100644 --- a/src/Leanings/penchants.jl +++ b/src/Leanings/penchants.jl @@ -88,7 +88,7 @@ function mean_observed_ρ(ρs, κs, ns_ec, L, weighted = false) end if weighted - ρ / L + return ρ / L else return ρ / ct end From 0310e53ba9dfe861a7787c89aad1ea7f84120ef4 Mon Sep 17 00:00:00 2001 From: Kristian Haaga Date: Wed, 1 Dec 2021 06:14:31 +0100 Subject: [PATCH 07/29] Add data requirement to docstring --- src/Leanings/penchants.jl | 22 +++++++++++++++------- 1 file changed, 15 insertions(+), 7 deletions(-) diff --git a/src/Leanings/penchants.jl b/src/Leanings/penchants.jl index 89e56940f..48a5ffae0 100644 --- a/src/Leanings/penchants.jl +++ b/src/Leanings/penchants.jl @@ -130,18 +130,26 @@ end Compute the *mean observed leaning* `ρ̄` (McCracken & Weigel, 2016)[^McCrackenWeigel2016] between scalar time series `x` and `y`, based on the standard `l`-assignment of the `l`-assignment ``\\{C, E\\} = \\{x, y\\}, \\{x_{t-l}, y_t\\}``. -If `ρ̄ > 0`, then the probability that `x` drives `y` is higher than the probability -that `y` drives `x`. Vice versa, if `ρ̄ > 0`, then the probability that `y` drives `x` is higher than the probability -that `x` drives `y`. - If `weighted == true`, then each penchant is weighted by the number of times it appears in the data, which yields the *weighted mean observed leaning*. This the default behavior, since "...[the weighted mean observed leaning] accounts for the frequency of observed cause-effect pairs within the data, which is assumed to be a predictor of causal influence" (McCracken & Weigel, 2016). -This implementation *does not* take care of the tolerance domains discussed in the paper. -Pre-process your time series, using appropriate binning or symbolization schemes, -so that `x` and `y` both consists of discrete symbols. +## Intepretation + +If `ρ̄ > 0`, then the probability that `x` drives `y` is higher than the probability +that `y` drives `x`. Vice versa, if `ρ̄ > 0`, then the probability that `y` drives `x` is higher than the probability +that `x` drives `y`. + +## Data requirements + +Both `x` and `y` *must* be discrete, and the number discrete states can't be too high +relative to the number of points in `x` and `y`. Otherwise, probabilities cannot be reliably estimated. +If the function returns `NaN`, then you probably either haven't discretized your time series, +or the partition is too fine-grained given the numbers of points available. + +Hence, this implementation *does not* take care of the tolerance domains discussed in the paper. +Pre-process your time series using appropriate binning or symbolization schemes. [^McCrackenWeigel2016]: McCracken, J. M., & Weigel, R. S. (2016). Nonparametric causal inference for bivariate time series. Physical Review E, 93(2), 022207. """ From 92b561e24aa492c888328195a5177ca92b465840 Mon Sep 17 00:00:00 2001 From: Kristian Haaga Date: Wed, 1 Dec 2021 06:30:21 +0100 Subject: [PATCH 08/29] `weighted` is a keyword --- src/Leanings/penchants.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/Leanings/penchants.jl b/src/Leanings/penchants.jl index 48a5ffae0..b61406291 100644 --- a/src/Leanings/penchants.jl +++ b/src/Leanings/penchants.jl @@ -153,7 +153,8 @@ Pre-process your time series using appropriate binning or symbolization schemes. [^McCrackenWeigel2016]: McCracken, J. M., & Weigel, R. S. (2016). Nonparametric causal inference for bivariate time series. Physical Review E, 93(2), 022207. """ -function mean_leaning(x, y, l = 1, weighted = true) +function mean_leaning(x, y, l = 1; weighted = true) return mean_observed_penchant(x, y, l, weighted = weighted) - mean_observed_penchant(y, x, l, weighted = weighted) -end \ No newline at end of file +end + From 944816a2d14e1f4343779eadc286711698533dee Mon Sep 17 00:00:00 2001 From: Kristian Haaga Date: Wed, 1 Dec 2021 06:40:33 +0100 Subject: [PATCH 09/29] Update test --- test/methods/test_leanings.jl | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/test/methods/test_leanings.jl b/test/methods/test_leanings.jl index 1f5fea0e1..659f754f3 100644 --- a/test/methods/test_leanings.jl +++ b/test/methods/test_leanings.jl @@ -1,11 +1,11 @@ -x = [0, 0, 1, 0, 0, 1, 0, 0, 1, 0] -y = [0, 0, 0, 1, 0, 0, 1, 0, 0, 1] - @testset "Leanings" begin - @test mean_observed_penchant(x, y, 1, weighted = false) ≈ 1.0 - @test mean_observed_penchant(y, x, 1, weighted = false) ≈ 1/7 - @test mean_observed_penchant(x, y, 1, weighted = true) ≈ 1.0 - @test mean_observed_penchant(y, x, 1, weighted = true) ≈ 3/63 + x = [0, 0, 1, 0, 0, 1, 0, 0, 1, 0] + y = [0, 0, 0, 1, 0, 0, 1, 0, 0, 1] + @test penchant(x, y, 1, weighted = false) ≈ 1.0 + @test penchant(y, x, 1, weighted = false) ≈ 1/7 + + @test penchant(x, y, 1, weighted = true) ≈ 1.0 + @test penchant(y, x, 1, weighted = true) ≈ 3/63 end \ No newline at end of file From ab7447ca3d2a541f83762f6492991aca529ea57e Mon Sep 17 00:00:00 2001 From: Kristian Haaga Date: Wed, 1 Dec 2021 06:40:51 +0100 Subject: [PATCH 10/29] Simplify and shorten function names --- src/Leanings/penchants.jl | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/Leanings/penchants.jl b/src/Leanings/penchants.jl index b61406291..e4e16a4b7 100644 --- a/src/Leanings/penchants.jl +++ b/src/Leanings/penchants.jl @@ -1,5 +1,5 @@ -export mean_observed_penchant, mean_leaning +export penchant, lean using Statistics @@ -95,7 +95,7 @@ function mean_observed_ρ(ρs, κs, ns_ec, L, weighted = false) end """ - mean_observed_penchant(x, y, l, weighted = false) → ρ̄ + penchant(x, y, l; weighted = false) → ρ̄ Computes the *mean observed penchant* of the `l`-assignment ``\\{C, E\\} = \\{x, y\\}, \\{x_{t-l}, y_t\\}``. @@ -103,7 +103,7 @@ If `weighted == true`, then compute the *weighted mean observed penchant*. [^McCrackenWeigel2016]: McCracken, J. M., & Weigel, R. S. (2016). Nonparametric causal inference for bivariate time series. Physical Review E, 93(2), 022207. """ -function mean_observed_penchant(x, y, l = 1; weighted = false) +function penchant(x, y, l = 1; weighted = false) n = length(x) ns_c, ns_e, ns_ec = penchant_counts(x, y, l) @@ -125,7 +125,7 @@ end """ - mean_leaning(x, y, l = 1, weighted = true) → ρ̄ + lean(x, y, l = 1; weighted = true) → ρ̄ Compute the *mean observed leaning* `ρ̄` (McCracken & Weigel, 2016)[^McCrackenWeigel2016] between scalar time series `x` and `y`, based on the standard `l`-assignment of the `l`-assignment ``\\{C, E\\} = \\{x, y\\}, \\{x_{t-l}, y_t\\}``. @@ -153,8 +153,8 @@ Pre-process your time series using appropriate binning or symbolization schemes. [^McCrackenWeigel2016]: McCracken, J. M., & Weigel, R. S. (2016). Nonparametric causal inference for bivariate time series. Physical Review E, 93(2), 022207. """ -function mean_leaning(x, y, l = 1; weighted = true) - return mean_observed_penchant(x, y, l, weighted = weighted) - - mean_observed_penchant(y, x, l, weighted = weighted) +function lean(x, y, l = 1; weighted = true) + return penchant(x, y, l, weighted = weighted) - + penchant(y, x, l, weighted = weighted) end From 1b8d401155a78debb86e4ffdea5049d1ca7cb92c Mon Sep 17 00:00:00 2001 From: Kristian Haaga Date: Wed, 1 Dec 2021 06:53:55 +0100 Subject: [PATCH 11/29] Docstring fix --- src/Leanings/penchants.jl | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/Leanings/penchants.jl b/src/Leanings/penchants.jl index e4e16a4b7..26d2f6b0f 100644 --- a/src/Leanings/penchants.jl +++ b/src/Leanings/penchants.jl @@ -97,10 +97,13 @@ end """ penchant(x, y, l; weighted = false) → ρ̄ -Computes the *mean observed penchant* of the `l`-assignment ``\\{C, E\\} = \\{x, y\\}, \\{x_{t-l}, y_t\\}``. +Computes the *mean observed penchant* (McCracken & Weigel, 2016)[^McCrackenWeigel2016] +of the `l`-assignment ``\\{\\text{cause}, \\text{effect}\\} =\\{x_{t-l}, y_t\\}``. If `weighted == true`, then compute the *weighted mean observed penchant*. +See also data requirements discussed in [`lean`](@ref). + [^McCrackenWeigel2016]: McCracken, J. M., & Weigel, R. S. (2016). Nonparametric causal inference for bivariate time series. Physical Review E, 93(2), 022207. """ function penchant(x, y, l = 1; weighted = false) @@ -128,7 +131,7 @@ end lean(x, y, l = 1; weighted = true) → ρ̄ Compute the *mean observed leaning* `ρ̄` (McCracken & Weigel, 2016)[^McCrackenWeigel2016] between scalar time series `x` and `y`, -based on the standard `l`-assignment of the `l`-assignment ``\\{C, E\\} = \\{x, y\\}, \\{x_{t-l}, y_t\\}``. +based on the standard `l`-assignment of the `l`-assignment ``\\{\\text{cause}, \\text{effect}\\} =\\{x_{t-l}, y_t\\}``. If `weighted == true`, then each penchant is weighted by the number of times it appears in the data, which yields the *weighted mean observed leaning*. This the default behavior, since From 28920fc746f48993f423d0917293a4f36fe381ca Mon Sep 17 00:00:00 2001 From: Kristian Haaga Date: Wed, 1 Dec 2021 06:55:46 +0100 Subject: [PATCH 12/29] Change function names in docs --- docs/src/leanings.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/src/leanings.md b/docs/src/leanings.md index 99587e144..543b12144 100644 --- a/docs/src/leanings.md +++ b/docs/src/leanings.md @@ -4,6 +4,6 @@ Leanings are probability-based exploratory causal inference measures that rely o implementations of the *penchants and leanings* algoriths from McCracken and Weigel (2016). ```@docs -mean_observed_penchant -mean_leaning +penchant +lean ``` \ No newline at end of file From 7c5d6cc9a32193cde3a00fb0ca6c0aa38f9be7d9 Mon Sep 17 00:00:00 2001 From: Kristian Haaga Date: Wed, 1 Dec 2021 06:56:09 +0100 Subject: [PATCH 13/29] Update version --- Project.toml | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index 0f3f19575..7464bda9c 100644 --- a/Project.toml +++ b/Project.toml @@ -41,10 +41,9 @@ Neighborhood = "0.2.3" Reexport = "0.2, 1" Requires = "^1" SimpleDiffEq = "^1" -StaticArrays = "^1" StatsBase = "^0.33" TimeseriesSurrogates = "^1.2" -TransferEntropy = "^1.3.2" +TransferEntropy = "^1.3.3" Wavelets = "0.9" julia = "^1.6" From c790973cba7ac5c6b8d80b6b1f91c4255026c222 Mon Sep 17 00:00:00 2001 From: Kristian Haaga Date: Wed, 1 Dec 2021 10:09:44 +0100 Subject: [PATCH 14/29] Exclude failing tests for now (they are not relevant for this PR) T --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index 1c3c61e7a..ad534b9d5 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -6,7 +6,7 @@ import DynamicalSystemsBase: DiscreteDynamicalSystem, ContinuousDynamicalSystem include("methods/test_smeasure.jl") include("methods/test_joint_distance_distribution.jl") -include("methods/test_predictive_asymmetry.jl") +#include("methods/test_predictive_asymmetry.jl") include("methods/test_crossmapping.jl") include("methods/test_leanings.jl") From 78ea28f54825fccf16927cb572d8874619e925df Mon Sep 17 00:00:00 2001 From: Kristian Haaga Date: Wed, 1 Dec 2021 10:10:03 +0100 Subject: [PATCH 15/29] Rename variables for clarity --- src/Leanings/penchants.jl | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/src/Leanings/penchants.jl b/src/Leanings/penchants.jl index 26d2f6b0f..62530902a 100644 --- a/src/Leanings/penchants.jl +++ b/src/Leanings/penchants.jl @@ -34,7 +34,7 @@ function penchant_counts(x, y, l = 1) penchants = get_penchants(x, y) @inbounds for penchant in penchants - sx, sy = penchant[1], penchant[2] + state_cause, state_effect = penchant[1], penchant[2] # Number of times the assumed cause has appeared n_c = 0 @@ -51,16 +51,15 @@ function penchant_counts(x, y, l = 1) effect = y[t] cause = x[t-l] - if effect == sx && cause == sy + if effect == state_effect && cause == state_cause n_ec += 1 end - - if effect == sx + if effect == state_effect n_e += 1 end - if cause == sy + if cause == state_cause n_c += 1 end end From 16c06d60b1993f6b7e830d5614037f9cabd08772 Mon Sep 17 00:00:00 2001 From: Kristian Haaga Date: Wed, 1 Dec 2021 10:10:12 +0100 Subject: [PATCH 16/29] Documentation improvements --- src/Leanings/penchants.jl | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/src/Leanings/penchants.jl b/src/Leanings/penchants.jl index 62530902a..d784dcaff 100644 --- a/src/Leanings/penchants.jl +++ b/src/Leanings/penchants.jl @@ -129,20 +129,19 @@ end lean(x, y, l = 1; weighted = true) → ρ̄ -Compute the *mean observed leaning* `ρ̄` (McCracken & Weigel, 2016)[^McCrackenWeigel2016] between scalar time series `x` and `y`, -based on the standard `l`-assignment of the `l`-assignment ``\\{\\text{cause}, \\text{effect}\\} =\\{x_{t-l}, y_t\\}``. +Compute the *mean observed leaning* `ρ̄` (McCracken & Weigel, 2016)[^McCrackenWeigel2016] between scalar time series `x` and `y`. + +If `ρ̄ > 0`, then the probability that `x` drives `y` is higher than the probability +that `y` drives `x`. Vice versa, if `ρ̄ > 0`, then the probability that `y` drives `x` is higher than the probability +that `x` drives `y`. + +## Weighting If `weighted == true`, then each penchant is weighted by the number of times it appears in the data, which yields the *weighted mean observed leaning*. This the default behavior, since "...[the weighted mean observed leaning] accounts for the frequency of observed cause-effect pairs within the data, which is assumed to be a predictor of causal influence" (McCracken & Weigel, 2016). -## Intepretation - -If `ρ̄ > 0`, then the probability that `x` drives `y` is higher than the probability -that `y` drives `x`. Vice versa, if `ρ̄ > 0`, then the probability that `y` drives `x` is higher than the probability -that `x` drives `y`. - ## Data requirements Both `x` and `y` *must* be discrete, and the number discrete states can't be too high From 7b74daf37abe3724d0f1b311ceb91038d7916c9a Mon Sep 17 00:00:00 2001 From: Kristian Haaga Date: Wed, 1 Dec 2021 10:15:31 +0100 Subject: [PATCH 17/29] Specify possible ranges for penchants and leanings --- src/Leanings/penchants.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Leanings/penchants.jl b/src/Leanings/penchants.jl index d784dcaff..cc91af65e 100644 --- a/src/Leanings/penchants.jl +++ b/src/Leanings/penchants.jl @@ -96,7 +96,7 @@ end """ penchant(x, y, l; weighted = false) → ρ̄ -Computes the *mean observed penchant* (McCracken & Weigel, 2016)[^McCrackenWeigel2016] +Computes the *mean observed penchant* `ρ̄ ∈ [-1, 1]` (McCracken & Weigel, 2016)[^McCrackenWeigel2016] of the `l`-assignment ``\\{\\text{cause}, \\text{effect}\\} =\\{x_{t-l}, y_t\\}``. If `weighted == true`, then compute the *weighted mean observed penchant*. @@ -129,7 +129,7 @@ end lean(x, y, l = 1; weighted = true) → ρ̄ -Compute the *mean observed leaning* `ρ̄` (McCracken & Weigel, 2016)[^McCrackenWeigel2016] between scalar time series `x` and `y`. +Compute the *mean observed leaning* `λ ∈ [-2, 2]` (McCracken & Weigel, 2016)[^McCrackenWeigel2016] between scalar time series `x` and `y`. If `ρ̄ > 0`, then the probability that `x` drives `y` is higher than the probability that `y` drives `x`. Vice versa, if `ρ̄ > 0`, then the probability that `y` drives `x` is higher than the probability From d916cd70ec5fed55bdb9123b1cde2f6b9f1be84e Mon Sep 17 00:00:00 2001 From: Kristian Haaga Date: Wed, 1 Dec 2021 10:15:39 +0100 Subject: [PATCH 18/29] Add test for lean --- test/methods/test_leanings.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/test/methods/test_leanings.jl b/test/methods/test_leanings.jl index 659f754f3..a3990f3ec 100644 --- a/test/methods/test_leanings.jl +++ b/test/methods/test_leanings.jl @@ -8,4 +8,6 @@ @test penchant(x, y, 1, weighted = true) ≈ 1.0 @test penchant(y, x, 1, weighted = true) ≈ 3/63 + + @test lean(x, y, 1, weighted = true) ≈ 60/63 end \ No newline at end of file From bfec698155aeae316ca45b7a558790de3843da0e Mon Sep 17 00:00:00 2001 From: Kristian Haaga Date: Wed, 1 Dec 2021 10:29:49 +0100 Subject: [PATCH 19/29] Definition of causal penchant --- src/Leanings/penchants.jl | 31 +++++++++++++++++++++++++++---- 1 file changed, 27 insertions(+), 4 deletions(-) diff --git a/src/Leanings/penchants.jl b/src/Leanings/penchants.jl index cc91af65e..9001e2f92 100644 --- a/src/Leanings/penchants.jl +++ b/src/Leanings/penchants.jl @@ -96,12 +96,35 @@ end """ penchant(x, y, l; weighted = false) → ρ̄ -Computes the *mean observed penchant* `ρ̄ ∈ [-1, 1]` (McCracken & Weigel, 2016)[^McCrackenWeigel2016] -of the `l`-assignment ``\\{\\text{cause}, \\text{effect}\\} =\\{x_{t-l}, y_t\\}``. +Computes the *mean observed penchant* `ρ̄` (McCracken & Weigel, 2016)[^McCrackenWeigel2016] +of the `l`-assignment ``\\{\\text{C}, \\text{E}\\} = \\{\\text{cause}, \\text{effect}\\} =\\{x_{t-l}, y_t\\}``. If `weighted == true`, then compute the *weighted mean observed penchant*. -See also data requirements discussed in [`lean`](@ref). +## Definition + +The *causal penchant* , or causal tendency, is a causal indicator defined as +``\\rho_{CE} = P(E|C) - P(E|\\bar{C})``, where ``\\rho_{CE} \\in [-1, 1]`` (McCracken & Weigel, 2016). +If ``\\rho_{CE} > 0``, then ``C`` causes or drives ``E``, and if +``\\rho_{CE} \\leq 0`, then the direction of influence is undetermined. + +A direct formula for ``\\rho_{CE}`` can be obtained using Bayes' theorem. + +```math +P(E|C) = P(C|E) \\dfrac{P(E)}{P(C)} +``` +Using the definitions of probability complements, one arrives at the following +expression (see the original paper for a detailed derivation): + +```math +\\rho_{CE} = P(E|C) \\left[1 + \\dfrac{P(C)}{1-P(C)} \\right] - \\dfrac{P(E)}{1 - P(C)}. +``` + +Applying appropriate discretization schemes, these probabilities can be estimated +directly from time series using simple counting, which makes the method fast and +well-suited for exploratory causal inference analysis. + +See also notes on data requirements discussed in [`lean`](@ref). [^McCrackenWeigel2016]: McCracken, J. M., & Weigel, R. S. (2016). Nonparametric causal inference for bivariate time series. Physical Review E, 93(2), 022207. """ @@ -157,5 +180,5 @@ Pre-process your time series using appropriate binning or symbolization schemes. function lean(x, y, l = 1; weighted = true) return penchant(x, y, l, weighted = weighted) - penchant(y, x, l, weighted = weighted) -end +end From 2d199903a2f2e85e018b95a3d7569cf12a81d0a9 Mon Sep 17 00:00:00 2001 From: Kristian Haaga Date: Wed, 1 Dec 2021 10:45:44 +0100 Subject: [PATCH 20/29] Rephrase docs --- src/Leanings/penchants.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Leanings/penchants.jl b/src/Leanings/penchants.jl index 9001e2f92..3e8436515 100644 --- a/src/Leanings/penchants.jl +++ b/src/Leanings/penchants.jl @@ -124,7 +124,7 @@ Applying appropriate discretization schemes, these probabilities can be estimate directly from time series using simple counting, which makes the method fast and well-suited for exploratory causal inference analysis. -See also notes on data requirements discussed in [`lean`](@ref). +See also [`lean`](@ref) and data requirements discussed therein. [^McCrackenWeigel2016]: McCracken, J. M., & Weigel, R. S. (2016). Nonparametric causal inference for bivariate time series. Physical Review E, 93(2), 022207. """ From 4c9fbeb935118a88f3fc6924816ee818a8b5375a Mon Sep 17 00:00:00 2001 From: Kristian Haaga Date: Wed, 1 Dec 2021 10:47:24 +0100 Subject: [PATCH 21/29] All tests involving Dataset construction fails, reintroduce the file --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index ad534b9d5..1c3c61e7a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -6,7 +6,7 @@ import DynamicalSystemsBase: DiscreteDynamicalSystem, ContinuousDynamicalSystem include("methods/test_smeasure.jl") include("methods/test_joint_distance_distribution.jl") -#include("methods/test_predictive_asymmetry.jl") +include("methods/test_predictive_asymmetry.jl") include("methods/test_crossmapping.jl") include("methods/test_leanings.jl") From 096fabab88b4e4b24c2d387b216ae8a061af35a2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20Agas=C3=B8ster=20Haaga?= Date: Wed, 1 Dec 2021 15:36:27 +0100 Subject: [PATCH 22/29] Trigger new CI run --- test/methods/test_leanings.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/test/methods/test_leanings.jl b/test/methods/test_leanings.jl index a3990f3ec..8cd1cb682 100644 --- a/test/methods/test_leanings.jl +++ b/test/methods/test_leanings.jl @@ -1,6 +1,5 @@ @testset "Leanings" begin - x = [0, 0, 1, 0, 0, 1, 0, 0, 1, 0] y = [0, 0, 0, 1, 0, 0, 1, 0, 0, 1] @test penchant(x, y, 1, weighted = false) ≈ 1.0 From f7c35911d7d4e999a318c7fe35bbd7731533b792 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20Agas=C3=B8ster=20Haaga?= Date: Wed, 1 Dec 2021 15:39:52 +0100 Subject: [PATCH 23/29] Require DelayEmbeddings 2.0.2 (breaks on Julia 1.7 otherwise) --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 7464bda9c..683d8f348 100644 --- a/Project.toml +++ b/Project.toml @@ -29,7 +29,7 @@ Wavelets = "29a6e085-ba6d-5f35-a997-948ac2efa89a" [compat] ChaosTools = "^2.0" -DelayEmbeddings = "^2.0" +DelayEmbeddings = "^2.0.2" Distances = "^0.10" Distributions = "^0.24, 0.25" DynamicalSystemsBase = "^2.0" From e21791f51a458da0ac4019d0565c7c4668ebbbb7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20Agas=C3=B8ster=20Haaga?= Date: Wed, 1 Dec 2021 15:41:55 +0100 Subject: [PATCH 24/29] Change heading in menu --- docs/src/leanings.md | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/docs/src/leanings.md b/docs/src/leanings.md index 543b12144..c95a1409f 100644 --- a/docs/src/leanings.md +++ b/docs/src/leanings.md @@ -1,9 +1,10 @@ -# Leanings +# Penchants and leanings -Leanings are probability-based exploratory causal inference measures that rely on the estimation of probabilities directly from the data. Here, we provide -implementations of the *penchants and leanings* algoriths from McCracken and Weigel (2016). +Penchants and leanings are probability-based exploratory causal inference measures that rely on the estimation of probabilities directly from the data. Here, we provide implementations of the *penchants and leanings* +algoriths from McCracken and Weigel (2016). ```@docs penchant lean -``` \ No newline at end of file +``` + From 2d29f3cc42b1996a517671060997667b04bb721b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20Agas=C3=B8ster=20Haaga?= Date: Wed, 1 Dec 2021 15:42:03 +0100 Subject: [PATCH 25/29] Latex typo --- src/Leanings/penchants.jl | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/Leanings/penchants.jl b/src/Leanings/penchants.jl index 3e8436515..728e18107 100644 --- a/src/Leanings/penchants.jl +++ b/src/Leanings/penchants.jl @@ -106,7 +106,7 @@ If `weighted == true`, then compute the *weighted mean observed penchant*. The *causal penchant* , or causal tendency, is a causal indicator defined as ``\\rho_{CE} = P(E|C) - P(E|\\bar{C})``, where ``\\rho_{CE} \\in [-1, 1]`` (McCracken & Weigel, 2016). If ``\\rho_{CE} > 0``, then ``C`` causes or drives ``E``, and if -``\\rho_{CE} \\leq 0`, then the direction of influence is undetermined. +``\\rho_{CE} \\leq 0``, then the direction of influence is undetermined. A direct formula for ``\\rho_{CE}`` can be obtained using Bayes' theorem. @@ -152,11 +152,12 @@ end lean(x, y, l = 1; weighted = true) → ρ̄ -Compute the *mean observed leaning* `λ ∈ [-2, 2]` (McCracken & Weigel, 2016)[^McCrackenWeigel2016] between scalar time series `x` and `y`. +Compute the *mean observed leaning* `λ ∈ [-2, 2]` (McCracken & Weigel, 2016)[^McCrackenWeigel2016] +between scalar time series `x` and `y`. If `ρ̄ > 0`, then the probability that `x` drives `y` is higher than the probability -that `y` drives `x`. Vice versa, if `ρ̄ > 0`, then the probability that `y` drives `x` is higher than the probability -that `x` drives `y`. +that `y` drives `x`. Vice versa, if `ρ̄ > 0`, then the probability that `y` drives `x` +is higher than the probability that `x` drives `y`. ## Weighting From 9bdb54b77fa577280248e3b61932c983560b3566 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20Agas=C3=B8ster=20Haaga?= Date: Wed, 1 Dec 2021 16:21:33 +0100 Subject: [PATCH 26/29] Add penchants and leanings to overview --- docs/src/index.md | 4 ++++ docs/src/leanings.md | 2 +- 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/docs/src/index.md b/docs/src/index.md index 2c927bfe3..42281bbcd 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -35,6 +35,10 @@ We also provide estimators for generalized entropy and mutual information, thoug - Generalized entropy ([`genentropy`](@ref)) and fast histograms ([`binhist`](@ref)) - Mutual information ([`mutualinfo`](@ref)) +### Exploratory methods + +- [Penchants and leanings](@ref leanings). + ## Contributions/issues Do you want additional methods or example systems to be implemented? Make a PR to the diff --git a/docs/src/leanings.md b/docs/src/leanings.md index c95a1409f..b75ca0aa0 100644 --- a/docs/src/leanings.md +++ b/docs/src/leanings.md @@ -1,4 +1,4 @@ -# Penchants and leanings +# [Penchants and leanings](@id leanings) Penchants and leanings are probability-based exploratory causal inference measures that rely on the estimation of probabilities directly from the data. Here, we provide implementations of the *penchants and leanings* algoriths from McCracken and Weigel (2016). From 235957acc5dd24a40aaee4c35da36a5feb9a6138 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20Agas=C3=B8ster=20Haaga?= Date: Wed, 1 Dec 2021 17:35:49 +0100 Subject: [PATCH 27/29] `weighted` should be true by default for both `lean` and `penchant`` --- src/Leanings/penchants.jl | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/Leanings/penchants.jl b/src/Leanings/penchants.jl index 728e18107..662792d7b 100644 --- a/src/Leanings/penchants.jl +++ b/src/Leanings/penchants.jl @@ -128,7 +128,7 @@ See also [`lean`](@ref) and data requirements discussed therein. [^McCrackenWeigel2016]: McCracken, J. M., & Weigel, R. S. (2016). Nonparametric causal inference for bivariate time series. Physical Review E, 93(2), 022207. """ -function penchant(x, y, l = 1; weighted = false) +function penchant(x, y, l = 1; weighted = true) n = length(x) ns_c, ns_e, ns_ec = penchant_counts(x, y, l) @@ -179,7 +179,6 @@ Pre-process your time series using appropriate binning or symbolization schemes. [^McCrackenWeigel2016]: McCracken, J. M., & Weigel, R. S. (2016). Nonparametric causal inference for bivariate time series. Physical Review E, 93(2), 022207. """ function lean(x, y, l = 1; weighted = true) - return penchant(x, y, l, weighted = weighted) - - penchant(y, x, l, weighted = weighted) + return penchant(x, y, l, weighted = weighted) - penchant(y, x, l, weighted = weighted) end From 9042531665cc54275e2b7370b9cb65ad3759f0a1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20Agas=C3=B8ster=20Haaga?= Date: Wed, 1 Dec 2021 18:30:10 +0100 Subject: [PATCH 28/29] example --- src/Leanings/example.jl | 57 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 57 insertions(+) create mode 100644 src/Leanings/example.jl diff --git a/src/Leanings/example.jl b/src/Leanings/example.jl new file mode 100644 index 000000000..6c8201ff9 --- /dev/null +++ b/src/Leanings/example.jl @@ -0,0 +1,57 @@ +using CausalityTools, DynamicalSystems, Plots, Statistics, Distributions; gr() + +npts, Ttr = 10000, 5000 +x, y = columns(trajectory(henon2(c_xy = 0.0), npts - 1, Ttr = Ttr)) +xr, yr = rand(Uniform(-1, 1), npts), rand(Uniform(-1, 1), npts) +plot() +plot!(x, label = "x") +plot!(y, label = "y") + +import Entropies: symbolize, SymbolicPermutation + +m = 2; τ = -1 +n = 5 +weighted = false +nreps = 50 +ls = 0:10 +leans_wn = zeros(nreps, length(ls)) +leans_henon_c0 = zeros(nreps, length(ls)) +leans_henon_c2 = zeros(nreps, length(ls)) + +function partition(x, n::Int) + xmin = minimum(x) # the reference point + Δx = (maximum(x) - xmin) / n + [floor(Int, (xᵢ - xmin) / Δx) for xᵢ in x] +end + +for i = 1:nreps + xr, yr = rand(Uniform(-1, 1), npts), rand(Uniform(-1, 1), npts) + X = partition(xr, n)#symbolize(xr, SymbolicPermutation(m = m, τ = τ)) + Y = partition(yr, n)#symbolize(yr, SymbolicPermutation(m = m, τ = τ)) + for (j, l) in enumerate(ls) + leans_wn[i, j] = lean(X, Y, l, weighted = weighted) + end + + x, y = columns(trajectory(henon2(c_xy = 0.0), npts - 1, Ttr = Ttr)) + X = partition(x, n)#symbolize(x, SymbolicPermutation(m = m, τ = τ)) + Y = partition(y, n)#symbolize(y, SymbolicPermutation(m = m, τ = τ)) + for (j, l) in enumerate(ls) + leans_henon_c0[i, j] = lean(X, Y, l, weighted = weighted) + end + + x, y = columns(trajectory(henon2(c_xy = 2.0), npts - 1, Ttr = Ttr)) + X = partition(x, n)#symbolize(x, SymbolicPermutation(m = m, τ = τ)) + Y = partition(y, n)#symbolize(y, SymbolicPermutation(m = m, τ = τ)) + for (j, l) in enumerate(ls) + leans_henon_c2[i, j] = lean(X, Y, l, weighted = weighted) + end +end + +mean_leans_wn = dropdims(mean(leans_wn, dims = 1), dims = 1) +mean_leans_henon_c0 = dropdims(mean(leans_henon_c0, dims = 1), dims = 1) +mean_leans_henon_c2 = dropdims(mean(leans_henon_c2, dims = 1), dims = 1) + +plot(xlabel = "l", ylabel = "λx→y", ylims = (-1.1, 1.1)) +#plot!(ls, mean_leans_wn, label = "white noise", marker = :star) +plot!(ls, mean_leans_henon_c0, label = "henon (c_xy = 0.0)", marker = :circle) +#plot!(ls, mean_leans_henon_c2, label = "henon (c_xy = 2.0)", marker = :star5) \ No newline at end of file From ce542deb39d26083c01848844fc6d6d1f63c5de8 Mon Sep 17 00:00:00 2001 From: Kristian Haaga Date: Mon, 6 Dec 2021 15:58:28 +0100 Subject: [PATCH 29/29] Add example system --- .../discretemaps/nonlinear_sindriver_2d.jl | 34 +++++++++++++++++++ 1 file changed, 34 insertions(+) create mode 100644 src/example_systems/discretemaps/nonlinear_sindriver_2d.jl diff --git a/src/example_systems/discretemaps/nonlinear_sindriver_2d.jl b/src/example_systems/discretemaps/nonlinear_sindriver_2d.jl new file mode 100644 index 000000000..ead84e60a --- /dev/null +++ b/src/example_systems/discretemaps/nonlinear_sindriver_2d.jl @@ -0,0 +1,34 @@ +function eom_nonlinear_sindriver_2d(dx, x, p, n) + a, b, c, t, Δt = (p...,) + x, y = x[1], x[2] + 𝒩 = Normal(0, 1) + + dx[1] = sin(t) + dx[2] = a*x * (1 - b*x) + c*rand(𝒩) + p[end-1] += Δt # update t + + return +end + +""" + nonlinear_sindriver_2d(;u₀ = 0, a = 1.0, b = 1.0, c = 2.0, Δt = 1) + +A discrete, nonlinear 2D system from McCracken & Weigel (2016)[^McCrackenWeigel2016], +governed by the following equations of motion: + +```math +\\begin{aligned} +x_{t+1} = \\sin(t) \\\\ +y(t+1) &= a x_t (1 - B x_t) + C \\eta_t +\\end{aligned} +``` + +where ``\\eta_t \\sim \\mathcal{N}(\\mu = 0, \\sigma = 1)``, ``A, B, C \\in [0, 1]``. + +[^McCrackenWeigel2016]: McCracken, J. M., & Weigel, R. S. (2016). Nonparametric causal inference for bivariate time series. Physical Review E, 93(2), 022207. +""" +function nonlinear_sindriver_2d(;u₀ = 0, a = 1.0, b = 1.0, c = 2.0, Δt = (1/30*π)) + @assert a >= 0; @assert b >= 0; @assert c >= 0 + + DiscreteDynamicalSystem(eom_nonlinear_sindriver2d, [0, u₀], [a, b, c, 0, Δt]) +end \ No newline at end of file