Public Interface

estimate_TE_from_event_times(
  parameters::CoTETEParameters,
  target_events::Array{<:AbstractFloat},
  source_events::Array{<:AbstractFloat};
  conditioning_events::Array{<:AbstractFloat} = Float32[],
)

Estimate the TE from lists of raw event times.

Examples

This example demonstrates estimating the TE between uncoupled homogeneous Poisson processes. This is covered in section II A of our paper. We first create the source and target processes, each with 10 000 events and with rate 1, before running the estimator.

julia> source = sort(1e4*rand(Int(1e4)));

julia> target = sort(1e4*rand(Int(1e4)));

julia> parameters = CoTETE.CoTETEParameters(l_x = 1, l_y = 1);

julia> TE = CoTETE.estimate_TE_from_event_times(parameters, target, source)
0.0

julia> abs(TE - 0) < 0.05 # For Doctesting purposes
true

We can also try increasing the length of the target and source history embeddings

julia> parameters = CoTETE.CoTETEParameters(l_x = 3, l_y = 3);

julia> TE = CoTETE.estimate_TE_from_event_times(parameters, target, source)
0.0

julia> abs(TE - 0) < 0.1 # For Doctesting purposes
true

Let's try some other options

julia> using Distances: Euclidean

julia> parameters = CoTETE.CoTETEParameters(l_x = 1,
                                            l_y = 2,
                                            k_global = 3,
                                            auto_find_start_and_num_events = false,
                                            start_event = 100,
                                            num_target_events = 5000,
                                            num_samples_ratio = 2.3,
                                            metric = Euclidean(),
                                            transform_to_uniform = true);

julia> TE = CoTETE.estimate_TE_from_event_times(parameters, target, source)
0.0

julia> abs(TE - 0) < 0.1 # For Doctesting purposes
true

The next example applies the estimator to a more complex problem, specifically, the process described as example B in Spinney et. al.. The application of the estimator to this example is covered in section II B of our paper. We create the source process as before. Howevever, the target process is originally created as an homogeneous Poisson process with rate 10, before a thinning algorithm is applied to it, in order to provide the dependence on the source.

function thin_target(source, target, target_rate)
    	   start_index = 1
    	   while target[start_index] < source[1]
           	 start_index += 1
    	   end
    	   target = target[start_index:end]
	       new_target = Float64[]
    	   index_of_last_source = 1
    	   for event in target
               while index_of_last_source < length(source) && source[index_of_last_source + 1] < event
               	     index_of_last_source += 1
               end
               distance_to_last_source = event - source[index_of_last_source]
               lambda = 0.5 + 5exp(-50(distance_to_last_source - 0.5)^2) - 5exp(-50(-0.5)^2)
               if rand() < lambda/target_rate
               	  push!(new_target, event)
               end
           end
    	   return new_target
       end

julia> source = sort(1e4*rand(Int(1e4)));

julia> target = sort(1e4*rand(Int(1e5)));

julia> target = thin_target(source, target, 10);

julia> parameters = CoTETE.CoTETEParameters(l_x = 1, l_y = 1);

julia> TE = CoTETE.estimate_TE_from_event_times(parameters, target, source)
0.5076

julia> abs(TE - 0.5076) < 0.1 # For Doctesting purposes
true

We can also try extending the length of the target embeddings in order to better resolve this dependency (along with some other options)

julia> parameters = CoTETE.CoTETEParameters(l_x = 3,
                                            l_y = 1,
                                            transform_to_uniform = true,
                                            k_global = 7,
                                            num_samples_ratio = 5.0);

julia> TE = CoTETE.estimate_TE_from_event_times(parameters, target, source)
0.5076

julia> abs(TE - 0.5076) < 0.1 # For Doctesting purposes
true

function estimate_TE_and_p_value_from_event_times(
    parameters::CoTETEParameters,
    target_events::Array{<:AbstractFloat},
    source_events::Array{<:AbstractFloat};
    conditioning_events::Array{<:AbstractFloat} = Float32[],
    return_surrogate_TE_values::Bool = false,
)

calculate the TE and the $p$ value of it being statistically different from zero.

Examples

This example demonstrates estimating the TE and $p$ value between uncoupled homogeneous Poisson processes. As the true value of the TE is zero, we expect the $p$ value to be uniformly disributed between zero and one.

We first create the source and target processes, each with 1 000 events and with rate 1, before running the estimator and the surrogate generation procedure.

julia> source = sort(1e3*rand(Int(1e3)));

julia> target = sort(1e3*rand(Int(1e3)));

julia> parameters = CoTETE.CoTETEParameters(l_x = 1, l_y = 1);

julia> TE, p = CoTETE.estimate_TE_and_p_value_from_event_times(parameters, target, source)
(0.0, 0.5)

julia> p > 0.05 # For Doctesting purposes. Should fail every now and then.
true

Lets try some other parameters

julia> parameters = CoTETE.CoTETEParameters(l_x = 1, l_y = 1, transform_to_uniform = true, k_perm = 20, sampling_method = "jittered_target");

julia> TE, p = CoTETE.estimate_TE_and_p_value_from_event_times(parameters, target, source)
(0.0, 0.5)

julia> p > 0.05 # For Doctesting purposes. Should fail every now and then.
true

This second example shows using this method on coupled processes for which the true value of the TE is nonzero. As there is a strong coupling between the source and target, we expect the $p$ value to be close to 0. The application of the estimator to this example is covered in section II B of our paper. See the above examples for estimate_TE_from_event_times for more details as well as the implementation of the thinning algorithm.

We create the source process as before. Howevever, the target process is originally created as an homogeneous Poisson process with rate 10, before the thinning algorithm is applied to it, in order to provide the dependence on the source.

DocTestSetup = quote
       function thin_target(source, target, target_rate)
    	   start_index = 1
    	   while target[start_index] < source[1]
           	 start_index += 1
    	   end
    	   target = target[start_index:end]
	       new_target = Float64[]
    	   index_of_last_source = 1
    	   for event in target
               while index_of_last_source < length(source) && source[index_of_last_source + 1] < event
               	     index_of_last_source += 1
               end
               distance_to_last_source = event - source[index_of_last_source]
               lambda = 0.5 + 5exp(-50(distance_to_last_source - 0.5)^2) - 5exp(-50(-0.5)^2)
               if rand() < lambda/target_rate
               	  push!(new_target, event)
               end
           end
    	   return new_target
       end
end

julia> source = sort(1e3*rand(Int(1e3)));

julia> target = sort(1e3*rand(Int(1e4)));

julia> target = thin_target(source, target, 10);

julia> parameters = CoTETE.CoTETEParameters(l_x = 1, l_y = 1);

julia> TE, p = CoTETE.estimate_TE_and_p_value_from_event_times(parameters, target, source)
(0.5, 0.01)

julia> p < 0.05 # For Doctesting purposes. Should fail very rarely.
true

Lets try with some other parameters

julia> parameters = CoTETE.CoTETEParameters(l_x = 3,
                                            l_y = 1,
                                            transform_to_uniform = true,
                                            k_global = 7,
                                            num_samples_ratio = 5.0,
                                            surrogate_num_samples_ratio = 5.0,
                                            sampling_method = "jittered_target",
                                            jittered_sampling_noise = 1.0);

julia> TE, p, locals = CoTETE.estimate_TE_and_p_value_from_event_times(parameters, target, source, return_locals = true);

julia> p < 0.05 # For Doctesting purposes. Should fail very rarely.
true

julia> abs(locals[1]) > 1e-10 # For Doctesting. Check that the locals aren't all zero.
true

estimate_AIS_from_event_times(
  parameters::CoTETEParameters,
  target_events::Array{<:AbstractFloat},
)

Estimate the Active Information Storage (AIS) from lists of raw event times.

See this thesis for a description of AIS.

Examples

This example estimates the AIS on an homogeneous Poisson process. The true value of the AIS on such a process is zero.

julia> target = sort(1e4*rand(Int(1e4)));

julia> parameters = CoTETE.CoTETEParameters(l_x = 1);

julia> AIS = CoTETE.estimate_AIS_from_event_times(parameters, target);

julia> abs(AIS - 0) < 0.05 # For Doctesting purposes
true

estimate_AIS_and_p_value_from_event_times(
  parameters::CoTETEParameters,
  target_events::Array{<:AbstractFloat};
  return_surrogate_AIS_values::Bool = false,
)

Estimate the Active Information Storage (AIS) along with the $p$ value of the AIS being different from 0 from lists of raw event times.

See this thesis for a description of AIS.

Examples

This example estimates the AIS and $p$ value on an homogeneous Poisson process. The true value of the AIS on such a process is zero. We expect the $p$ value to be uniformly distributed between zero and 1.

julia> target = sort(1e3*rand(Int(1e3)));

julia> parameters = CoTETE.CoTETEParameters(l_x = 1);

julia> AIS, p = CoTETE.estimate_AIS_and_p_value_from_event_times(parameters, target)
(0.0, 0.5)

julia> p > 0.05 # For Doctesting purposes. Should fail from time to time
true

This next example estimates the AIS for a process where we know that the AIS must be nonzero. This process has an event occurring every one time unit, with a bit of noise added to the event times.

julia> target = sort(cumsum(ones(Int(1e3))) .+ 1e-2*randn(Int(1e3)));

julia> parameters = CoTETE.CoTETEParameters(l_x = 1);

julia> AIS, p = CoTETE.estimate_AIS_and_p_value_from_event_times(parameters, target)
(1.0, 0.01)

julia> p < 0.05 # For Doctesting purposes. Should fail from time to time
true

struct CoTETEParameters
    l_x::Integer = 0
    l_y::Integer = 0
    l_z::Integer = 0
    auto_find_start_and_num_events::Bool = true
    num_target_events_cap::Integer = -1
    start_event::Integer = 1
    num_target_events::Integer = 0
    num_samples_ratio::AbstractFloat = 1.0
    sampling_method::String = "random_uniform"
    jittered_sampling_noise::AbstractFloat = 5.0
    k_global::Integer = 5
    metric::Metric = Cityblock()
    kraskov_noise_level::AbstractFloat = 1e-8
    transform_to_uniform::Bool = false
    num_average_samples::Int = -1
    num_surrogates::Integer = 100
    surrogate_num_samples_ratio::AbstractFloat = 1.0
    k_perm::Integer = 10
    use_exclusion_windows::Bool = true
    add_dummy_exclusion_windows::Bool = false
end

• l_x::Integer: The number of intervals in the target process to use in the history embeddings. Corresponds to $l_X$ in ^[1].
• l_y::Integer: The number of intervals in the source process to use in the history embeddings. Corresponds to $l_Y$ in ^[1].
• l_z::Integer = 0: The number of intervals in the single conditioning process to use in the history embeddings. Corresponds to $l_{Z_1}$ in ^[1].
  Single conditioning process

  Note that, although the framework developed in our paper considers an arbitrary number of extra conditioning processes, at present the framework can only handle a single such process. This will change in future releases.

• auto_find_start_and_num_events::Bool = true: When set to true, the start event will be set to the first event for which there are sufficient preceding events in all processes such that the embeddings can be constructed. The number of target events will be set such that all time between this first event and the last target event is included.
• num_target_events_cap::Integer = -1: Upper limit on the number of target events to use.
• start_event::Integer = 1: only used when auto_find_start_and_num_events = false. The index of the event in the target process from which to start the analysis.
• num_target_events::Integer = 0: only used when auto_find_start_and_num_events = false. The TE will be calculated on the time series from the timestamp of the start_event-th event of the target process to the timestamp of the start_event + num_target_events-th event of the target process.
• num_samples_ratio::AbstractFloat = 1.0: Controls the number of samples used to estimate the probability density of histories unconditional of the occurrence of events in the target process. This number of samples will be num_samples_ratio * num_target_events. Corresponds to $N_U/N_X$ in ^[1].
• sampling_method::String = "random_uniform": Method with which to place the random sampling points. Setting it to "random_uniform" will place the samples uniformly randomly between the first and last target events. "fixed_interval" will space them at a constant interval. "jittered_target" will copy target spikes and add noise to their timestamps. This last method was developed for this paper. See section IV F of that paper for a discussion of this method.
• jittered_sampling_noise::AbstractFloat = 5.0: Width of the uniform jitter added to the target spike times used in resampling when sampling_method is set to "jittered_target".
• k_global::Integer = 5: The number of nearest neighbours to consider in initial searches.
• metric::Metric = Cityblock(): The metric to use for nearest neighbour and radius searches.
• kraskov_noise_level::AbstractFloat = 1e-8: Adds a little noise to each value in the embeddings, as suggested by Kraskov
• transform_to_uniform::Bool = false: Independently transform each dimension of the embeddings to be uniformly distributed.
• num_average_samples::Int = -1: Number of target events over which to estimate the log densities, the average of which will be the estimate of the TE. Setting this to -1 will result in all events being used. Note that num_target_events target events will still be used in the estimation of the densities.
• num_surrogates::Integer = 100: The number of surrogate processes to generate (and estimate the TE on) when finding a $p$ value.
• surrogate_num_samples_ratio::AbstractFloat = 1.0: Controls the number of samples used to to construct the alternate set of history embeddings used by our local permutation scheme. This number of samples will be surrogate_num_samples_ratio * num_target_events. Corresponds to $N_{U, \textrm{surrogate}}/N_X$ in ^[1].
• k_perm::Integer = 5: The number of neighbouring source embeddings from which to randomly select a replacement embedding in the local permutation scheme.
• use_exclusion_windows::Bool = true: Whether or not to use dynamic exclusion windows. An explanation of these windows can be found in the Methods subsection `Handling dynamic correlations' in ^[1].
• add_dummy_exclusion_windows::Bool = false If set to true, will add extra dynamic exclusion windows to the non-surrogate embeddings, simulating the windows they would have gotten if they were surrogates.

^[1] Shorten, D. P., Spinney, R. E., Lizier, J.T. (2020). Estimating Transfer Entropy in Continuous Time Between Neural Spike Trains or Other Event-Based Data. bioRxiv 2020.06.16.154377.

^[2] Spinney, R. E., Prokopenko, M., & Lizier, J. T. (2017). Transfer entropy in continuous time, with applications to jump and neural spiking processes. Physical Review E, 95(3), 032319.