IntrinsicTimescales

A Julia package for estimation of timescales from time series data.

Features

• Standard techniques for INT calculation
• Approximate Bayesian Computation (ABC) for parameter inference
• ADVI for variational inference
• Multiple model types:
  • Single timescale
  • Single timescale with oscillations
  • Models supporting missing data
• Summary statistics using periodogram, Welch (from DSP.jl) and Lomb-Scargle (from LombScargle.jl):
  • Autocorrelation function (ACF)
  • Power spectral density (PSD)

Submodules

• Models: Abstract model types and interfaces
• ABC: Approximate Bayesian Computation algorithms
• TuringBackend: Turing.jl integration for ADVI
• SummaryStats: ACF and PSD implementations
• Distances: Distance metrics for ABC
• Utils: Utility functions for analysis
• OrnsteinUhlenbeck: OU process generation using DifferentialEquations.jl
• OneTimescale: Single timescale model
• OneTimescaleAndOsc: Single timescale with oscillations
• OneTimescaleWithMissing: Single timescale with missing data
• OneTimescaleAndOscWithMissing: Single timescale and oscillations with missing data
• Plotting: Plotting functions for results

AbstractTimescaleModel

Abstract type representing models for timescale inference. All concrete model implementations should subtype this.

BaseModel <: AbstractTimescaleModel

Base model structure for timescale inference using various methods.

Fields

• data: Input time series data
• time: Time points corresponding to the data
• data_sum_stats: Pre-computed summary statistics of the data
• fitmethod::Symbol: Fitting method to use. Options: :abc, :advi, :acw
• summary_method::Symbol: Summary statistic type. Options: :psd (power spectral density) or :acf (autocorrelation)
• lags_freqs::AbstractVector{<:Real}: Lags (for ACF) or frequencies (for PSD) at which to compute summary statistics
• prior: Prior distributions for parameters. Can be Vector{Distribution}, single Distribution, or "informed_prior"
• acwtypes::Union{Vector{Symbol}, Symbol}: ACW analysis types (e.g., :ACW50, :ACW0, :ACWe, :tau, :knee)
• distance_method::Symbol: Distance metric type. Options: :linear or :logarithmic
• dt::Real: Time step between observations
• T::Real: Total time span of the data
• numTrials::Real: Number of trials/iterations
• data_mean::Real: Mean of the input data
• data_sd::Real: Standard deviation of the input data

check_acwtypes(acwtypes, possible_acwtypes)

Validate the ACW analysis types against allowed options.

Arguments

• acwtypes: Symbol or Vector of Symbols specifying ACW analysis types
• possible_acwtypes: Vector of allowed ACW analysis types

Returns

• Validated vector of ACW types

Throws

• ErrorException: If invalid ACW types are provided

check_inputs(fitmethod, summary_method)

Validate the fitting method and summary statistic choices.

Arguments

• fitmethod: Symbol specifying the fitting method
• summary_method: Symbol specifying the summary statistic type

Throws

• ArgumentError: If invalid options are provided

check_model_inputs(data, time, fit_method, summary_method, prior, acwtypes, distance_method)

Validate inputs for timescale model construction.

Arguments

• data: Input time series data
• time: Time points corresponding to the data
• fit_method: Fitting method (:abc, :advi)
• summary_method: Summary statistic type (:psd or :acf)
• prior: Prior distribution(s) for parameters
• acwtypes: Types of ACW analysis
• distance_method: Distance metric type (:linear or :logarithmic)

Throws

• ArgumentError: If any inputs are invalid or incompatible

distance_function(model::AbstractTimescaleModel, summary_stats, summary_stats_synth)

Compute distance between two sets of summary statistics.

Arguments

• model: Model instance
• summary_stats: First set of summary statistics
• summary_stats_synth: Second set of summary statistics (typically from synthetic data)

Returns

• Distance value according to model.distance_method

draw_theta(model::AbstractTimescaleModel)

Draw parameter values from the model's prior distributions.

Returns

• Array of proposed model parameters sampled from their respective priors

generate_data(model::AbstractTimescaleModel, theta)

Generate synthetic data using the forward model with given parameters.

Arguments

• model: Model instance
• theta: Array of model parameters

Returns

• Synthetic dataset with same structure as the original data

generate_data_and_reduce(model::AbstractTimescaleModel, theta)

Combined function to generate synthetic data and compute distance from observed data. This is a convenience function commonly used in ABC algorithms.

Arguments

• model: Model instance
• theta: Array of model parameters

Returns

• Distance value between synthetic and observed summary statistics

fit(model::AbstractTimescaleModel, param_dict=nothing)

Fit the timescale model using the specified fitting method.

Arguments

• model: The timescale model instance to fit
• param_dict: Optional dictionary of fitting parameters. If not provided, default parameters will be used.

Returns

For ADVI fitting method:

• samples: Array of posterior samples
• map_estimate: Maximum a posteriori estimate of parameters
• vi_result: Full variational inference result object

For ABC fitting method:

• samples: Array of accepted parameter samples
• weights: Importance weights for the samples
• distances: Distances between simulated and observed summary statistics

summary_stats(model::AbstractTimescaleModel, data)

Compute summary statistics (PSD or ACF) from the data.

Arguments

• model: Model instance
• data: Input data (original or synthetic)

Returns

• Array of summary statistics computed according to model.summary_method

ABCResults

Container for ABC results to standardize plotting interface.

Fields

• theta_history::Vector{Matrix{Float64}}: History of parameter values across iterations
• epsilon_history::Vector{Float64}: History of epsilon values
• acc_rate_history::Vector{Float64}: History of acceptance rates
• weights_history::Vector{Vector{Float64}}: History of weights
• final_theta::Matrix{Float64}: Final accepted parameter values
• final_weights::Vector{Float64}: Final weights
• MAP::Vector{Float64}: Maximum A Posteriori (MAP) estimate of parameters

abc_results(output_record::Vector{NamedTuple})

Construct an ABCResults struct from PMC-ABC output records.

Arguments

• output_record::Vector{NamedTuple}: Vector of named tuples containing PMC-ABC iteration results. Each tuple must contain:
  • theta_accepted: Accepted parameter values
  • epsilon: Epsilon threshold value
  • n_accepted: Number of accepted samples
  • n_total: Total number of samples
  • weights: Importance weights

Returns

• ABCResults: Struct to contain ABC results. See the documentation for ABCResults for more details.

basic_abc(model::Models.AbstractTimescaleModel; kwargs...)

Perform basic ABC rejection sampling. The algorithm stops either when max_iter is reached or when min_accepted samples have been accepted.

Arguments

• model::Models.AbstractTimescaleModel: Model to perform inference on
• epsilon::Float64: Acceptance threshold for distance between simulated and observed data
• max_iter::Integer: Maximum number of iterations to perform
• min_accepted::Integer: Minimum number of accepted samples required before stopping
• pmc_mode::Bool=false: Whether to use PMC proposal distribution instead of prior
• weights::Array{Float64}: Importance weights for PMC sampling (only used if pmc_mode=true)
• theta_prev::Array{Float64}: Previous parameters for PMC sampling (only used if pmc_mode=true)
• tau_squared::Array{Float64}: Covariance matrix for PMC sampling (only used if pmc_mode=true)
• show_progress::Bool=true: Whether to show progress bar with acceptance count and speed

Returns

NamedTuple containing:

• samples::Matrix{Float64}: Matrix (maxiter × nparams) of all proposed parameters
• isaccepted::Vector{Bool}: Boolean mask of accepted samples for first n_total iterations
• theta_accepted::Matrix{Float64}: Matrix (naccepted × nparams) of accepted parameters
• distances::Vector{Float64}: Vector of distances for first n_total iterations
• n_accepted::Int: Number of accepted samples
• n_total::Int: Total number of iterations performed
• epsilon::Float64: Acceptance threshold used
• weights::Vector{Float64}: Uniform weights (ones) for accepted samples
• tau_squared::Matrix{Float64}: Zero matrix (nparams × nparams) for basic ABC
• eff_sample::Int: Effective sample size (equals n_accepted in basic ABC)

Implementation Details

1. Draws parameters either from prior (basic mode) or PMC proposal (pmc_mode)
2. Generates synthetic data and computes distance to observed data
3. Accepts parameters if distance ≤ epsilon
4. Stops when either maxiter reached or minaccepted samples accepted
5. Returns uniform weights and zero covariance matrix in basic mode

calc_weights(theta_prev, theta, tau_squared, weights, prior)

Calculate importance weights for PMC-ABC algorithm.

Arguments

• theta_prev::Union{Vector{Float64}, Matrix{Float64}}: Previously accepted parameters. For multiple parameters, each row is a sample and each column is a parameter
• theta::Union{Vector{Float64}, Matrix{Float64}}: Current parameters in same format as theta_prev
• tau_squared::Matrix{Float64}: Covariance matrix for the proposal distribution
• weights::Vector{Float64}: Previous iteration's importance weights
• prior::Union{Vector, dist.Distribution}: Prior distribution(s). Can be single distribution or vector of distributions

Returns

• Vector{Float64}: Normalized importance weights (sum to 1)

compute_adaptive_alpha(iteration::Integer, current_acc_rate::Float64, target_acc_rate::Float64; kwargs...)

Compute an adaptive step size (alpha) for epsilon adjustment in ABC, based on iteration progress and distance from target acceptance rate.

Arguments

Required Arguments

• iteration::Integer: Current iteration number
• current_acc_rate::Float64: Current acceptance rate
• target_acc_rate::Float64: Target acceptance rate to achieve

Optional Keyword Arguments

Bounds Parameters

• alpha_max::Float64=0.9: Maximum allowed alpha value
• alpha_min::Float64=0.1: Minimum allowed alpha value
• total_iterations::Integer=100: Total number of iterations planned

Adaptation Parameters

• acc_rate_far::Float64=2.0: Relative difference threshold for "far from target"
• acc_rate_close::Float64=0.2: Relative difference threshold for "close to target"
• alpha_far_mult::Float64=1.5: Multiplier for alpha when far from target
• alpha_close_mult::Float64=0.5: Multiplier for alpha when close to target

Returns

• Float64: Adaptive alpha value between alpha_min and alpha_max

Implementation Details

1. Computes base alpha using linear decay between max and min:

   • base_alpha = alpha_max * (1 - progress) + alpha_min * progress

   where progress = iteration/total_iterations

2. Adjusts base alpha based on relative difference from target:

   • acc_rate_diff = |current_acc_rate - target_acc_rate|/target_acc_rate

3. Final alpha selection:

   • If acc_rate_diff > acc_rate_far: More aggressive adaptation alpha = min(alpha_max, base_alpha * alpha_far_mult)
   • If acc_rate_diff < acc_rate_close: More conservative adaptation alpha = max(alpha_min, base_alpha * alpha_close_mult)
   • Otherwise: Use base alpha alpha = base_alpha

draw_theta_pmc(model, theta_prev, weights, tau_squared; jitter::Float64=1e-5)

Draw new parameter values using the PMC proposal distribution.

Arguments

• model: Model instance
• theta_prev: Previously accepted parameters
• weights: Importance weights from previous iteration
• tau_squared: Covariance matrix for proposal distribution
• jitter::Float64=1e-5: Small value added to covariance diagonal for numerical stability

Returns

• Vector of proposed parameters

effective_sample_size(w::Vector{Float64})

Calculate effective sample size (ESS) from importance weights.

Arguments

• w::Vector{Float64}: Vector of importance sampling weights (need not be normalized)

Returns

• Float64: Effective sample size

Details

The effective sample size is always less than or equal to the actual number of samples. It reaches its maximum (equal to sample size) when all weights are equal, and approaches its minimum (1) when one weight dominates all others.

find_MAP(theta_accepted::AbstractArray{Float64}, N::Integer=10000)

Find Maximum A Posteriori (MAP) estimates with grid search.

Arguments

• theta_accepted::AbstractArray{Float64}: Matrix of accepted samples from the final step of ABC
• N::Integer=10000: Number of random grid points to evaluate for each parameter

Returns

• theta_map::Vector{Float64}: MAP estimate of the parameters

get_param_dict_abc()

Get default parameter dictionary for ABC algorithm.

Returns

Dictionary containing default values for all ABC parameters including:

Basic ABC Parameters

• :epsilon_0 => 1.0: Initial epsilon threshold
• :max_iter => 10000: Maximum iterations per step
• :min_accepted => 100: Minimum number of accepted samples
• :steps => 30: Maximum PMC steps
• :sample_only => false: If true, only perform sampling without adaptation

Acceptance Rate Parameters

• :minAccRate => 0.01: Minimum acceptance rate before early stopping
• :target_acc_rate => 0.01: Target acceptance rate
• :target_epsilon => 1e-4: Target epsilon for early stopping

Display Parameters

• :show_progress => true: Show progress bar
• :verbose => true: Print detailed progress information

Numerical Stability Parameters

• :jitter => 1e-6: Small value added to covariance matrix

Epsilon Selection Parameters

• :distance_max => 10.0: Maximum valid distance
• :quantile_lower => 25.0: Lower quantile for epsilon bounds
• :quantile_upper => 75.0: Upper quantile for epsilon bounds
• :quantile_init => 50.0: Initial quantile
• :acc_rate_buffer => 0.1: Allowed deviation from target rate

Adaptive Alpha Parameters

• :alpha_max => 0.9: Maximum adaptation rate
• :alpha_min => 0.1: Minimum adaptation rate
• :acc_rate_far => 2.0: Threshold for "far from target"
• :acc_rate_close => 0.2: Threshold for "close to target"
• :alpha_far_mult => 1.5: Multiplier when far from target
• :alpha_close_mult => 0.5: Multiplier when close to target

Early Stopping Parameters

• :convergence_window => 5: Steps to check for convergence
• :theta_rtol => 1e-2: Relative tolerance for convergence
• :theta_atol => 1e-3: Absolute tolerance for convergence

MAP Estimation Parameters

• :N => 10000: Number of grid points for MAP estimation

Example

params = get_param_dict_abc()
params[:epsilon_0] = 0.5  # Modify initial epsilon
params[:max_iter] = 5000  # Reduce maximum iterations

pmc_abc(model::Models.AbstractTimescaleModel; kwargs...)

Perform Population Monte Carlo Approximate Bayesian Computation (PMC-ABC) inference.

Arguments

Basic ABC parameters

• model::Models.AbstractTimescaleModel: Model to perform inference on
• epsilon_0::Real=1.0: Initial epsilon threshold for acceptance
• max_iter::Integer=10000: Maximum number of iterations per step
• min_accepted::Integer=100: Minimum number of accepted samples required
• steps::Integer=10: Maximum number of PMC steps to perform
• sample_only::Bool=false: If true, only perform sampling without adaptation

Acceptance rate parameters

• minAccRate::Float64=0.01: Minimum acceptance rate before early stopping
• target_acc_rate::Float64=0.01: Target acceptance rate for epsilon adaptation
• target_epsilon::Float64=5e-3: Target epsilon value for early stopping

Display parameters

• show_progress::Bool=true: Whether to show progress bar
• verbose::Bool=true: Whether to print detailed progress information

Numerical stability parameters

• jitter::Float64=1e-6: Small value added to covariance matrix for stability

Epsilon selection parameters

• distance_max::Float64=10.0: Maximum distance to consider valid
• quantile_lower::Float64=25.0: Lower quantile for epsilon adjustment
• quantile_upper::Float64=75.0: Upper quantile for epsilon adjustment
• quantile_init::Float64=50.0: Initial quantile when no acceptance rate
• acc_rate_buffer::Float64=0.1: Buffer around target acceptance rate

Adaptive alpha parameters

• alpha_max::Float64=0.9: Maximum adaptation rate
• alpha_min::Float64=0.1: Minimum adaptation rate
• acc_rate_far::Float64=2.0: Threshold for "far from target" adjustment
• acc_rate_close::Float64=0.2: Threshold for "close to target" adjustment
• alpha_far_mult::Float64=1.5: Multiplier for alpha when far from target
• alpha_close_mult::Float64=0.5: Multiplier for alpha when close to target

Early stopping parameters

• convergence_window::Integer=3: Number of steps to check for convergence
• theta_rtol::Float64=1e-2: Relative tolerance for parameter convergence
• theta_atol::Float64=1e-3: Absolute tolerance for parameter convergence

Returns

ABCResults: A struct containing:

• ABCResults.theta_history: Parameter value history across iterations
• ABCResults.epsilon_history: Epsilon value history
• ABCResults.acc_rate_history: Acceptance rate history
• ABCResults.weight_history: Weight history
• ABCResults.theta_final: Final parameter values
• ABCResults.weights_final: Final weights
• ABCResults.theta_map: MAP estimate

Early Stopping Conditions

The algorithm stops and returns results if any of these conditions are met:

1. Acceptance rate falls below minAccRate
2. Parameters converge within tolerances over convergence_window steps
3. Epsilon falls below target_epsilon
4. Maximum number of steps reached

Implementation Details

1. First step uses basic ABC with prior sampling
2. Subsequent steps use PMC proposal with adaptive epsilon
3. Epsilon is adjusted based on acceptance rates and distance quantiles
4. Covariance and weights are updated each step unless sample_only=true
5. Parameter convergence is checked using both relative and absolute tolerances

select_epsilon(distances::Vector{Float64}, current_epsilon::Float64; kwargs...)

Adaptively select the epsilon threshold for ABC based on acceptance rates and distance distribution. Uses a combination of quantile-based bounds and adaptive step sizes to adjust epsilon towards achieving the target acceptance rate.

Arguments

Required Arguments

• distances::Vector{Float64}: Vector of distances from ABC simulations
• current_epsilon::Float64: Current epsilon threshold value

Optional Keyword Arguments

Acceptance Rate Parameters

• target_acc_rate::Float64=0.01: Target acceptance rate to achieve
• current_acc_rate::Float64=0.0: Current acceptance rate
• acc_rate_buffer::Float64=0.1: Allowed deviation from target acceptance rate

Iteration Parameters

• iteration::Integer=1: Current iteration number
• total_iterations::Integer=100: Total number of iterations planned

Distance Processing Parameters

• distance_max::Float64=10.0: Maximum valid distance (larger values filtered out)
• quantile_lower::Float64=25.0: Lower quantile for epsilon bounds
• quantile_upper::Float64=75.0: Upper quantile for epsilon bounds
• quantile_init::Float64=50.0: Initial quantile for first iteration

Adaptive Step Size Parameters

• alpha_max::Float64=0.9: Maximum adaptation rate
• alpha_min::Float64=0.1: Minimum adaptation rate
• acc_rate_far::Float64=2.0: Threshold for "far from target" adjustment
• acc_rate_close::Float64=0.2: Threshold for "close to target" adjustment
• alpha_far_mult::Float64=1.5: Multiplier for alpha when far from target
• alpha_close_mult::Float64=0.5: Multiplier for alpha when close to target

Returns

• Float64: New epsilon value

Implementation Details

1. Filters out NaN and distances larger than distance_max
2. Computes quantile-based bounds for epsilon adjustment
3. Uses adaptive alpha value based on iteration and acceptance rate (see compute_adaptive_alpha)
4. For first iteration (iteration=1):
   • Returns initial quantile of valid distances
5. For subsequent iterations:
   • If acceptance rate too high: decreases epsilon by (1-alpha)
   • If acceptance rate too low: increases epsilon by (1+alpha)
   • Keeps epsilon unchanged if within buffer of target rate
6. Always constrains new epsilon between quantile bounds

Notes

• Returns current epsilon if no valid distances are found
• Uses computeadaptivealpha for step size calculation
• Adjustments are proportional to distance from target acceptance rate

weighted_covar(x::Matrix{Float64}, w::Vector{Float64})

Calculate weighted covariance matrix.

Arguments

• x::Matrix{Float64}: Matrix of observations where each row is an observation and each column is a variable
• w::Vector{Float64}: Vector of weights corresponding to each observation (row of x)

Returns

• Matrix{Float64}: Weighted covariance matrix of size (nvariables × nvariables)

ACW

Module providing autocorrelation width (ACW) calculations for time series analysis, including:

• ACW-0 (zero-crossing)
• ACW-50 (50% decay)
• ACW-euler (1/e decay)
• Exponential decay timescale (tau)
• Knee frequency estimation

ACWResults

Structure holding ACW analysis inputs and results.

Fields

• fs::Real: Sampling frequency
• acw_results: Computed ACW values (type depends on number of ACW types requested)
• acwtypes::Union{Vector{<:Symbol}, Symbol}: Types of ACW computed
• n_lags::Union{Int, Nothing}: Number of lags used for ACF calculation
• freqlims::Union{Tuple{Real, Real}, Nothing}: Frequency limits used for spectral analysis
• acf::Union{AbstractVector{<:Real}, AbstractArray{<:Real}, Nothing}: Autocorrelation function
• psd::Union{AbstractVector{<:Real}, AbstractArray{<:Real}, Nothing}: Power spectral density
• freqs::Union{AbstractVector{<:Real}, AbstractArray{<:Real}, Nothing}: Frequency vector for PSD
• lags::Union{AbstractVector{<:Real}, AbstractArray{<:Real}, Nothing}: Lag vector for ACF
• x_dim::Union{Int, Nothing}: Dimension index corresponding to x-axis (lags/freqs)

Notes

• Supported ACW types: :acw0, :acw50, :acweuler, :auc, :tau, :knee
• Results order matches input acwtypes order
• If only one ACW type is requested, acw_results is a scalar; otherwise it's a vector

acw(data, fs; acwtypes=possible_acwtypes, n_lags=nothing, freqlims=nothing, time=nothing, 
    dims=ndims(data), return_acf=true, return_psd=true, average_over_trials=false,
    trial_dims::Int=setdiff([1, 2], dims)[1], skip_zero_lag::Bool=false, max_peaks::Int=1, oscillation_peak::Bool=true,
    allow_variable_exponent::Bool=false, parallel::Bool=false)

Compute various timescale measures for time series data. For detailed documentaion, see https://duodenum96.github.io/IntrinsicTimescales.jl/stable/acw/.

Arguments

• data::AbstractArray{<:Real}: Input time series data
• fs::Real: Sampling frequency
• acwtypes::Union{Vector{Symbol}, Symbol}=[:acw0, :acw50, :acweuler, :auc, :tau, :knee]: Types of ACW to compute.

Supported ACW types:

• :acw0 - Time to first zero crossing
• :acw50 - Time to 50% decay
• :acweuler - Time to 1/e decay
• :auc - Area under curve of ACF before ACW0
• :tau - Exponential decay timescale
• :knee - Knee frequency from spectral analysis
• n_lags::Union{Int, Nothing}=nothing: Number of lags for ACF calculation. If not specified, uses 1.1 * ACW0.
• freqlims::Union{Tuple{Real, Real}, Nothing}=nothing: Frequency limits for spectral analysis. If not specified, uses full frequency range.
• time::Union{Vector{Real}, Nothing}=nothing: Time vector. This is required for Lomb-Scargle method in the case of missing data.
• dims::Int=ndims(data): Dimension along which to compute ACW (Dimension of time)
• return_acf::Bool=true: Whether to return the ACF
• return_psd::Bool=true: Whether to return the PSD
• average_over_trials::Bool=false: Whether to average the ACF or PSD over trials
• trial_dims::Int=setdiff([1, 2], dims)[1]: Dimension along which to average the ACF or PSD over trials (Dimension of trials)
• skip_zero_lag::Bool=false: Whether to skip the zero lag for fitting an exponential decay function. Used only for :tau.
• max_peaks::Int=1: Maximum number of oscillatory peaks to fit in spectral analysis
• oscillation_peak::Bool=true: Whether to fit an oscillation peak in the spectral analysis
• allow_variable_exponent::Bool=false: Whether to allow variable exponent in spectral fitting
• parallel::Bool=false: Whether to use parallel computation

Returns

• ACWResults: Structure containing computed ACW measures and intermediate results

Fields of the ACWResults structure:

• fs::Real: Sampling frequency
• acw_results: Computed ACW values (type depends on number of ACW types requested)
• acwtypes::Union{Vector{<:Symbol}, Symbol}: Types of ACW computed
• n_lags::Union{Int, Nothing}: Number of lags used for ACF calculation
• freqlims::Union{Tuple{Real, Real}, Nothing}: Frequency limits used for spectral analysis
• acf::Union{AbstractVector{<:Real}, AbstractArray{<:Real}, Nothing}: Autocorrelation function
• psd::Union{AbstractVector{<:Real}, AbstractArray{<:Real}, Nothing}: Power spectral density
• freqs::Union{AbstractVector{<:Real}, AbstractArray{<:Real}, Nothing}: Frequency vector for PSD
• lags::Union{AbstractVector{<:Real}, AbstractArray{<:Real}, Nothing}: Lag vector for ACF
• x_dim::Union{Int, Nothing}: Dimension index corresponding to x-axis (lags/freqs)

ADVIResults{T<:Real}

Container for ADVI (Automatic Differentiation Variational Inference) results.

Fields

• samples::AbstractArray{T}: Matrix of posterior samples
• MAP::AbstractVector{T}: Maximum a posteriori estimates
• variational_posterior::Any: Turing variational posterior object containing full inference results

create_turing_model(model, data_sum_stats; σ_prior=Exponential(1))

Create a Turing probabilistic model for variational inference.

Arguments

• model: Model instance containing prior distributions and data generation methods
• data_sum_stats: Summary statistics of the observed data
• σ_prior=Exponential(1): Prior distribution for the uncertainty parameter σ

Returns

• Turing model object ready for inference

Notes

The created model includes:

• Parameter sampling from truncated priors (positive values only)
• Data generation using the model's forward simulation
• Likelihood computation using Normal distribution

fit_vi(model; n_samples=4000, n_iterations=10, n_elbo_samples=20, 
       optimizer=AutoForwardDiff())

Perform variational inference using ADVI (Automatic Differentiation Variational Inference).

Arguments

• model: Model instance to perform inference on
• n_samples::Int=4000: Number of posterior samples to draw
• n_iterations::Int=10: Number of ADVI iterations
• n_elbo_samples::Int=20: Number of samples for ELBO estimation
• optimizer=AutoForwardDiff(): Optimization algorithm for ADVI

Returns

• ADVIResults: Container with inference results including:
  • samples_matrix: Matrix of posterior samples
  • MAP: Maximum a posteriori parameter estimates
  • variational_posterior: Turing variational posterior object containing full inference results

Notes

Uses Turing.jl's ADVI implementation for fast approximate Bayesian inference. The model is automatically constructed with appropriate priors and likelihood.

get_param_dict_advi()

Get default parameter dictionary for ADVI (Automatic Differentiation Variational Inference) algorithm.

Returns

Dictionary containing default values for ADVI parameters including:

• n_samples: Number of posterior samples to draw (default: 4000)
• n_iterations: Number of ADVI iterations (default: 50)
• n_elbo_samples: Number of samples for ELBO estimation (default: 20)
• autodiff: Automatic differentiation backend (default: AutoForwardDiff())

SummaryStats

Module for computing various summary statistics from time series data. Includes functions for:

• Autocorrelation (FFT and time-domain methods)
• Power spectral density (periodogram and Welch methods)
• Cross-correlation
• Special handling for missing data (NaN values)

_comp_psd_lombscargle(times, data, frequency_grid)

Internal function to compute Lomb-Scargle periodogram for a single time series.

Arguments

• times: Time points vector (without NaN)
• data: Time series data (without NaN)
• frequency_grid: Pre-computed frequency grid

Returns

• power: Lomb-Scargle periodogram values

Notes

• Uses LombScargle.jl for core computation
• Assumes data has been pre-processed and doesn't contain NaN values

acf_statsmodels(x::Vector{T}; kwargs...) where {T <: Real}

Julia implementation of statsmodels.tsa.stattools.acf function. Only for testing.

Arguments

• x: Time series data vector
• adjusted=false: Use n-k denominators if true
• nlags=nothing: Number of lags (default: min(10*log10(n), n-1))
• qstat=false: Return Ljung-Box Q-statistics
• isfft=false: Use FFT method
• alpha=nothing: Confidence level for intervals
• bartlett_confint=false: Use Bartlett's formula
• missing_handling="conservative": NaN handling method

Returns

• Vector of autocorrelation values

Notes

• Supports multiple missing data handling methods:
  • "none": No checks
  • "raise": Error on NaN
  • "conservative": NaN-aware computations
  • "drop": Remove NaN values

Estimate autocovariances. Translated to Julia from statsmodels.tsa.stattools

Parameters

x : array_like Time series data. Must be 1d. adjusted : bool, default False If True, then denominators is n-k, otherwise n. demean : bool, default True If True, then subtract the mean x from each element of x. fft : bool, default True If True, use FFT convolution. This method should be preferred for long time series. missing : str, default "none" A string in ["none", "raise", "conservative", "drop"] specifying how the NaNs are to be treated. "none" performs no checks. "raise" raises an exception if NaN values are found. "drop" removes the missing observations and then estimates the autocovariances treating the non-missing as contiguous. "conservative" computes the autocovariance using nan-ops so that nans are removed when computing the mean and cross-products that are used to estimate the autocovariance. When using "conservative", n is set to the number of non-missing observations. nlag : {int, None}, default None Limit the number of autocovariances returned. Size of returned array is nlag + 1. Setting nlag when fft is False uses a simple, direct estimator of the autocovariances that only computes the first nlag + 1 values. This can be much faster when the time series is long and only a small number of autocovariances are needed.

Returns

ndarray The estimated autocovariances.

References

.. [1] Parzen, E., 1963. On spectral analysis with missing observations and amplitude modulation. Sankhya: The Indian Journal of Statistics, Series A, pp.383-392.

bat_autocorr(x::AbstractVector{<:Real})

Compute the normalized autocorrelation function of a 1D time series using FFT. Returns a vector containing the autocorrelation values.

bat_autocorr(x::AbstractMatrix{<:Real})

Compute the normalized autocorrelation function of a 2D time series using FFT. data is a matrix with dimensions (nseries × ntimepoints)

returns a matrix with dimensions (nseries × nlags)

bat_integrated_autocorr_len(
    v::AbstractVectorOfSimilarVectors{<:Real};
    c::Integer = 5, tol::Integer = 50, strict = true
)

Estimate the integrated autocorrelation length of variate series v.

• c: Step size for window search.

• tol: Minimum number of autocorrelation times needed to trust the estimate.

• strict: Throw exception if result is not trustworthy

This estimate uses the iterative procedure described on page 16 of Sokal's notes to determine a reasonable window size.

Ported to Julia from the emcee Python package, under MIT License. Original authors Dan Foreman-Mackey et al.

bat_integrated_autocorr_weight(
    samples::DensitySampleVector;
    c::Integer = 5, tol::Integer = 50, strict = true
)

Estimate the integrated autocorrelation weight of samples.

Uses bat_integrated_autocorr_len.

comp_ac_fft(data::AbstractArray{T}; dims::Int=ndims(data), n_lags::Integer=size(data, dims), parallel::Bool=false) where {T <: Real}

Compute autocorrelation using FFT along specified dimension.

Arguments

• data: Array of time series data
• dims: Dimension along which to compute autocorrelation (defaults to last dimension)
• n_lags: Number of lags to compute (defaults to size of data along specified dimension)
• parallel: Whether to use parallel computation

Returns

Array with autocorrelation values, the specified dimension becomes the dimension of lags while the other dimensions denote ACF values

comp_ac_fft(data::Vector{T}; n_lags::Real=length(data)) where {T <: Real}

Compute autocorrelation using FFT method.

Arguments

• data: Input time series vector
• n_lags: Number of lags to compute (defaults to length of data)

Returns

• Vector of autocorrelation values from lag 0 to n_lags-1

Notes

• Uses FFT for efficient computation
• Pads data to next power of 2 for FFT efficiency
• Normalizes by variance (first lag)

comp_ac_time(data::AbstractArray{T}, max_lag::Integer; dims::Int=ndims(data), parallel::Bool=false) where {T <: Real}

Compute autocorrelation in time domain along specified dimension.

Arguments

• data: Array of time series data
• max_lag: Maximum lag to compute
• dims: Dimension along which to compute autocorrelation (defaults to last dimension)
• parallel=false: Whether to use parallel computation

Returns

Array with autocorrelation values, the specified dimension becomes the dimension of lags while the other dimensions denote ACF values

comp_ac_time(data::AbstractArray{T}, max_lag::Integer; dims::Int=ndims(data)) where {T <: Real}

Compute autocorrelation in time domain along specified dimension.

Arguments

• data: Array of time series data
• max_lag: Maximum lag to compute
• dims: Dimension along which to compute autocorrelation (defaults to last dimension)

Returns

Array with autocorrelation values, the specified dimension becomes the dimension of lags while the other dimensions denote ACF values

comp_ac_time_missing(data::AbstractArray{T}; kwargs...) where {T <: Real}

Compute autocorrelation for data with missing values.

Arguments

• data: Time series data (may contain NaN)
• dims=ndims(data): Dimension along which to compute
• n_lags=size(data,dims): Number of lags to compute
• parallel=false: Whether to use parallel computation

Returns

• Array of autocorrelation values

Notes

• Handles missing data using missing="conservative" approach of

statsmodels.tsa.stattools.acf. See https://www.statsmodels.org/dev/generated/statsmodels.tsa.stattools.acf.html for details.

comp_cc(data1::AbstractArray{T}, data2::AbstractArray{T}, max_lag::Integer;
       dims::Int=ndims(data1)) where {T <: Real}

Compute cross-correlation between two arrays along specified dimension.

Arguments

• data1: First array of time series data
• data2: Second array of time series data
• max_lag: Maximum lag to compute
• dims: Dimension along which to compute cross-correlation (defaults to last dimension)

Returns

Array with cross-correlation values, reduced along specified dimension

comp_psd(x::AbstractArray{T}, fs::Real; kwargs...) where {T <: Real}

Compute power spectral density using periodogram or welch method.

Arguments

• x: Time series data (time × channels)
• fs: Sampling frequency
• dims=ndims(x): Dimension along which to compute PSD
• method="periodogram": Method to use ("periodogram" or "welch")
• window=dsp.hamming: Window function
• n=div(size(x,dims),8): Window size for Welch method
• noverlap=div(n,2): Overlap for Welch method
• parallel=false: Whether to use parallel computation

Returns

• power: Power spectral density values (excludes DC component)
• freqs: Corresponding frequencies (excludes DC component)

comp_psd_adfriendly(x::AbstractArray{<:Real}, fs::Real; dims::Int=ndims(x), parallel::Bool=false)

Compute power spectral density using an automatic differentiation (AD) friendly implementation.

Arguments

• x: Time series data
• fs: Sampling frequency
• dims=ndims(x): Dimension along which to compute PSD
• parallel=false: Whether to use parallel computation

Returns

• power: Power spectral density values
• freqs: Corresponding frequencies

comp_psd_lombscargle(times, data, nanmask, dt; dims=ndims(data))

Compute Lomb-Scargle periodogram for data with missing values.

Arguments

• times: Time points vector
• data: Time series data (may contain NaN)
• nanmask: Boolean mask indicating NaN positions
• dt: Time step
• dims=ndims(data): Dimension along which to compute

Returns

• power: Lomb-Scargle periodogram values
• frequency_grid: Corresponding frequencies

Notes

• Handles irregular sampling due to missing data
• Uses frequency grid based on shortest valid time series
• Automatically determines appropriate frequency range

prepare_lombscargle(times, data, nanmask)

Prepare data for Lomb-Scargle periodogram computation by handling missing values.

Arguments

• times: Time points vector
• data: Time series data (may contain NaN)
• nanmask: Boolean mask indicating NaN positions

Returns

• valid_times: Time points with NaN values removed
• valid_data: Data points with NaN values removed
• frequency_grid: Suggested frequency grid for analysis

Distances

Module providing distance metrics for comparing summary statistics in ABC inference. Currently implements linear (L2) and logarithmic distances.

linear_distance(data, synth_data)

Compute mean squared error (MSE) between summary statistics.

Arguments

• data::Union{AbstractArray, Real}: Observed data summary statistics
• synth_data::Union{AbstractArray, Real}: Simulated data summary statistics

Returns

• Float64: Mean squared difference between data and synth_data

Notes

• Handles both scalar and array inputs
• For arrays, computes element-wise differences before averaging
• Useful for comparing summary statistics on linear scales

logarithmic_distance(data, synth_data)

Compute mean squared distance between logarithms of summary statistics.

Arguments

• data::Union{AbstractArray, Real}: Observed data summary statistics
• synth_data::Union{AbstractArray, Real}: Simulated data summary statistics

Returns

• Float64: Mean squared difference between log(data) and log(synth_data)

Notes

• Handles both scalar and array inputs
• For arrays, computes element-wise log differences before averaging
• Useful for comparing summary statistics spanning multiple orders of magnitude
• Assumes all values are positive

Utils

Module providing utility functions for time series analysis, including:

• Exponential decay fitting
• Oscillation peak detection
• Knee frequency estimation
• Lorentzian fitting
• ACF width calculations

acw0(lags, acf; dims=ndims(acf), parallel=false)

Compute the ACW0 (autocorrelation width at zero crossing) along specified dimension.

Arguments

• lags::AbstractVector{T}: Vector of lag values
• acf::AbstractArray{T}: Array of autocorrelation values
• dims::Int=ndims(acf): Dimension along which to compute ACW0
• parallel=false: Whether to use parallel computation

Returns

• First lag where autocorrelation crosses zero

Notes

• Alternative measure of characteristic timescale
• More sensitive to noise than ACW50

acw50(lags, acf; dims=ndims(acf), parallel=false)

Compute the ACW50 (autocorrelation width at 50%) along specified dimension.

Arguments

• lags::AbstractVector{T}: Vector of lag values
• acf::AbstractArray{T}: Array of autocorrelation values
• dims::Int=ndims(acf): Dimension along which to compute ACW50
• parallel=false: Whether to use parallel computation

Returns

• First lag where autocorrelation falls below 0.5

Notes

• Used for estimating characteristic timescales
• Related to tau by: tau = -acw50/log(0.5)

acw_romberg(dt, acf; dims=ndims(acf), parallel=false)

Calculate the area under the curve of ACF using Romberg integration.

Arguments

• dt::Real: Time step
• acf::AbstractVector: Array of autocorrelation values
• dims::Int=ndims(acf): Dimension along which to compute
• parallel=false: Whether to use parallel computation

Returns

• AUC of ACF

Notes

• Returns only the integral value, discarding the error estimate

acweuler(lags, acf; dims=ndims(acf), parallel=false)

Compute the ACW at 1/e (≈ 0.368) along specified dimension.

Arguments

• lags::AbstractVector{T}: Vector of lag values
• acf::AbstractVector{S}: Array of autocorrelation values
• dims::Int=ndims(acf): Dimension along which to compute
• parallel=false: Whether to use parallel computation

Returns

• First lag where autocorrelation falls below 1/e

Notes

• For exponential decay, equals the timescale parameter tau

expdecay(tau, lags)

Compute exponential decay function.

Arguments

• tau::Real: Timescale parameter
• lags::AbstractVector: Time lags

Returns

• Vector of exp(-t/tau) values

Notes

• Used for fitting autocorrelation functions
• Assumes exponential decay model: acf = exp(-t/tau)

expdecay_3_parameters(p, lags)

Compute exponential decay function with amplitude and offset parameters. This is used in acw with the setting skip_zero_lag=true.

Arguments

• p::AbstractVector: Parameters [A, tau, B] where:
  • A: Amplitude parameter
  • tau: Timescale parameter
  • B: Offset parameter
• lags::AbstractVector: Time lags

Returns

• Vector of A*(exp(-t/tau) + B) values

find_knee_frequency(psd, freqs; dims=ndims(psd), min_freq=freqs[1], max_freq=freqs[end], constrained=false, allow_variable_exponent=false, parallel=false)

Find knee frequency by fitting Lorentzian to power spectral density.

Arguments

• psd::AbstractArray{T}: Power spectral density values
• freqs::Vector{T}: Frequency values
• dims::Int=ndims(psd): Dimension along which to compute
• min_freq::T=freqs[1]: Minimum frequency to consider
• max_freq::T=freqs[end]: Maximum frequency to consider
• constrained::Bool=false: Whether to use constrained optimization
• allow_variable_exponent::Bool=false: Whether to allow variable exponent (PLE)
• parallel=false: Whether to use parallel computation

Returns

• Vector of the fit for the equation amp/(1 + (f/knee)^{exponent}).

If allowvariableexponent=false, assumes exponent=2 and returns [amplitude, kneefrequency]. If true, returns [amplitude, kneefrequency, exponent].

Notes

• Uses Lorentzian fitting with NonlinearSolve.jl or Optimization.jl
• Initial guess for amplitude based on low frequency power
• Initial guess for knee at half-power point
• Initial guess for exponent is 2.0 when variable exponent allowed

find_oscillation_peak(psd, freqs; min_freq=5.0/1000.0, max_freq=50.0/1000.0, min_prominence_ratio=0.1)

Find dominant oscillatory peak in power spectral density.

Arguments

• psd::AbstractVector: Power spectral density values
• freqs::AbstractVector: Frequency values
• min_freq::Real=5.0/1000.0: Minimum frequency to consider
• max_freq::Real=50.0/1000.0: Maximum frequency to consider
• min_prominence_ratio::Real=0.1: Minimum peak prominence as fraction of max PSD

Returns

• Frequency of most prominent peak, or NaN if no significant peak found

Notes

• Uses peak prominence for robustness
• Filters peaks by minimum prominence threshold
• Returns NaN if no peaks meet criteria

fit_expdecay(lags, acf; dims=ndims(acf), parallel=false)

Fit exponential decay to autocorrelation function.

Arguments

• lags::AbstractVector{T}: Time lags
• acf::AbstractArray{T}: Autocorrelation values
• dims::Int=ndims(acf): Dimension along which to fit
• parallel=false: Whether to use parallel computation

Returns

• Fitted timescale parameter(s)

Notes

• Uses NonlinearSolve.jl with FastShortcutNLLSPolyalg
• Initial guess based on ACW50

fit_expdecay_3_parameters(lags, acf; parallel=false)

Fit a 3-parameter exponential decay function to autocorrelation data ( A*(exp(-t/tau) + B) ). Excludes lag 0 from fitting.

Arguments

• lags::AbstractVector{T}: Time lags
• acf::AbstractVector{T}: Autocorrelation values
• parallel=false: Whether to use parallel computation

Returns

• Fitted timescale parameter (tau)

Notes

• Initial guess: A=0.5, tau from ACW50, B=0.0

fit_gaussian(psd, freqs, initial_peak; min_freq=freqs[1], max_freq=freqs[end])

Fit Gaussian to power spectral density around a peak.

Arguments

• psd::AbstractVector{<:Real}: Power spectral density values
• freqs::AbstractVector{<:Real}: Frequency values
• initial_peak::Real: Initial guess for center frequency
• min_freq::Real: Minimum frequency to consider
• max_freq::Real: Maximum frequency to consider

Returns

• Vector{Float64}: Fitted parameters [amplitude, centerfreq, stddev]

Notes

• Uses initial peak location from findoscillationpeak

fooof_fit(psd, freqs; dims=ndims(psd), min_freq=freqs[1], max_freq=freqs[end], 
          oscillation_peak=true, max_peaks=3, allow_variable_exponent=false, constrained=false, parallel=false)

Perform FOOOF-style fitting of power spectral density. The default behavior is to fit a Lorentzian with PLE = 2. If allowvariableexponent=true, the function will fit a Lorentzian with variable PLE.

Arguments

• psd::AbstractArray{T}: Power spectral density values
• freqs::Vector{T}: Frequency values
• dims::Int=ndims(psd): Dimension along which to compute
• min_freq::T=freqs[1]: Minimum frequency to consider
• max_freq::T=freqs[end]: Maximum frequency to consider
• oscillation_peak::Bool=true: Whether to compute oscillation peaks
• max_peaks::Int=3: Maximum number of oscillatory peaks to fit
• return_only_knee::Bool=false: Whether to return only knee frequency
• allow_variable_exponent::Bool=false: Whether to allow variable exponent (PLE)
• constrained::Bool=false: Whether to use constrained optimization
• parallel=false: Whether to use parallel computation

Returns

If returnonlyknee=false:

• Tuple of (kneefrequency, oscillationparameters) where oscillationparameters is Vector of (centerfreq, amplitude, std_dev) for each peak

If returnonlyknee=true:

• knee_frequency only

Notes

• Implements iterative FOOOF-style fitting:
  1. Fit initial Lorentzian to PSD
  2. Find and fit Gaussian peaks iteratively
  3. Subtract all Gaussians from original PSD
  4. Refit Lorentzian to cleaned PSD
• Default behavior fits standard Lorentzian (PLE = 2)
• If allowvariableexponent=true, fits Lorentzian with variable PLE

gaussian(f, u)

Gaussian function for fitting oscillations.

Arguments

• f::AbstractVector: Frequency values
• u::Vector: Parameters [amplitude, centerfreq, stddev]

Returns

• Vector of Gaussian values: amp * exp(-(f-center)²/(2*std²))

get_slices(x::AbstractArray{T}; dims::Int=ndims(x)) where {T}

Get array slices along all dimensions except the specified one.

Arguments

• x::AbstractArray{T}: Input array
• dims::Int=ndims(x): Dimension to exclude when taking slices

Returns

• Iterator of array slices

Notes

• Returns slices along all dimensions except dims
• Each slice represents a view into the array with dims fixed

lorentzian(f, u)

Compute Lorentzian function values.

Arguments

• f::AbstractVector: Frequency values
• u::Vector: Parameters [amplitude, knee_frequency]

Returns

• Vector of Lorentzian values: amp/(1 + (f/knee)²)

lorentzian_initial_guess(psd, freqs; min_freq=freqs[1], max_freq=freqs[end])

Estimate initial parameters for Lorentzian fitting.

Arguments

• psd::AbstractVector{<:Real}: Power spectral density values
• freqs::AbstractVector{<:Real}: Frequency values
• min_freq::Real: Minimum frequency to consider
• max_freq::Real: Maximum frequency to consider

Returns

• Vector{Float64}: Initial guess for [amplitude, knee_frequency]

Notes

• Estimates amplitude from average power of low frequencies.
• Estimates knee frequency from half-power point.
• If allowvariableexponent=true, sets initial guess for exponent to 2.0.
• Used as starting point for nonlinear fitting

lorentzian_with_exponent(f, u)

Compute Lorentzian function values that allow variable exponent (PLE).

Arguments

• f::AbstractVector: Frequency values
• u::Vector: Parameters [amplitude, knee_frequency, exponent]

Returns

• Vector of Lorentzian values: amp/(1 + (f/knee)^exponent)

Notes

• Generalizes standard Lorentzian to allow variable power law exponent
• When exponent=2, reduces to standard Lorentzian

Residual function for expdecay du: residual u: parameters p: data

stack_and_reshape(slices; dims::Int)

Reconstruct an array from slices by stacking and reshaping them back to original dimensions.

Arguments

• slices: Iterator or collection of array slices (from get_slices)
• dims::Int: Dimension along which the original slicing was performed (should be same as dims fed into get_slices)

Returns

• Reconstructed array with proper dimension ordering

Example

# Split array into slices along dimension 2
x = rand(3, 4, 5)
slices = get_slices(x; dims=2)

# Reconstruct original array
x_reconstructed = stack_and_reshape(slices; dims=2)

OrnsteinUhlenbeck

Module for generating Ornstein-Uhlenbeck processes with various configurations. Uses DifferentialEquations.jl.

generate_ou_process(tau, true_D, dt, duration, num_trials; standardize=true, rng=Random.default_rng(), deq_seed=nothing)

Generate an Ornstein-Uhlenbeck process with a single timescale

Arguments

• tau::Union{Real, Vector{<:Real}}: Timescale(s) of the OU process
• true_D::Real: Target variance for scaling the process
• dt::Real: Time step size
• duration::Real: Total time length
• num_trials::Real: Number of trials/trajectories
• standardize::Bool=true: Whether to standardize output to match true_D
• rng::AbstractRNG=Random.default_rng(): Random number generator for initial conditions
• deq_seed::Integer=nothing: Random seed for DifferentialEquations.jl solver. If nothing, uses StochasticDiffEq.jl defaults. Note that for full replicability,

you need to set both rng and deq_seed.

Returns

• Matrix{Float64}: Generated OU process data with dimensions (numtrials, numtimesteps)

Notes

• Uses generateouprocess_sciml internally
• Returns NaN matrix if SciML solver fails
• Standardizes output to have specified variance if standardize=true

generate_ou_process_sciml(tau, true_D, dt, duration, num_trials, standardize; rng=Random.default_rng(), deq_seed=nothing)

Generate an Ornstein-Uhlenbeck process using DifferentialEquations.jl.

Arguments

• tau::Union{T, Vector{T}}: Timescale(s) of the OU process
• true_D::Real: Target variance for scaling
• dt::Real: Time step size
• duration::Real: Total time length
• num_trials::Integer: Number of trials/trajectories
• standardize::Bool: Whether to standardize output to match true_D
• rng::AbstractRNG=Random.default_rng(): Random number generator for initial conditions
• deq_seed::Union{Integer, Nothing}=nothing: Random seed for DifferentialEquations.jl solver. If nothing, uses StochasticDiffEq.jl defaults. Note that for full replicability,

you need to set both rng and deq_seed.

Returns

• Tuple{Matrix{Float64}, ODESolution}:
  • Scaled OU process data
  • Full SDE solution object

Notes

• Switches between static and dynamic arrays based on num_trials

Example:

tau = 1.0
true_D = 1.0
dt = 0.01
duration = 10.0
num_trials = 100

ou, _ = generate_ou_process_sciml(tau, true_D, dt, duration, num_trials, true)

# Reproducible example
deq_seed = 42
ou, _ = generate_ou_process_sciml(tau, true_D, dt, duration, num_trials, true, rng=Xoshiro(42), deq_seed=deq_seed)

generate_ou_with_oscillation(theta, dt, duration, num_trials, data_mean, data_sd; rng=Random.default_rng(), deq_seed=nothing)

Generate a one-timescale OU process with an additive oscillation.

Arguments

• theta::Vector{T}: Parameters [timescale, frequency, coefficient]
• dt::Real: Time step size
• duration::Real: Total time length
• num_trials::Integer: Number of trials
• data_mean::Real: Target mean value
• data_sd::Real: Target standard deviation
• rng::AbstractRNG=Random.default_rng(): Random number generator for initial conditions
• deq_seed::Union{Integer, Nothing}=nothing: Random seed for DifferentialEquations.jl solver. If nothing, uses StochasticDiffEq.jl defaults. Note that for full replicability,

you need to set both rng and deq_seed.

Returns

• Matrix{Float64}: Generated data with dimensions (numtrials, numtimesteps)

Notes

• Coefficient is bounded between 0 and 1
• Combines OU process with sinusoidal oscillation
• Standardizes and scales output to match target mean and standard deviation
• Returns NaN matrix if SciML solver fails

OneTimescale

Module for inferring a single timescale from time series data using the Ornstein-Uhlenbeck process.

OneTimescaleModel <: AbstractTimescaleModel

Model for inferring a single timescale from time series data using the Ornstein-Uhlenbeck process. We recommend using the one_timescale_model constructor function rather than creating directly.

Fields

• data::AbstractArray{<:Real}: Input time series data
• time::AbstractVector{<:Real}: Time points corresponding to the data
• fit_method::Symbol: Fitting method (:abc or :advi)
• summary_method::Symbol: Summary statistic type (:psd or :acf)
• lags_freqs: Lags (for ACF) or frequencies (for PSD)
• prior: Prior distribution(s) for parameters
• distance_method::Symbol: Distance metric type (:linear or :logarithmic)
• data_sum_stats: Pre-computed summary statistics
• dt::Real: Time step between observations
• T::Real: Total time span
• numTrials::Real: Number of trials/iterations
• data_mean::Real: Mean of input data
• data_sd::Real: Standard deviation of input data
• freqlims: Frequency limits for PSD analysis
• n_lags: Number of lags for ACF
• freq_idx: Boolean mask for frequency selection
• dims::Int: Dimension along which to compute statistics
• distance_combined::Bool: Whether to use combined distance metric
• weights::Vector{Real}: Weights for combined distance
• data_tau::Union{Real, Nothing}: Pre-computed timescale
• u0::Union{Vector{Real}, Nothing}: Initial parameter guess

Models.distance_function(model::OneTimescaleModel, sum_stats, data_sum_stats)

Calculate the distance between summary statistics of simulated and observed data.

Arguments

• model::OneTimescaleModel: Model instance
• sum_stats: Summary statistics from simulated data
• data_sum_stats: Summary statistics from observed data

Returns

• Distance value based on model.distancemethod (:linear or :logarithmic) or combined distance if model.distancecombined is true

Notes

If distance_combined is true:

• For ACF: Combines ACF distance with fitted exponential decay timescale distance
• For PSD: Combines PSD distance with knee frequency timescale distance

Models.generate_data(model::OneTimescaleModel, theta)

Generate synthetic data from the Ornstein-Uhlenbeck process with given timescale.

Arguments

• model::OneTimescaleModel: Model instance containing simulation parameters
• theta: Vector containing single timescale parameter (τ)

Returns

• Synthetic time series data with same dimensions as model.data

Notes

Uses the model's stored parameters:

• data_sd: Standard deviation for the OU process
• dt: Time step
• T: Total time span
• numTrials: Number of trials/trajectories

int_fit(model::OneTimescaleModel, param_dict::Dict=Dict())

Perform inference using the specified fitting method.

Arguments

• model::OneTimescaleModel: Model instance
• param_dict::Dict=Dict(): Dictionary of algorithm parameters. If empty, uses defaults.

Returns

For ABC method:

• posterior_samples: Matrix of accepted parameter samples
• posterior_MAP: Maximum a posteriori estimate
• abc_record: Full record of ABC iterations

For ADVI method:

• ADVIResults: Container with samples, MAP estimates, variances, and full chain

Notes

• For ABC: Uses Population Monte Carlo ABC with adaptive epsilon selection
• For ADVI: Uses Automatic Differentiation Variational Inference via Turing.jl
• Parameter dictionary can be customized for each method (see getparamdict_abc())

Models.summary_stats(model::OneTimescaleModel, data)

Compute summary statistics (ACF or PSD) from time series data.

Arguments

• model::OneTimescaleModel: Model instance specifying summary statistic type
• data: Time series data to analyze

Returns

For ACF (summary_method = :acf):

• Mean autocorrelation function up to n_lags

For PSD (summary_method = :psd):

• Mean power spectral density within specified frequency range

Notes

• ACF is computed using FFT-based method
• PSD is computed and filtered according to model.freq_idx
• Throws ArgumentError if summary_method is invalid

combined_distance(model::OneTimescaleModel, simulation_summary, data_summary,
                 weights, data_tau, simulation_tau)

Compute combined distance metric between simulated and observed data.

Arguments

• model: OneTimescaleModel instance
• simulation_summary: Summary statistics from simulation
• data_summary: Summary statistics from observed data
• weights: Weights for combining distances
• data_tau: Timescale from observed data
• simulation_tau: Timescale from simulation

Returns

• Weighted combination of summary statistic distance and timescale distance

one_timescale_model(data, time, fit_method; kwargs...)

Construct a OneTimescaleModel for time series analysis.

Arguments

• data: Input time series data
• time: Time points corresponding to the data
• fit_method: Fitting method to use (:abc or :advi)

Keyword Arguments

• summary_method=:acf: Summary statistic type (:psd or :acf)
• data_sum_stats=nothing: Pre-computed summary statistics
• lags_freqs=nothing: Custom lags or frequencies
• prior=nothing: Prior distribution(s) for parameters
• n_lags=nothing: Number of lags for ACF
• distance_method=nothing: Distance metric type
• dt=time[2]-time[1]: Time step
• T=time[end]: Total time span
• numTrials=size(data,1): Number of trials
• data_mean=mean(data): Data mean
• data_sd=std(data): Data standard deviation
• freqlims=nothing: Frequency limits for PSD
• freq_idx=nothing: Frequency selection mask
• dims=ndims(data): Analysis dimension
• distance_combined=false: Use combined distance
• weights=[0.5, 0.5]: Distance weights
• data_tau=nothing: Pre-computed timescale
• u0=nothing: Initial parameter guess

Returns

• OneTimescaleModel: Model instance configured for specified analysis method

Notes

Two main usage patterns:

1. ACF-based inference: summary_method=:acf, fit_method=:abc/:advi
2. PSD-based inference: summary_method=:psd, fit_method=:abc/:advi

OneTimescaleAndOscModel <: AbstractTimescaleModel

Model for inferring a single timescale and oscillation from time series data using the Ornstein-Uhlenbeck process. Parameters: [tau, freq, coeff] representing timescale, oscillation frequency, and oscillation coefficient.

Fields

• data::AbstractArray{<:Real}: Input time series data
• time::AbstractVector{<:Real}: Time points corresponding to the data
• fit_method::Symbol: Fitting method (:abc or :advi)
• summary_method::Symbol: Summary statistic type (:psd or :acf)
• lags_freqs: Lags (for ACF) or frequencies (for PSD)
• prior: Prior distribution(s) for parameters
• distance_method::Symbol: Distance metric type (:linear or :logarithmic)
• data_sum_stats: Pre-computed summary statistics
• dt::Real: Time step between observations
• T::Real: Total time span
• numTrials::Real: Number of trials/iterations
• data_mean::Real: Mean of input data
• data_sd::Real: Standard deviation of input data
• freqlims: Frequency limits for PSD analysis
• n_lags: Number of lags for ACF
• freq_idx: Boolean mask for frequency selection
• dims::Int: Dimension along which to compute statistics
• distance_combined::Bool: Whether to use combined distance metric
• weights::Vector{Real}: Weights for combined distance
• data_tau::Union{Real, Nothing}: Pre-computed timescale
• data_osc::Union{Real, Nothing}: Pre-computed oscillation frequency

Models.distance_function(model::OneTimescaleAndOscModel, sum_stats, data_sum_stats)

Calculate the distance between summary statistics of simulated and observed data.

Arguments

• model::OneTimescaleAndOscModel: Model instance
• sum_stats: Summary statistics from simulated data
• data_sum_stats: Summary statistics from observed data

Returns

• Distance value based on model.distancemethod (:linear or :logarithmic) or combined distance if model.distancecombined is true

Notes

If distance_combined is true:

• For ACF: Combines ACF distance with fitted exponential decay timescale distance
• For PSD: Combines PSD distance with knee frequency timescale distance and oscillation frequency distance

Models.generate_data(model::OneTimescaleAndOscModel, theta::AbstractVector{<:Real})

Generate synthetic data from the Ornstein-Uhlenbeck process with oscillation.

Arguments

• model::OneTimescaleAndOscModel: Model instance containing simulation parameters
• theta::AbstractVector{<:Real}: Vector containing parameters [tau, freq, coeff]

Returns

• Synthetic time series data with same dimensions as model.data

Notes

Uses the model's stored parameters:

• dt: Time step
• T: Total time span
• numTrials: Number of trials/trajectories
• data_mean: Mean of the process
• data_sd: Standard deviation of the process

int_fit(model::OneTimescaleAndOscModel, param_dict::Dict=Dict())

Perform inference using the specified fitting method.

Arguments

• model::OneTimescaleAndOscModel: Model instance
• param_dict::Dict=Dict(): Dictionary of algorithm parameters. If empty, uses defaults.

Returns

For ABC method:

• abc_record: Complete record of ABC iterations including posterior samples and MAP estimate

For ADVI method:

• result: Results from variational inference containing samples, MAP estimates, and variational posterior

Notes

• For ABC: Uses Population Monte Carlo ABC with adaptive epsilon selection
• For ADVI: Uses Automatic Differentiation Variational Inference via Turing.jl
• Parameter dictionary can be customized for each method (see getparamdictabc() and getparamdictadvi())

Models.summary_stats(model::OneTimescaleAndOscModel, data)

Compute summary statistics (ACF or PSD) from time series data.

Arguments

• model::OneTimescaleAndOscModel: Model instance specifying summary statistic type
• data: Time series data to analyze

Returns

For ACF (summary_method = :acf):

• Mean autocorrelation function up to n_lags

For PSD (summary_method = :psd):

• Mean power spectral density within specified frequency range

Notes

• ACF is computed using FFT-based method
• PSD is computed using AD-friendly implementation
• Throws ArgumentError if summary_method is invalid

combined_distance(model::OneTimescaleAndOscModel, simulation_summary, data_summary,
                 weights, distance_method, data_tau, simulation_tau, 
                 data_osc, simulation_osc)

Compute combined distance metric between simulated and observed data.

Arguments

• model: OneTimescaleAndOscModel instance
• simulation_summary: Summary statistics from simulation
• data_summary: Summary statistics from observed data
• weights: Weight vector for combining distances. For ACF: [w1, w2]. For PSD: [w1, w2, w3]
• distance_method: Distance metric type (:linear or :logarithmic)
• data_tau: Timescale from observed data
• simulation_tau: Timescale from simulation
• data_osc: Oscillation frequency from observed data (used for PSD only)
• simulation_osc: Oscillation frequency from simulation (used for PSD only)

Returns

For ACF:

• Weighted combination: weights[1] * summarydistance + weights[2] * timescaledistance

For PSD:

• Weighted combination: weights[1] * summarydistance + weights[2] * timescaledistance + weights[3] * oscillation_distance

one_timescale_and_osc_model(data, time, fit_method; kwargs...)

Construct a OneTimescaleAndOscModel for time series with oscillations. See https://duodenum96.github.io/IntrinsicTimescales.jl/stable/onetimescaleand_osc/ for details and complete examples.

Arguments

• data: Input time series data
• time: Time points corresponding to the data
• fit_method: Fitting method to use (:abc or :advi)

Keyword Arguments

• summary_method=:psd: Summary statistic type (:psd or :acf)
• data_sum_stats=nothing: Pre-computed summary statistics
• lags_freqs=nothing: Custom lags or frequencies
• prior=nothing: Prior distribution(s) for parameters
• n_lags=nothing: Number of lags for ACF
• distance_method=nothing: Distance metric type
• dt=time[2]-time[1]: Time step
• T=time[end]: Total time span
• numTrials=size(data,1): Number of trials
• data_mean=mean(data): Data mean
• data_sd=std(data): Data standard deviation
• freqlims=nothing: Frequency limits for PSD
• freq_idx=nothing: Frequency selection mask
• dims=ndims(data): Analysis dimension
• distance_combined=false: Use combined distance
• weights=[0.5, 0.5]: Distance weights for combined distance
• data_tau=nothing: Pre-computed timescale
• data_osc=nothing: Pre-computed oscillation frequency

Returns

• OneTimescaleAndOscModel: Model instance configured for specified analysis method

Notes

Four main usage patterns:

1. ACF-based ABC/ADVI: summary_method=:acf, fit_method=:abc/:advi
2. PSD-based ABC/ADVI: summary_method=:psd, fit_method=:abc/:advi

OneTimescaleWithMissing

Module for handling time series analysis with missing data. Uses specialized methods for handling NaN values:

• For ACF: Uses compactime_missing (equivalent to statsmodels.tsa.statstools.acf with missing="conservative")
• For PSD: Uses Lomb-Scargle periodogram to handle irregular sampling

OneTimescaleWithMissingModel <: AbstractTimescaleModel

Model for inferring a single timescale from time series data with missing values.

Fields

• data::AbstractArray{<:Real}: Input time series data (may contain NaN)
• time::AbstractVector{<:Real}: Time points corresponding to the data
• fit_method::Symbol: Fitting method (:abc or :advi)
• summary_method::Symbol: Summary statistic type (:psd or :acf)
• lags_freqs: Lags (for ACF) or frequencies (for PSD)
• prior: Prior distribution(s) for parameters
• distance_method::Symbol: Distance metric type (:linear or :logarithmic)
• data_sum_stats: Pre-computed summary statistics
• dt::Real: Time step between observations
• T::Real: Total time span
• numTrials::Real: Number of trials/iterations
• data_mean::Real: Mean of input data (excluding NaN)
• data_sd::Real: Standard deviation of input data (excluding NaN)
• freqlims: Frequency limits for PSD analysis
• n_lags: Number of lags for ACF
• freq_idx: Boolean mask for frequency selection
• dims::Int: Dimension along which to compute statistics
• distance_combined::Bool: Whether to use combined distance metric
• weights::Vector{Real}: Weights for combined distance
• data_tau::Union{Real, Nothing}: Pre-computed timescale
• missing_mask::AbstractArray{Bool}: Boolean mask indicating NaN positions

Models.distance_function(model::OneTimescaleWithMissingModel, sum_stats, data_sum_stats)

Calculate the distance between summary statistics of simulated and observed data.

Arguments

• model::OneTimescaleWithMissingModel: Model instance
• sum_stats: Summary statistics from simulated data
• data_sum_stats: Summary statistics from observed data

Returns

• Distance value based on model.distancemethod (:linear or :logarithmic) or combined distance if model.distancecombined is true

Notes

If distance_combined is true:

• For ACF: Combines ACF distance with fitted exponential decay timescale distance
• For PSD: Combines PSD distance with knee frequency timescale distance

Models.generate_data(model::OneTimescaleWithMissingModel, theta)

Generate synthetic data from the Ornstein-Uhlenbeck process and apply missing data mask.

Arguments

• model::OneTimescaleWithMissingModel: Model instance containing simulation parameters
• theta: Vector containing single timescale parameter (τ)

Returns

• Synthetic time series data with NaN values at positions specified by model.missing_mask

Notes

1. Generates complete OU process data
2. Applies missing data mask from original data
3. Returns data with same missing value pattern as input

int_fit(model::OneTimescaleWithMissingModel, param_dict=Dict())

Perform inference using the specified fitting method.

Arguments

• model::OneTimescaleWithMissingModel: Model instance
• param_dict=nothing: Optional dictionary of algorithm parameters. If nothing, uses defaults.

Returns

For ABC method:

• posterior_samples: Matrix of accepted parameter samples
• posterior_MAP: Maximum a posteriori estimate
• abc_record: Full record of ABC iterations

For ADVI method:

• ADVIResults: Container with samples, MAP estimates, variances, and full chain

Notes

• For ABC: Uses Population Monte Carlo ABC with adaptive epsilon selection
• For ADVI: Uses Automatic Differentiation Variational Inference via Turing.jl
• Parameter dictionary can be customized for each method (see getparamdict_abc())

Models.summary_stats(model::OneTimescaleWithMissingModel, data)

Compute summary statistics (ACF or PSD) from time series data with missing values.

Arguments

• model::OneTimescaleWithMissingModel: Model instance specifying summary statistic type
• data: Time series data to analyze (may contain NaN)

Returns

For ACF (summary_method = :acf):

• Mean autocorrelation function up to n_lags, computed with missing data handling

For PSD (summary_method = :psd):

• Mean Lomb-Scargle periodogram within specified frequency range

Notes

• ACF uses compactime_missing for proper handling of NaN values
• PSD uses Lomb-Scargle periodogram for irregular sampling

combined_distance(model, simulation_summary, data_summary, weights, data_tau, simulation_tau)

Calculate a weighted combination of summary statistic distance and timescale distance.

Arguments

• model::OneTimescaleWithMissingModel: Model instance
• simulation_summary: Summary statistics from simulated data
• data_summary: Summary statistics from observed data
• weights: Weight vector for combining distances
• data_tau: Timescale extracted from observed data
• simulation_tau: Timescale extracted from simulated data

Returns

• Combined weighted distance value

Notes

Uses the model's distance_method (:linear or :logarithmic) for summary statistic comparison and linear distance for timescale comparison.

one_timescale_with_missing_model(data, time, fit_method; kwargs...)

Construct a OneTimescaleWithMissingModel for time series analysis with missing data. See https://duodenum96.github.io/IntrinsicTimescales.jl/stable/onetimescalewith_missing/ for details and complete examples.

Arguments

• data: Input time series data (may contain NaN)
• time: Time points corresponding to the data
• fit_method: Fitting method to use (:abc or :advi)

Keyword Arguments

• summary_method=:acf: Summary statistic type (:psd or :acf)
• data_sum_stats=nothing: Pre-computed summary statistics
• lags_freqs=nothing: Custom lags or frequencies
• prior=nothing: Prior distribution(s) for parameters
• n_lags=nothing: Number of lags for ACF
• distance_method=nothing: Distance metric type
• dt=time[2]-time[1]: Time step
• T=time[end]: Total time span
• numTrials=size(data,1): Number of trials
• data_mean=nanmean(data): Data mean (excluding NaN)
• data_sd=nanstd(data): Data standard deviation (excluding NaN)
• freqlims=nothing: Frequency limits for PSD
• freq_idx=nothing: Frequency selection mask
• dims=ndims(data): Analysis dimension
• distance_combined=false: Use combined distance
• weights=[0.5, 0.5]: Distance weights
• data_tau=nothing: Pre-computed timescale

Returns

• OneTimescaleWithMissingModel: Model instance configured for specified analysis method

Notes

Four main usage patterns:

1. ACF-based ABC/ADVI: summary_method=:acf, fit_method=:abc/:advi
2. PSD-based ABC/ADVI: summary_method=:psd, fit_method=:abc/:advi

OneTimescaleAndOscWithMissing

Module for handling time series analysis with both oscillations and missing data. Uses specialized methods for handling NaN values:

• For ACF: Uses compactime_missing for proper handling of gaps
• For PSD: Uses Lomb-Scargle periodogram for irregular sampling

OneTimescaleAndOscWithMissingModel <: AbstractTimescaleModel

Model for inferring a single timescale from time series data with missing values. Parameters: [tau, freq, coeff] representing timescale, oscillation frequency, and oscillation coefficient.

Fields

• data::AbstractArray{<:Real}: Input time series data (may contain NaN)
• time::AbstractVector{<:Real}: Time points corresponding to the data
• fit_method::Symbol: Fitting method (:abc, :advi)
• summary_method::Symbol: Summary statistic type (:psd or :acf)
• lags_freqs: Lags (for ACF) or frequencies (for PSD)
• prior: Prior distribution(s) for parameters
• distance_method::Symbol: Distance metric type (:linear or :logarithmic)
• data_sum_stats: Pre-computed summary statistics
• dt::Real: Time step between observations
• T::Real: Total time span
• numTrials::Real: Number of trials/iterations
• data_mean::Real: Mean of input data (excluding NaN)
• data_sd::Real: Standard deviation of input data (excluding NaN)
• freqlims: Frequency limits for PSD analysis
• n_lags: Number of lags for ACF
• freq_idx: Boolean mask for frequency selection
• dims::Int: Dimension along which to compute statistics
• distance_combined::Bool: Whether to use combined distance metric
• weights::Vector{Real}: Weights for combined distance
• data_tau::Union{Real, Nothing}: Pre-computed timescale
• data_osc::Union{Real, Nothing}: Pre-computed oscillation frequency
• missing_mask::AbstractArray{Bool}: Boolean mask indicating NaN positions

Models.distance_function(model::OneTimescaleAndOscWithMissingModel, sum_stats, data_sum_stats)

Compute distance between simulated and observed summary statistics.

Arguments

• model::OneTimescaleAndOscWithMissingModel: Model instance specifying distance configuration
• sum_stats: Summary statistics from simulation
• data_sum_stats: Summary statistics from observed data

Returns

• Distance value based on model configuration

Notes

• If distancecombined=true: Uses combineddistance with both the estimated parameters and the full summary statistics.
• If distance_combined=false: Uses simple distance based on full summary statistics.
• Distance methods: :linear (Euclidean) or :logarithmic (log-space Euclidean)

Models.generate_data(model::OneTimescaleAndOscWithMissingModel, theta)

Generate synthetic data from the Ornstein-Uhlenbeck process with oscillation and apply missing data mask.

Arguments

• model::OneTimescaleAndOscWithMissingModel: Model instance containing simulation parameters
• theta: Vector containing parameters [tau, freq, coeff]

Returns

• Synthetic time series data with oscillations and NaN values at positions specified by model.missing_mask

Notes

1. Generates complete OU process data with oscillation
2. Applies missing data mask from original data
3. Returns data with same missing value pattern as input

int_fit(model::OneTimescaleAndOscWithMissingModel, param_dict=Dict())

Perform inference using the specified fitting method.

Arguments

• model::OneTimescaleAndOscWithMissingModel: Model instance
• param_dict=Dict(): Dictionary of algorithm parameters. If empty, uses defaults.

Returns

For ABC method:

• abc_record: Complete ABC record containing accepted samples, distances, and convergence info

For ADVI method:

• ADVIResults: Container with samples, MAP estimates, variances, and convergence information

Notes

• For ABC: Uses Population Monte Carlo ABC with adaptive epsilon selection
• For ADVI: Uses Automatic Differentiation Variational Inference via Turing.jl
• Default parameters available via getparamdictabc() and getparamdictadvi()

Models.summary_stats(model::OneTimescaleAndOscWithMissingModel, data)

Compute summary statistics (ACF or PSD) from time series data with missing values.

Arguments

• model::OneTimescaleAndOscWithMissingModel: Model instance specifying summary statistic type
• data: Time series data to analyze (may contain NaN)

Returns

For ACF (summary_method = :acf):

• Mean autocorrelation function up to n_lags, computed with missing data handling

For PSD (summary_method = :psd):

• Mean Lomb-Scargle periodogram within specified frequency range

Notes

• ACF uses compactime_missing for proper handling of NaN values
• PSD uses Lomb-Scargle periodogram for irregular sampling

combined_distance(model::OneTimescaleAndOscWithMissingModel, simulation_summary, data_summary,
                 weights, distance_method, data_tau, simulation_tau, data_osc, simulation_osc)

Compute combined distance metric between simulated and observed data.

Arguments

• model: OneTimescaleAndOscWithMissingModel instance
• simulation_summary: Summary statistics from simulation
• data_summary: Summary statistics from observed data
• weights: Weights for combining distances (length 2 for ACF, length 3 for PSD)
• distance_method: Distance metric type (:linear or :logarithmic)
• data_tau: Timescale from observed data
• simulation_tau: Timescale from simulation
• data_osc: Oscillation frequency from observed data (used only for PSD)
• simulation_osc: Oscillation frequency from simulation (used only for PSD)

Returns

For ACF:

• Weighted combination: weights[1] * summarydistance + weights[2] * taudistance

For PSD:

• Weighted combination: weights[1] * summarydistance + weights[2] * taudistance + weights[3] * osc_distance

Notes

• Ensure weights vector has correct length: 2 for ACF, 3 for PSD method

one_timescale_and_osc_with_missing_model(data, time, fit_method; kwargs...)

Construct a OneTimescaleAndOscWithMissingModel for inferring timescales and oscillations from data with missing values.

Arguments

• data: Input time series data (may contain NaN values)
• time: Time points corresponding to the data
• fit_method: Fitting method to use (:abc or :advi)

Keyword Arguments

• summary_method=:psd: Summary statistic type (:psd or :acf)
• data_sum_stats=nothing: Pre-computed summary statistics
• lags_freqs=nothing: Custom lags or frequencies
• prior=nothing: Prior distribution(s) for parameters
• n_lags=nothing: Number of lags for ACF
• distance_method=nothing: Distance metric type (:linear or :logarithmic)
• dt=time[2]-time[1]: Time step between observations
• T=time[end]: Total time span
• data_mean=nanmean(data): Data mean (excluding NaN values)
• data_sd=nanstd(data): Data standard deviation (excluding NaN values)
• freqlims=nothing: Frequency limits for PSD analysis (defaults to (0.5/1000, 100/1000) kHz)
• freq_idx=nothing: Frequency selection mask
• dims=ndims(data): Dimension along which to compute statistics
• numTrials=size(data, setdiff([1,2], dims)[1]): Number of trials/iterations
• distance_combined=false: Whether to use combined distance metric
• weights=[0.5, 0.5]: Weights for combined distance (length 2 for ACF, length 3 for PSD)
• data_tau=nothing: Pre-computed timescale for combined distance
• data_osc=nothing: Pre-computed oscillation frequency for combined distance

Returns

• OneTimescaleAndOscWithMissingModel: Model instance configured for specified analysis method

Notes

Main usage patterns:

1. ACF-based analysis: summary_method=:acf - Uses autocorrelation with missing data handling
2. PSD-based analysis: summary_method=:psd - Uses Lomb-Scargle periodogram for irregular sampling

Key differences from regular model:

• Automatically detects and handles NaN values in input data
• Uses nanmean/nanstd for robust statistics computation
• Applies specialized missing-data algorithms for summary statistics
• Preserves missing data pattern in synthetic data generation

For combined distance (distance_combined=true):

• ACF method: Combines summary distance with timescale distance (2 weights)
• PSD method: Combines summary, timescale, and oscillation distances (3 weights)

acwplot(acwresults::ACWResults; only_acf::Bool=false, only_psd::Bool=false, show::Bool=true)

Placeholder function for plotting autocorrelation function (ACF) and power spectral density (PSD) results.

This function is implemented in the IntPlottingExt extension package, which is automatically loaded when Plots.jl is available. The extension provides comprehensive visualization capabilities for ACW analysis results.

Functionality (available when Plots.jl is loaded):

• Plots autocorrelation function and/or power spectral density from ACWResults
• Supports plotting either ACF alone, PSD alone, or both in a subplot layout
• ACF plots show individual traces in color with mean overlaid in black
• PSD plots use logarithmic scales for both axes
• Handles up to 2-dimensional ACF data

Arguments (when extension is loaded):

• acwresults::ACWResults: Container with ACF and/or PSD analysis results
• only_acf::Bool=false: If true, plot only the autocorrelation function
• only_psd::Bool=false: If true, plot only the power spectral density
• show::Bool=true: Whether to display the plot immediately

Requirements:

• Requires Plots.jl to be loaded for the extension to activate
• Install with: using Plots or import Plots

See IntPlottingExt documentation for complete details.

posterior_predictive(container::Union{ABCResults, ADVIResults}, model::Models.AbstractTimescaleModel; show::Bool=true, n_samples::Int=100)

Placeholder function for posterior predictive check plotting.

This function is implemented in the IntPlottingExt extension package, which is automatically loaded when Plots.jl is available. The extension provides posterior predictive checking visualization for both ABC and ADVI inference results.

Functionality (available when Plots.jl is loaded):

• Creates posterior predictive check plots showing data vs. model predictions
• Works with both ABCResults and ADVIResults containers
• Shows posterior predictive intervals with uncertainty bands
• Automatically handles ACF and PSD summary statistics with appropriate scaling
• Supports logarithmic scaling when the distance method is logarithmic

Arguments:

• container::Union{ABCResults, ADVIResults}: Results container from inference
• model::Models.AbstractTimescaleModel: Model used for inference
• show::Bool=true: Whether to display the plot immediately
• n_samples::Int=100: Number of posterior samples to use for prediction

Requirements:

• Requires Plots.jl to be loaded for the extension to activate
• Install with: using Plots or import Plots

See IntPlottingExt documentation for complete details.