|
| 1 | +# Pairwise asymmetric inference |
| 2 | + |
| 3 | +```@docs |
| 4 | +pai |
| 5 | +``` |
| 6 | + |
| 7 | +## Example: nonlinear system |
| 8 | + |
| 9 | +Let's try to reproduce figure 8 in McCracken & Weigel (2014)[^McCracken2014]. We'll start by defining the their example B (equations 6-7). This system consists of two |
| 10 | +variables ``X`` and ``Y``, where ``X`` drives ``Y``. |
| 11 | + |
| 12 | +```@example pai_ex |
| 13 | +using CausalityTools, DynamicalSystems, Plots, StatsBase, Statistics, Distributions; gr() |
| 14 | +
|
| 15 | +function eom_nonlinear_sindriver(dx, x, p, n) |
| 16 | + a, b, c, t, Δt = (p...,) |
| 17 | + x, y = x[1], x[2] |
| 18 | + 𝒩 = Normal(0, 1) |
| 19 | + |
| 20 | + dx[1] = sin(t) |
| 21 | + dx[2] = a*x * (1 - b*x) + c*rand(𝒩) |
| 22 | + p[end-1] += 1 # update t |
| 23 | +
|
| 24 | + return |
| 25 | +end |
| 26 | +
|
| 27 | +function nonlinear_sindriver(;u₀ = rand(2), a = 1.0, b = 1.0, c = 2.0, Δt = 1) |
| 28 | + DiscreteDynamicalSystem(eom_nonlinear_sindriver, u₀, [a, b, c, 0, Δt]) |
| 29 | +end |
| 30 | +
|
| 31 | +# Create a system of nonidentical logistic maps where coupling from variable x to variable y |
| 32 | +# is stronger than vice versa. |
| 33 | +sys = nonlinear_sindriver(a = 1.0, b = 1.0, c = 2.0) |
| 34 | +npts = 100 |
| 35 | +orbit = trajectory(sys, npts, Ttr = 10000) |
| 36 | +x, y = columns(orbit); |
| 37 | +plot(xlabel = "Time step", ylabel = "Value") |
| 38 | +# Standardize and plot data |
| 39 | +plot!((x .- mean(x)) ./ std(x), label = "x") |
| 40 | +plot!((y .- mean(y)) ./ std(y), label = "y") |
| 41 | +savefig("pai_ts.svg") # hide |
| 42 | +``` |
| 43 | + |
| 44 | + |
| 45 | + |
| 46 | +Now, let's generate such time series for many different values of the parameters `a` and `b`, and compute PAI for fixed `p = 2.0`. This will replicate the upper right panel of figure 8 in the original paper. |
| 47 | + |
| 48 | +```@example pai_ex |
| 49 | +as = 0.25:0.25:4.0 |
| 50 | +bs = 0.25:0.25:4.0 |
| 51 | +
|
| 52 | +pai_xys = zeros(length(as), length(bs)) |
| 53 | +pai_yxs = zeros(length(as), length(bs)) |
| 54 | +c = 2.0 |
| 55 | +npts = 2000 |
| 56 | +d, τ = 2, 1 |
| 57 | +for (i, a) in enumerate(as) |
| 58 | + for (j, b) in enumerate(bs) |
| 59 | + s = nonlinear_sindriver(a = a, b = a, c = c) |
| 60 | + X, Y = columns(trajectory(s, npts, Ttr = 10000)) |
| 61 | + # Use the segment bootstrap estimator, take the mean of 50 reps over segments of length L = 200 |
| 62 | + pai_xys[i, j] = pai(X, Y, d, τ, :segment, L = 200, nreps = 50) |> mean |
| 63 | + pai_yxs[i, j] = pai(Y, X, d, τ, :segment, L = 200, nreps = 50) |> mean |
| 64 | + end |
| 65 | +end |
| 66 | +``` |
| 67 | + |
| 68 | +Now that we have computed the PAI in both directions, we define a measure of directionality as the difference between PAI in the ``X \to Y`` direction and in the ``Y \to X`` direction, so that if ``X`` drives ``Y``, then ``\Delta < 0``. |
| 69 | + |
| 70 | +```@example pai_ex |
| 71 | +Δ = pai_xys .- pai_yxs |
| 72 | +
|
| 73 | +clr = cgrad(:magma, categorical = true) |
| 74 | +plot(xlabel = "a", ylabel = "b") |
| 75 | +yticks!((1:length(as), string.(as))) |
| 76 | +xticks!((1:length(bs), string.(bs))) |
| 77 | +heatmap!(Δ, c = clr, logscale = true) |
| 78 | +savefig("heatmap_pai.svg") # hide |
| 79 | +``` |
| 80 | + |
| 81 | + |
| 82 | + |
| 83 | +As expected, ``\Delta < 0`` for all parameter combinations, implying that ``X`` "PAI drives" ``Y``. |
0 commit comments