| Title: | Empirical Discovery of Differential Equations from Time Series Data |
|---|---|
| Description: | A comprehensive toolkit for discovering differential and difference equations from empirical time series data using symbolic regression. The package implements a complete workflow from data preprocessing (including Total Variation Regularized differentiation for noisy economic data), visual exploration of dynamical structure, and symbolic equation discovery via genetic algorithms. It leverages a high-performance 'Julia' backend ('SymbolicRegression.jl') to provide industrial-grade robustness, physics-informed constraints, and rigorous out-of-sample validation. Designed for economists, physicists, and researchers studying dynamical systems from observational data. |
| Authors: | José Mauricio Gómez Julián [aut, cre]
|
| Maintainer: | José Mauricio Gómez Julián <[email protected]> |
| License: | GPL (>= 3) |
| Version: | 0.1.14 |
| Built: | 2026-09-26 09:58:03 UTC |
| Source: | https://github.com/IsadoreNabi/EmpiricalDynamics |
Creates a concise text summary of the analysis.
analysis_summary(results, verbose = TRUE)analysis_summary(results, verbose = TRUE)
results |
Analysis results list. |
verbose |
Include additional details. |
Character string with summary.
Examines how fixed points change as a parameter varies.
analyze_bifurcations( equation, variable, parameter, param_range = c(-5, 5), n_param = 50, z_range = c(-10, 10), exogenous_values = list(), n_draws = 200L, interval = 0.9 )analyze_bifurcations( equation, variable, parameter, param_range = c(-5, 5), n_param = 50, z_range = c(-10, 10), exogenous_values = list(), n_draws = 200L, interval = 0.9 )
equation |
Fitted equation object: a |
variable |
Name of the main variable. |
parameter |
Name of the parameter to vary. Must be a fitted
coefficient or a variable of the equation, and must differ from
|
param_range |
Range for parameter values. |
n_param |
Number of parameter values to test. |
z_range |
Range for searching fixed points. |
exogenous_values |
Fixed values for other variables. |
n_draws |
Number of posterior draws to sweep at each parameter value
when |
interval |
Width of the posterior interval reported in
|
The parameter must actually appear in the equation, and is varied
according to what it is: a fitted coefficient is overridden through
coefficient_values of analyze_fixed_points, and a
free variable of the expression is bound through
exogenous_values. A name that is neither is an error – sweeping
it would evaluate the same equation n_param times and label the
identical copies with the parameter's name.
When one name could denote both: for lm/glm equations a
plain term's coefficient label coincides with the variable's name, and
the name is taken as the variable (labels that are never variables, like
"I(Z)", denote their coefficient); for nls and
symbolic_equation equations a coefficient name is a symbol of the
formula itself and is taken as the coefficient.
When equation carries posterior draws, the sweep is run once per
draw and the answer is a distribution of fixed points at every parameter
value rather than a point. This is the interventional reading of the
sweep – the coordinate is SET to each value in every draw, the remaining
uncertainty left as estimated – and not the conditional one, which would
keep only the draws that happen to agree with the value and would empty out
at exactly the extreme parameter values a bifurcation diagram exists to
visit. $data then carries a draw column and $summary
holds the branch-by-branch reduction.
Object of class "bifurcation_analysis". $data holds the
fixed points found; $status holds one row per requested
parameter value – "ok", "no_fixed_points" or
"error" with the error's message – so the requested grid can
always be reconstructed and a failure is never mistaken for an
absence; $mode records whether the parameter was varied as a
"coefficient" or as an "exogenous" variable.
$posterior says whether the sweep carries uncertainty. When it
does, $data gains a draw column and holds one block of rows
per draw per parameter value; $status gains n_draws and
prob_n_fixed_points, and its n_fixed_points is the count
most draws agreed on rather than the number of rows; and $summary
holds one row per parameter value per branch, with
fixed_point_lower and fixed_point_upper bounding a
$interval posterior interval. $summary is NULL for a
point estimate, where $data already is the answer.
# dZ = r*Z - Z^2: transcritical bifurcation at r = 0 set.seed(1) data <- data.frame(Z = seq(-1, 3, length.out = 80)) data$dZ <- 2 * data$Z - data$Z^2 + 0.01 * stats::rnorm(80) fit <- stats::nls(dZ ~ r * Z - Z^2, data = data, start = list(r = 1)) bif <- analyze_bifurcations(fit, variable = "Z", parameter = "r", param_range = c(-1, 1), n_param = 11, z_range = c(-3, 3)) print(bif)# dZ = r*Z - Z^2: transcritical bifurcation at r = 0 set.seed(1) data <- data.frame(Z = seq(-1, 3, length.out = 80)) data$dZ <- 2 * data$Z - data$Z^2 + 0.01 * stats::rnorm(80) fit <- stats::nls(dZ ~ r * Z - Z^2, data = data, start = list(r = 1)) bif <- analyze_bifurcations(fit, variable = "Z", parameter = "r", param_range = c(-1, 1), n_param = 11, z_range = c(-3, 3)) print(bif)
Finds and characterizes fixed points of the discovered equation.
analyze_fixed_points( equation, variable, range = c(-10, 10), n_grid = 100, exogenous_values = list(), coefficient_values = list(), n_draws = 200L )analyze_fixed_points( equation, variable, range = c(-10, 10), n_grid = 100, exogenous_values = list(), coefficient_values = list(), n_draws = 200L )
equation |
Fitted equation object ( |
variable |
Name of the main variable. |
range |
Numeric vector of length 2 specifying search range. |
n_grid |
Number of grid points for initial search. |
exogenous_values |
Named list of fixed values for exogenous variables. |
coefficient_values |
Named list of deliberate overrides for fitted coefficients; every name must match a fitted coefficient. |
n_draws |
Number of posterior draws to sweep when |
For lm and glm equations the right-hand side is evaluated
through the model's own predict() machinery (on the response scale
for glm), so the fitted coefficients enter exactly as fitted –
I() terms, factors and link functions included. For
symbolic_equation and nls equations the expression is
evaluated with the fitted coefficients bound by name.
A fitted object that carries posterior draws – an rstanarm fit, for
instance, whose class is stanreg, glm, lm – is swept through its
draws instead of through its point estimate, and the result then carries
one block of rows per draw with a draw column identifying them; see
the n_draws argument and the "posterior" attribute of the
return value. Objects are recognised by whether they can produce a draws
matrix covering their own coefficients, not by class name.
When coefficient_values is used, the substitution is verified rather
than assumed: the requested override must actually move the predictions
somewhere on the search grid. An object whose predict() does not read
$coefficients – which is true of every fit that stores its estimate
elsewhere and merely summarises it into $coefficients – is refused
with an "ed_substitution_ignored" condition instead of being swept
into a table of identical copies. A coefficient whose term is identically
zero over the search range is refused as "ed_substitution_inert".
Every name consulted by the equation must be accounted for: the search
variable, a fitted coefficient, or a value declared in
exogenous_values. A free symbol without a value, an unnamed
declaration, a declaration that collides with the search variable or with
a fitted coefficient – each of these is reported instead of being
silently ignored or silently overwriting the fit. To vary a fitted
coefficient deliberately, use coefficient_values.
Data frame of fixed points with stability classification. The
values any exogenous variables were held at travel with the result as
its "exogenous_values" attribute, because a fixed point found
under fixed exogenous values is conditional on that fixing; the
deliberate coefficient overrides travel as "coefficient_values".
Both attributes are attached whether or not any fixed point was found.
The "posterior" attribute says whether the sweep was run over
posterior draws; when it is TRUE the frame has a draw
column, holds one block of rows per draw, and "n_draws" gives how
many draws were swept. ed_consensus_fixed_points() reduces such a
frame to one row per branch.
# dZ = 2*Z - Z^2: fixed points at 0 (unstable, f' = +2) # and at 2 (stable, f' = -2) data <- data.frame(Z = seq(0.1, 3, length.out = 50)) data$dZ <- 2 * data$Z - data$Z^2 model <- stats::lm(dZ ~ I(Z) + I(Z^2) + 0, data = data) analyze_fixed_points(model, variable = "Z", range = c(-1, 3))# dZ = 2*Z - Z^2: fixed points at 0 (unstable, f' = +2) # and at 2 (stable, f' = -2) data <- data.frame(Z = seq(0.1, 3, length.out = 50)) data$dZ <- 2 * data$Z - data$Z^2 model <- stats::lm(dZ ~ I(Z) + I(Z^2) + 0, data = data) analyze_fixed_points(model, variable = "Z", range = c(-1, 3))
Records researcher hypotheses based on visual exploration, to be used as constraints or guides in subsequent symbolic search.
annotate_hypotheses(data, hypotheses)annotate_hypotheses(data, hypotheses)
data |
Data frame (with exploration results as attribute). |
hypotheses |
Character vector of hypotheses. |
Data frame with hypotheses attached as attribute.
# Toy example data <- data.frame(Z = 1:10) data <- annotate_hypotheses(data, c( "Z exhibits logistic saturation around Z=100", "Effect of X appears linear" ))# Toy example data <- data.frame(Z = 1:10) data <- annotate_hypotheses(data, c( "Z exhibits logistic saturation around Z=100", "Effect of X appears linear" ))
Computes bootstrap confidence intervals for equation parameters.
bootstrap_parameters( equation, data, response = NULL, derivative_col = NULL, n_boot = 500, conf_level = 0.95, block_size = NULL, weights = NULL )bootstrap_parameters( equation, data, response = NULL, derivative_col = NULL, n_boot = 500, conf_level = 0.95, block_size = NULL, weights = NULL )
equation |
Fitted equation object. |
data |
Original data. |
response |
Name of response column. |
derivative_col |
Alias for response. |
n_boot |
Number of bootstrap samples. |
conf_level |
Confidence level (default 0.95). |
block_size |
Block size for block bootstrap (time series). |
weights |
Optional vector of observation weights, one per row of
|
A discovered equation carries its constants as literals and has no
coefficients to resample. It is therefore parameterised first, exactly as
in cross_validate, and what is resampled are the parameters
that stand in for those literals.
Data frame with parameter estimates and confidence intervals.
Checks whether the discovered equation exhibits expected qualitative features, and refuses to be asked for anything it does not check.
check_qualitative_behavior( equation, data, variable, expected_features = list(), n_grid = 100, boundary_tol = ed_derivative_step(), n_draws = 200L )check_qualitative_behavior( equation, data, variable, expected_features = list(), n_grid = 100, boundary_tol = ed_derivative_step(), n_draws = 200L )
equation |
Fitted equation object. |
data |
Original data. Must contain |
variable |
Main variable name. |
expected_features |
Named list of expected qualitative features. The recognised names, and nothing else, are:
|
n_grid |
Number of grid points the field is read on, over the
observed range, when |
boundary_tol |
Magnitude below which the field is treated as zero when
its sign is read. The default is |
n_draws |
Number of posterior draws to sweep when |
Every name in expected_features must be one this function honours;
an unrecognised name is an error rather than a silent omission. It has to
be: a dropped name lowers the number of checks reported without saying
which check went missing, so 2/2 passed cannot be told apart from
2/3 passed with a typo in the third. Names are matched exactly,
never by the partial matching $ performs on lists.
The features are read over two different ranges, deliberately. Fixed points
are searched over the observed range of variable widened by 20\
a root just outside the data is not missed. Boundedness and monotonicity are
read over the observed range itself, unwidened, because they are statements
about the interval the data actually visited: widening it would ask about
states that were never observed, and on the textbook logistic
dZ/dt = r Z (1 - Z/K) the widened range dips below zero, where the
field points away, and would answer that the bounded system is unbounded.
What bounded claims. TRUE means the observed range is
invariant: the field points inward at both ends, f(min) > 0 and
f(max) < 0, so no trajectory starting inside the observed range ever
leaves it. In one dimension that is exact, and it is decided by algebra –
there is no horizon to choose, no trajectory to simulate and no seed. It is
stronger than global boundedness and does not follow from it: the
cubic -Z(Z - 0.5)(Z - 1.5) is globally bounded, yet from the low end
of a range observed above 0.5 the trajectory leaves that range on its way to
the fixed point at zero, and this check reports FALSE. An end that is
itself a fixed point within boundary_tol leaves the flow neither
entering nor leaving, and the check declines to decide rather than deciding
against.
What monotonicity claims. The sign of the field over the
observed range – "positive", "negative", or "none"
when it is neither strictly one nor strictly the other. It is therefore the
monotonicity in time of every trajectory that stays in the observed range,
and "none" is the answer whenever the range contains a fixed point.
Under a posterior. When equation carries posterior draws,
every feature is decided per draw and reported as the verdict most draws
agreed on, together with the fraction that held it – the same reduction
ed_consensus_fixed_points applies to the number of fixed
points. A tie breaks toward the weaker claim: not bounded before bounded,
undecidable before either.
Object of class "qualitative_check": $passed holds one
named logical per requested feature – always as many as were requested –
$messages the line-by-line verdicts, $checks$fixed_points
the roots found, $checks$field the ends of the observed range and
the values of the field there, and $posterior whether the checks
were run over posterior draws.
# Logistic truth: the observed range is invariant, and the field vanishes # inside it, so the trajectory is not monotone over the whole range. d <- data.frame(Z = seq(0.05, 1.9, length.out = 200)) d$dZ <- 1.2 * d$Z * (1 - d$Z / 1.5) fit <- stats::lm(dZ ~ I(Z) + I(Z^2) + 0, data = d) check_qualitative_behavior(fit, d, "Z", expected_features = list(bounded = TRUE, monotonicity = "none"))# Logistic truth: the observed range is invariant, and the field vanishes # inside it, so the trajectory is not monotone over the whole range. d <- data.frame(Z = seq(0.05, 1.9, length.out = 200)) d$dZ <- 1.2 * d$Z * (1 - d$Z / 1.5) fit <- stats::lm(dZ ~ I(Z) + I(Z^2) + 0, data = d) check_qualitative_behavior(fit, d, "Z", expected_features = list(bounded = TRUE, monotonicity = "none"))
Creates a publication-ready table of estimated coefficients.
coefficient_table( equation, bootstrap_results = NULL, format = c("data.frame", "latex", "markdown", "html"), caption = "Estimated Coefficients", label = "tab:coefficients" )coefficient_table( equation, bootstrap_results = NULL, format = c("data.frame", "latex", "markdown", "html"), caption = "Estimated Coefficients", label = "tab:coefficients" )
equation |
Fitted equation object. |
bootstrap_results |
Results from |
format |
Output format: "data.frame", "latex", "markdown", "html". |
caption |
Table caption. |
label |
LaTeX label for referencing. |
Formatted table.
Applies multiple differentiation methods to the same data and produces a comparison plot.
compare_differentiation_methods( Z, t = NULL, methods = c("tvr", "savgol", "spline", "finite_diff"), plot = TRUE )compare_differentiation_methods( Z, t = NULL, methods = c("tvr", "savgol", "spline", "finite_diff"), plot = TRUE )
Z |
Numeric vector of observations. |
t |
Numeric vector of time points. |
methods |
Character vector of methods to compare. |
plot |
Produce comparison plot? |
List of derivative vectors from each method.
Produces a comparison of drift estimates from ordinary least squares versus iterative GLS.
compare_estimation_methods(ols_model, gls_model, data)compare_estimation_methods(ols_model, gls_model, data)
ols_model |
SDE model estimated with OLS |
gls_model |
SDE model estimated with iterative GLS |
data |
Data frame used for estimation |
Invisibly returns comparison statistics
Computes metrics comparing simulated trajectories to observed data.
compare_trajectories( simulation, observed_data, time_col = "time", var_col = NULL )compare_trajectories( simulation, observed_data, time_col = "time", var_col = NULL )
simulation |
Trajectory simulation object. |
observed_data |
Data frame with observed values. |
time_col |
Name of time column. |
var_col |
Name of variable column to compare. |
Data frame with comparison metrics.
Main dispatcher function for numerical differentiation. Supports multiple methods appropriate for different data characteristics.
compute_derivative( Z, t = NULL, method = c("tvr", "savgol", "spline", "finite_diff", "spectral"), ... )compute_derivative( Z, t = NULL, method = c("tvr", "savgol", "spline", "finite_diff", "spectral"), ... )
Z |
Numeric vector of observations. |
t |
Numeric vector of time points (NULL assumes dt=1). |
method |
Differentiation method: "tvr", "savgol", "spline", "finite_diff", "fd", or "spectral". |
... |
Additional arguments passed to the specific method. |
Available methods:
tvr: Total Variation Regularized differentiation (recommended for economic data with trends and shocks).
savgol: Savitzky-Golay filter (moderate noise, preserves peaks).
spline: Smoothing spline (high noise, prioritizes trend).
finite_diff or fd: Centered finite differences (low noise).
spectral: FFT-based (periodic data only).
Numeric vector of estimated derivatives with diagnostic attributes.
compute_derivative_tvr, suggest_differentiation_method
t <- seq(0, 10, by = 0.1) Z <- sin(t) + rnorm(length(t), sd = 0.1) # Finite differences (fast, no dependencies) dZ_fd <- compute_derivative(Z, t, method = "finite_diff") # Access the derivative vector for plotting plot(t, dZ_fd$derivative, type = "l", main = "Derivative Comparison") lines(t, cos(t), col = "red", lty = 2) # True derivative # TVR (requires CVXR) if (requireNamespace("CVXR", quietly = TRUE)) { dZ_tvr <- compute_derivative(Z, t, method = "tvr") }t <- seq(0, 10, by = 0.1) Z <- sin(t) + rnorm(length(t), sd = 0.1) # Finite differences (fast, no dependencies) dZ_fd <- compute_derivative(Z, t, method = "finite_diff") # Access the derivative vector for plotting plot(t, dZ_fd$derivative, type = "l", main = "Derivative Comparison") lines(t, cos(t), col = "red", lty = 2) # True derivative # TVR (requires CVXR) if (requireNamespace("CVXR", quietly = TRUE)) { dZ_tvr <- compute_derivative(Z, t, method = "tvr") }
Computes derivatives using centered finite differences.
compute_derivative_fd(Z, t = NULL, ...)compute_derivative_fd(Z, t = NULL, ...)
Z |
Numeric vector of observations. |
t |
Numeric vector of time points. |
... |
Additional arguments (ignored). |
List with derivative vector.
Computes derivatives using the Savitzky-Golay filter.
compute_derivative_savgol(Z, t = NULL, p = 3, n = NULL, m = 1, ...)compute_derivative_savgol(Z, t = NULL, p = 3, n = NULL, m = 1, ...)
Z |
Numeric vector of observations. |
t |
Numeric vector of time points. |
p |
Polynomial order (default 3). |
n |
Filter length (must be odd, default auto-selected). |
m |
Derivative order (default 1). |
... |
Additional arguments (ignored). |
List with derivative vector.
Computes derivatives using the Fourier transform.
compute_derivative_spectral(Z, t = NULL, ...)compute_derivative_spectral(Z, t = NULL, ...)
Z |
Numeric vector of observations. |
t |
Numeric vector of time points. |
... |
Additional arguments (ignored). |
List with derivative vector.
This method assumes the signal is periodic.
Computes derivatives by fitting a smoothing spline and differentiating it.
compute_derivative_spline(Z, t = NULL, spar = NULL, df = NULL, ...)compute_derivative_spline(Z, t = NULL, spar = NULL, df = NULL, ...)
Z |
Numeric vector of observations. |
t |
Numeric vector of time points. |
spar |
Smoothing parameter (NULL for automatic selection). |
df |
Degrees of freedom (alternative to spar). |
... |
Additional arguments (ignored). |
List with derivative vector.
Computes derivatives by solving a convex optimization problem that balances fidelity to the data against smoothness of the derivative. The problem is internally rescaled for numerical stability, and a cascading solver chain is used for maximum robustness.
compute_derivative_tvr(Z, t = NULL, lambda = "auto", solver = "clarabel", ...)compute_derivative_tvr(Z, t = NULL, lambda = "auto", solver = "clarabel", ...)
Z |
Numeric vector of observations. |
t |
Numeric vector of time points (NULL assumes dt=1). |
lambda |
Regularization parameter ("auto" for cross-validation selection). |
solver |
Preferred solver: "clarabel" (default, interior point), "scs" (conic ADMM), or "osqp" (QP ADMM). If the preferred solver fails or returns an inaccurate solution, other solvers are tried automatically. |
... |
Additional arguments (ignored). |
Solver status interpretation:
"optimal": The solver converged within its tolerances.
Results are reliable.
"optimal_inaccurate": The solver found a solution but did not
meet strict tolerance criteria. For TVR this rarely affects derivative
quality in a perceptible way, because internal rescaling keeps the
problem well-conditioned. If degradation is observed, try a different
preferred solver or adjust lambda manually.
TVR and stochastic differential equations:
TVR penalizes total variation of the derivative, which effectively smooths
out high-frequency components. In SDEs, the diffusion signature lives in
those high-frequency residuals. As a consequence, the better TVR estimates
the drift, the more it removes the information needed to recover the
diffusion coefficient from residuals. For SDE applications, consider using
estimate_diffusion_qv (quadratic variation) to estimate
diffusion directly from increments rather than from TVR residuals.
Object of class "tvr_derivative" (also a list with $derivative).
Convenience function to compute derivatives for all endogenous variables in a specified dataset.
compute_derivatives(data, method = "tvr", prefix = "d_", ...)compute_derivatives(data, method = "tvr", prefix = "d_", ...)
data |
Data frame with variable specifications (from |
method |
Differentiation method. |
prefix |
Prefix for derivative column names (default "d_"). |
... |
Additional arguments passed to |
Data frame with derivative columns added.
Calculates residuals from a fitted symbolic equation or SDE model.
compute_residuals(model, data, target = NULL)compute_residuals(model, data, target = NULL)
model |
A symbolic_equation or sde_model object |
data |
Data frame containing the variables |
target |
Name of target variable (auto-detected if NULL) |
Numeric vector of residuals
Combines a drift equation and diffusion model into a complete SDE: dZ = f(Z, X, Y) dt + g(Z, X, Y) dW
construct_sde( drift, diffusion = NULL, variable = NULL, refine_with_gls = FALSE, gls_max_iter = 10, gls_tolerance = 1e-04, data = NULL, target = NULL )construct_sde( drift, diffusion = NULL, variable = NULL, refine_with_gls = FALSE, gls_max_iter = 10, gls_tolerance = 1e-04, data = NULL, target = NULL )
drift |
Symbolic equation for the drift term f(.) |
diffusion |
Variance model for the diffusion term g(.) |
variable |
Name of the main state variable |
refine_with_gls |
Use iterative GLS to refine estimates? |
gls_max_iter |
Maximum iterations for GLS |
gls_tolerance |
Convergence tolerance for GLS |
data |
Data frame (required if refine_with_gls = TRUE) |
target |
Target variable name (required if refine_with_gls = TRUE) |
An object of class "sde_model"
Generates derived variables based on specified transformations that may be theoretically relevant.
create_transformations( data, transformations = NULL, variables = NULL, cols = NULL )create_transformations( data, transformations = NULL, variables = NULL, cols = NULL )
data |
Data frame. |
transformations |
List of formulas specifying transformations. |
variables |
Character vector of variable names (alternative interface). |
cols |
Alias for variables. |
Data frame with transformation columns added.
data <- data.frame(X = 1:10, Y = 10:1) # Simple interface data <- create_transformations(data, variables = c("X", "Y")) # Formula interface data <- create_transformations(data, transformations = list( ratios = ~ X/Y ))data <- data.frame(X = 1:10, Y = 10:1) # Simple interface data <- create_transformations(data, variables = c("X", "Y")) # Formula interface data <- create_transformations(data, transformations = list( ratios = ~ X/Y ))
Performs k-fold or block cross-validation to assess out-of-sample predictive performance of the discovered equation.
cross_validate( equation, data, response = NULL, derivative_col = NULL, k = 5, method = c("block", "random", "rolling"), block_size = NULL, horizon = 1, window = c("expanding", "rolling", "sliding"), train_size = NULL, refit_derivative = FALSE, diff_method = "tvr", weights = NULL, refit_engine = c("r", "julia"), verbose = TRUE )cross_validate( equation, data, response = NULL, derivative_col = NULL, k = 5, method = c("block", "random", "rolling"), block_size = NULL, horizon = 1, window = c("expanding", "rolling", "sliding"), train_size = NULL, refit_derivative = FALSE, diff_method = "tvr", weights = NULL, refit_engine = c("r", "julia"), verbose = TRUE )
equation |
Fitted equation object from |
data |
Data frame containing all variables. |
response |
Name of the response column (derivative). |
derivative_col |
Alias for response (for compatibility). |
k |
Number of folds for cross-validation. |
method |
CV method: |
block_size |
For |
horizon |
For |
window |
For |
train_size |
Length of the training window when |
refit_derivative |
Logical; whether to recompute derivatives for each fold (currently unused). |
diff_method |
Differentiation method if refitting (currently unused). |
weights |
Optional vector of observation weights, one per row of
|
refit_engine |
Which engine re-estimates the constants on each training
fold: |
verbose |
Print progress. |
A discovered equation arrives with its constants already fitted, as
numeric literals, and there is nothing in it left to estimate. Before the
folds begin, such an equation is passed through
parameterize_equation: its literals become free parameters
starting at the values the search found, and each fold re-estimates them on
its own training rows. An equation specified by the user, which already has
named free parameters, keeps taking its starting values from its own
coefficients. An equation with no constants at all is not refitted – there
is nothing to refit – and is scored by evaluating it on the held-out rows;
the returned object says so in refitted.
Object of class "cv_result" containing:
rmse |
Root mean squared error per fold |
mae |
Mean absolute error per fold |
r_squared |
R-squared per fold |
mean_rmse |
Average RMSE across folds |
sd_rmse |
Standard deviation of RMSE |
predictions |
List of predicted vs actual per fold |
fold_indices |
Indices used for each fold |
refitted |
Whether the equation was re-estimated on each fold |
parameterization |
The parameterised expression that was refitted, when
the equation was a discovered one, and |
refit_engine |
The engine that did the refitting |
stalled_folds |
Folds in which the Julia optimiser did not move off its starting constants, so the error reported for them is that of the equation as the search left it. Empty for the R engine. |
train_indices |
The rows each fold trained on. Only for
|
n_folds_ok, prob_converged, failures
|
How many folds refitted, what
fraction that is, and the error message of each one that did not, named
by fold. |
rmse_given_converged, mae_given_converged, r_squared_given_converged
|
The same errors over the folds that refitted.
Read them together with |
r_squared_pooled |
A ratio of sums rather than a mean of ratios. |
rmse_with_fallback, r_squared_pooled_with_fallback
|
The error of the procedure rather than of the model: a fold whose refit failed is scored with the prediction a caller would have been left with, the training mean, which contributes exactly zero to the pooled R-squared. |
ss_res, ss_tot
|
Per fold, so that a nearly degenerate benchmark is visible instead of being rounded into or out of existence. |
R-squared is against the mean of each fold's TRAINING rows, which
is the constant a caller would have predicted without fitting anything and
the only one available before the held-out rows are seen. Using the test
fold's own mean, as this function did before, benchmarks the model against
a constant that can only be computed after the answers are in; makes every
fold a different quantity, so averaging them means nothing; and is exactly
zero for a single-row test fold, which is what "rolling" with the
default horizon produces and which returned -Inf. Under the present
definition a model that predicts the training mean scores exactly 0 and a
perfect model exactly 1.
The k-fold figures are NA unless all k folds refitted.
Averaging over the folds that converged is conditioning on success. The
conditional figures are there under their own names.
The folds are not independent of one another when the response is
a numerically estimated derivative, which is the usual case in this
package. compute_derivative builds each dZ[t] from
neighbouring rows – a window for "savgol" and "spline",
the whole series for "tvr" – so a held-out row's response is
partly a function of training rows' predictors. This does not bias a
comparison between candidate equations, which all receive the same favour,
but it does mean the figure is not the error to expect on a fresh series,
and least of all on a forecast, where the rows after t are not
available to differentiate with. Recomputing the derivative inside each
fold is what refit_derivative is reserved for, and it is not
implemented.
# Toy example using lm data <- data.frame( time = 1:50, y = seq(1, 10, length.out = 50) + stats::rnorm(50, sd = 0.1) ) # Simple linear model as a proxy for a discovered equation model <- stats::lm(y ~ time, data = data) # Run cross-validation cv_res <- cross_validate( equation = model, data = data, response = "y", k = 3, method = "random" ) print(cv_res)# Toy example using lm data <- data.frame( time = 1:50, y = seq(1, 10, length.out = 50) + stats::rnorm(50, sd = 0.1) ) # Simple linear model as a proxy for a discovered equation model <- stats::lm(y ~ time, data = data) # Run cross-validation cv_res <- cross_validate( equation = model, data = data, response = "y", k = 3, method = "random" ) print(cv_res)
Defines custom mathematical operators for use in symbolic search.
define_custom_operators(...)define_custom_operators(...)
... |
Named functions to add as operators. |
List of operator definitions suitable for symbolic_search.
ops <- define_custom_operators( logistic = function(x, k = 1, x0 = 0) 1 / (1 + exp(-k * (x - x0))), threshold = function(x, c) ifelse(x > c, 1, 0) )ops <- define_custom_operators( logistic = function(x, k = 1, x0 = 0) 1 / (1 + exp(-k * (x - x0))), threshold = function(x, c) ifelse(x > c, 1, 0) )
Evaluates whether the sampling frequency is appropriate for capturing the dynamics of the phenomenon.
diagnose_sampling_frequency(Z, t = NULL)diagnose_sampling_frequency(Z, t = NULL)
Z |
Numeric vector of observations. |
t |
Numeric vector of time points. |
List with diagnostic information and recommendations.
Turns the per draw output of analyze_fixed_points on a
posterior-bearing fit into one row per branch, with an interval instead of
a number.
ed_consensus_fixed_points(fps, n_draws, interval = 0.9)ed_consensus_fixed_points(fps, n_draws, interval = 0.9)
fps |
Data frame returned by |
n_draws |
Number of draws that were swept, including those that found
no fixed point at all and therefore contribute no row to |
interval |
Width of the reported posterior interval. |
A posterior sweep produces a set of fixed points per draw, and the size of that set can differ between draws. That is not noise: it is the posterior disagreeing about how many fixed points there are, which near a saddle-node is the interesting fact rather than a nuisance to be averaged away.
The reduction is therefore explicitly conditional. The modal count k
is taken – ties broken toward the smaller count, so the summary never
claims more structure than the posterior supports – the draws attaining
k are kept, their roots are ordered along the line, and the j-th root
of each is one branch. Ordering along the line is well defined wherever the
count is constant, and stays continuous through a branch crossing because
the roots coincide there. prob_n_fixed_points reports the fraction of
draws the summary was computed from, so the interval is never read as
unconditional.
Data frame with one row per branch: fixed_point (the
posterior median), stability, eigenvalue,
branch, fixed_point_lower, fixed_point_upper,
prob_stability, n_fixed_points,
prob_n_fixed_points and n_draws. Zero rows when most draws
found no fixed point; the fraction that did is then in the
"prob_n_fixed_points" attribute.
Returns the step used by the central difference that
analyze_fixed_points() takes to get the eigenvalue of a fixed point.
ed_derivative_step()ed_derivative_step()
This is exported because it is a commitment to whoever consumes
this package, not an implementation detail. A caller that compares its
own analytic derivative against the one returned here has to know how much
the two are allowed to differ, and that separation is derived from this
step: the truncation of a central difference is . Without
the number the caller either skips the comparison or demands exactness,
and both are wrong.
It was a bare literal inside a function body until 0.1.13, and a consumer was reading it by deparsing that function and matching its source. That dependency broke silently in 0.1.12, when the literal moved into a helper: the scrape found nothing and the consumer's contrast lost its tolerance. Exporting it makes changing it an API change, which is where a change like that belongs.
A single positive finite numeric: the step.
ed_derivative_step()ed_derivative_step()
A clean, publication-ready theme for all diagnostic plots.
ed_theme(base_size = 11, base_family = "")ed_theme(base_size = 11, base_family = "")
base_size |
Base font size. |
base_family |
Base font family. |
A ggplot2 theme object.
Estimates the diffusion coefficient g(.) of an SDE directly from the
quadratic variation of increments, without relying on derivative residuals.
For an SDE dZ = f dt + g dW, the realized quadratic variation satisfies
, providing a
model-free estimator of .
estimate_diffusion_qv( Z, t, predictors, data = NULL, method = c("linear", "quadratic", "gam", "constant"), smooth_window = 21L )estimate_diffusion_qv( Z, t, predictors, data = NULL, method = c("linear", "quadratic", "gam", "constant"), smooth_window = 21L )
Z |
Numeric vector of the state variable. |
t |
Numeric vector of time points. |
predictors |
Character vector of predictor names or data frame. |
data |
Data frame containing the predictors (required if predictors is a character vector). |
method |
Modeling method for g^2: "linear", "quadratic", "gam", or "constant". |
smooth_window |
Size of the rolling median window for smoothing raw quadratic variation (must be odd, default 21). |
This method is recommended over residual-based diffusion estimation when TVR is used to compute derivatives. TVR smooths out high-frequency components that carry the diffusion signature, making residual-based estimation unreliable. Quadratic variation bypasses this by working directly with the raw increments.
The raw per-step estimates are noisy (each is
a single chi-squared realization), so a rolling median is applied before
fitting.
An object of class "variance_model" compatible with the SDE pipeline.
Estimates reasonable starting values for nonlinear least squares.
estimate_initial_values( expression, data, response = NULL, derivative_col = NULL, method = c("grid_search", "random", "heuristic"), n_tries = 100 )estimate_initial_values( expression, data, response = NULL, derivative_col = NULL, method = c("grid_search", "random", "heuristic"), n_tries = 100 )
expression |
Character string with the equation expression. |
data |
Data frame with variables. |
response |
Name of the response variable. |
derivative_col |
Alias for response (for compatibility). |
method |
Estimation method: "grid_search", "random", or "heuristic". |
n_tries |
Number of attempts for random method. |
Named list of initial values.
Refines drift and diffusion estimates using iterative Generalized Least Squares, which is more appropriate when heteroscedasticity is substantial.
estimate_sde_iterative( target, predictors, data, initial_drift = NULL, max_iter = 10, tol = 1e-04, diffusion_method = c("quadratic_variation", "residual"), Z = NULL, t = NULL, selection = c("blocked_cv", "deviance") )estimate_sde_iterative( target, predictors, data, initial_drift = NULL, max_iter = 10, tol = 1e-04, diffusion_method = c("quadratic_variation", "residual"), Z = NULL, t = NULL, selection = c("blocked_cv", "deviance") )
target |
Numeric vector of target values (derivatives) |
predictors |
Data frame of predictor variables |
data |
Full data frame |
initial_drift |
Initial drift equation (optional) |
max_iter |
Maximum number of iterations |
tol |
Convergence tolerance: both the relative change of the weighted deviance and the relative distance between successive fitted functions must fall below it. |
diffusion_method |
Method for estimating the final diffusion coefficient:
"quadratic_variation" (default) estimates g(.) directly from increments
(DeltaZ)^2/dt, which is robust when TVR is used for derivatives.
"residual" uses the classical residual-based approach via
|
Z |
Numeric vector of the state variable (required when diffusion_method = "quadratic_variation"). |
t |
Numeric vector of time points (required when diffusion_method = "quadratic_variation"). |
selection |
How the candidate returned is chosen among the iterations of
the trajectory. |
Convergence is declared only when both conditions hold:
the relative change of the weighted deviance and the relative distance
between successive fitted functions, each below tol. Either alone
has a failure mode the other covers – the objective can stop moving while
the search oscillates between two structures of equivalent fit – and a
criterion that cannot be evaluated (non-finite predictions) is reported as
such instead of standing in for a convergence that never arrives.
An object of class sde_model. Besides drift and
diffusion it carries what the loop did, because a convergence that
cannot be audited is a word rather than a fact:
converged, converged_at and stop_reason
("converged", "max_iter" or "non_finite") –
recorded rather than inferred from the iteration count.
history: one entry per iteration with its equation, the
weights it was fitted under, its fitted values, its weighted deviance
and what each convergence condition was worth.
selected_iteration and selection_scores: which
iteration was returned and the score of every one of them, all
re-computed under the FINAL weights. Iterations cannot be compared
under their own weights, because those weights change every round.
selection_is_in_sample: always TRUE, and said out
loud. Every candidate was discovered while looking at all of these
data, so the winner's score is optimistic; what this supports is the
ranking between candidates, not an out-of-sample error.
Functions for visually exploring the structure of dynamical systems before formal model fitting. These diagnostics should be used BEFORE any symbolic search to inform hypotheses about functional forms.
Generates a battery of diagnostic plots to explore the dynamical structure of the data and suggests potential functional forms.
explore_dynamics( data, target, predictors = NULL, time = NULL, n_bins = 10, include = "all" )explore_dynamics( data, target, predictors = NULL, time = NULL, n_bins = 10, include = "all" )
data |
Data frame containing the time series. |
target |
Name of the target variable (or its derivative). |
predictors |
Character vector of predictor variable names. |
time |
Name of the time column (auto-detected if NULL). |
n_bins |
Number of bins for conditional analysis. |
include |
Which blocks to run: |
A list containing:
suggestions: Character vector of suggested functional forms
statistics: A named list of diagnostic statistics, not a data
frame. Its entries are target_mean, target_sd,
target_range, trend_correlation, and then, for each
predictor, cor_<pred> with the correlation, nonlin_<pred>
with the name of the form the AIC preferred, fit_<pred> with that
form itself, and interaction_<a>_<b> with the p-value of a pair.
plots: List of ggplot objects (if available)
Every pair of predictors that is tested is recorded
in statistics under interaction_<a>_<b>, whether or not it
came out significant, so that a pair which was tested and came out flat
can be told from a pair that was never tested. The threshold used to
announce a pair in suggestions is 0.05 and is fixed, and no
correction for multiplicity is applied: with p predictors there are
p(p-1)/2 tests, so at six predictors the probability of at least one
spurious announcement under the null is about 0.54. The p-values are all
returned precisely so that a caller who needs a correction can apply the
one they can defend.
fit_<pred> carries the fit that won the
AIC comparison, so that a caller can ask questions of the curve and not
only of its name. It is a list with form (the same word as
nonlin_<pred>), degree (1, 2 or 3), coefficients,
predictor_range and aic (all three compared values).
The coefficients are in ascending powers of the predictor and are
degree + 1 long, so that the fitted value at x is
sum(coefficients * x^(0:degree)). An aliased term comes back as
NA in its own slot rather than shifting the powers. The range is the
observed range of the predictor, which is the interval over which anything
asked of the fitted curve is a statement about the data rather than an
extrapolation.
# Toy example data <- data.frame( time = 1:50, Z = sin(seq(0, 10, length.out = 50)), X = cos(seq(0, 10, length.out = 50)) ) data$dZ <- c(diff(data$Z)/diff(data$time), NA) data <- na.omit(data) result <- explore_dynamics(data, target = "dZ", predictors = c("Z", "X") ) print(result$suggestions)# Toy example data <- data.frame( time = 1:50, Z = sin(seq(0, 10, length.out = 50)), X = cos(seq(0, 10, length.out = 50)) ) data$dZ <- c(diff(data$Z)/diff(data$time), NA) data <- na.omit(data) result <- explore_dynamics(data, target = "dZ", predictors = c("Z", "X") ) print(result$suggestions)
Exports analysis results to various file formats.
export_results( results, output_dir, prefix = "empirical_dynamics", formats = c("rds", "csv") )export_results( results, output_dir, prefix = "empirical_dynamics", formats = c("rds", "csv") )
results |
Analysis results list. |
output_dir |
Output directory (required, no default to comply with CRAN policy). |
prefix |
File name prefix. |
formats |
Vector of formats: "rds", "csv", "json", "latex". |
List of paths to created files.
# Toy example tmp_dir <- tempdir() mock_results <- list( equation = stats::lm(mpg ~ wt, data = mtcars) ) # Export paths <- export_results(mock_results, output_dir = tmp_dir, formats = c("csv", "rds"))# Toy example tmp_dir <- tempdir() mock_results <- list( equation = stats::lm(mpg ~ wt, data = mtcars) ) # Export paths <- export_results(mock_results, output_dir = tmp_dir, formats = c("csv", "rds"))
Fits candidate probability distributions to residuals, optionally with parameters that depend on state variables.
fit_residual_distribution( residuals, candidates = c("normal", "t", "skew-normal"), conditional_on = NULL, data = NULL )fit_residual_distribution( residuals, candidates = c("normal", "t", "skew-normal"), conditional_on = NULL, data = NULL )
residuals |
Numeric vector of residuals |
candidates |
Character vector of distribution families to try |
conditional_on |
Formula for conditional parameters (optional) |
data |
Data frame (required if conditional_on specified) |
List with best fitting distribution and parameters
Fits a researcher-specified functional form, estimating only the parameters. Uses Levenberg-Marquardt algorithm for robustness.
fit_specified_equation( expression, data, response = NULL, derivative_col = NULL, start = NULL, method = c("LM", "nls", "optim"), weights = NULL, lower = -Inf, upper = Inf )fit_specified_equation( expression, data, response = NULL, derivative_col = NULL, start = NULL, method = c("LM", "nls", "optim"), weights = NULL, lower = -Inf, upper = Inf )
expression |
Character string specifying the equation (e.g., "a + b * Z"). |
data |
Data frame with predictor variables. |
response |
Name of the response/target column. |
derivative_col |
Alias for response (for compatibility). |
start |
List of starting values for parameters (auto-estimated if NULL). |
method |
Optimization method: "LM" (Levenberg-Marquardt, recommended), "nls" (standard), or "optim" (general optimization). |
weights |
Optional weight vector. |
lower |
Lower bounds for parameters (for "optim" method). |
upper |
Upper bounds for parameters (for "optim" method). |
An object of class "symbolic_equation" containing the fitted model.
# Toy example data <- data.frame(Z = seq(1, 10, length.out = 20)) data$dZ <- 0.5 * data$Z * (1 - data$Z / 20) + rnorm(20, sd = 0.01) # Fit logistic equation eq <- fit_specified_equation( expression = "r * Z * (1 - Z/K)", data = data, response = "dZ", start = list(r = 0.5, K = 20) ) print(eq)# Toy example data <- data.frame(Z = seq(1, 10, length.out = 20)) data$dZ <- 0.5 * data$Z * (1 - data$Z / 20) + rnorm(20, sd = 0.01) # Fit logistic equation eq <- fit_specified_equation( expression = "r * Z * (1 - Z/K)", data = data, response = "dZ", start = list(r = 0.5, K = 20) ) print(eq)
Creates a nicely formatted string representation of the equation.
format_equation( equation, format = c("text", "latex", "markdown"), precision = 4 )format_equation( equation, format = c("text", "latex", "markdown"), precision = 4 )
equation |
Fitted equation object. |
format |
Output format: "text", "latex", "markdown". |
precision |
Number of decimal places. |
Formatted string.
Creates a comprehensive report of the entire analysis workflow.
generate_report( results, output_file, format = c("markdown", "html", "latex"), title = "Empirical Dynamics Analysis Report", author = "EmpiricalDynamics", include_plots = TRUE )generate_report( results, output_file, format = c("markdown", "html", "latex"), title = "Empirical Dynamics Analysis Report", author = "EmpiricalDynamics", include_plots = TRUE )
results |
List containing analysis results with elements:
|
output_file |
Path for output file (required, no default to comply with CRAN policy). |
format |
Report format: "markdown", "html", "latex". |
title |
Report title. |
author |
Author name. |
include_plots |
Include diagnostic plots. |
Path to generated report.
# Toy example to demonstrate report generation # Using a temporary file to avoid writing to user's working directory tmp_file <- tempfile("report_example") # Mock results object mock_results <- list( data = data.frame(time = 1:10, Z = runif(10)), equation = stats::lm(Z ~ time, data = data.frame(time = 1:10, Z = runif(10))) ) # Generate report report_path <- generate_report(mock_results, output_file = tmp_file, format = "markdown") if(file.exists(report_path)) unlink(report_path)# Toy example to demonstrate report generation # Using a temporary file to avoid writing to user's working directory tmp_file <- tempfile("report_example") # Mock results object mock_results <- list( data = data.frame(time = 1:10, Z = runif(10)), equation = stats::lm(Z ~ time, data = data.frame(time = 1:10, Z = runif(10))) ) # Generate report report_path <- generate_report(mock_results, output_file = tmp_file, format = "markdown") if(file.exists(report_path)) unlink(report_path)
Returns a template script for running a complete analysis.
get_analysis_template(output_file = NULL)get_analysis_template(output_file = NULL)
output_file |
Path to save template. |
Template code as character string (invisibly).
# Save template to a temporary file tmp_file <- tempfile("analysis_template", fileext = ".R") get_analysis_template(tmp_file) # Clean up if (file.exists(tmp_file)) unlink(tmp_file)# Save template to a temporary file tmp_file <- tempfile("analysis_template", fileext = ".R") get_analysis_template(tmp_file) # Clean up if (file.exists(tmp_file)) unlink(tmp_file)
Returns all equations on the Pareto front as a list.
get_pareto_set(results)get_pareto_set(results)
results |
A symbolic_search_result object. |
List of symbolic_equation objects.
Returns information about the example datasets included with the package.
list_example_data()list_example_data()
A data.frame with dataset names and descriptions.
list_example_data()list_example_data()
Load one of the example datasets included with the package.
load_example_data(name)load_example_data(name)
name |
Name of the dataset to load. Available datasets:
|
A data.frame containing the time series data.
# Load logistic growth data if available if(requireNamespace("utils", quietly = TRUE)) { try({ data <- load_example_data("logistic_growth") head(data) }) }# Load logistic growth data if available if(requireNamespace("utils", quietly = TRUE)) { try({ data <- load_example_data("logistic_growth") head(data) }) }
Creates a table comparing multiple candidate equations.
model_comparison_table( equations, data, derivative_col, format = c("data.frame", "latex", "markdown"), caption = "Model Comparison" )model_comparison_table( equations, data, derivative_col, format = c("data.frame", "latex", "markdown"), caption = "Model Comparison" )
equations |
Named list of fitted equation objects. |
data |
Data for computing fit statistics. |
derivative_col |
Response variable column. |
format |
Output format. |
caption |
Table caption. |
Comparison table.
Estimates how the residual variance depends on state variables, used for constructing the diffusion term of an SDE.
model_conditional_variance( residuals, predictors, data = NULL, method = c("symbolic", "linear", "quadratic", "gam", "constant"), transform = c("absolute", "squared", "log_squared"), ... )model_conditional_variance( residuals, predictors, data = NULL, method = c("symbolic", "linear", "quadratic", "gam", "constant"), transform = c("absolute", "squared", "log_squared"), ... )
residuals |
Numeric vector of residuals |
predictors |
Formula or data frame of predictor variables (or vector of names) |
data |
Data frame (if predictors is a formula or vector of names) |
method |
Modeling method: "symbolic", "linear", "quadratic", "gam", or "constant" |
transform |
Transformation of residuals: "squared", "absolute", or "log_squared" |
... |
Additional arguments passed to the modeling function |
An object of class "variance_model" containing the fitted model
Functions for generating publication-ready outputs including LaTeX equations, comprehensive reports, and formatted summaries.
A symbolic search returns its equations with the constants already fitted, as
numeric literals: "(x * x) * 3.9305". Anything that has to re-estimate
the equation on a subsample – cross-validation fold by fold, a bootstrap
replicate – needs those constants back as free parameters, together with
starting values. This function performs that conversion: every numeric
literal becomes a named parameter and the literal itself becomes its starting
value, which is the best possible warm start because it is the optimum the
search already found on the full sample.
parameterize_equation( expression, data = NULL, prefix = "c", protect_exponents = TRUE )parameterize_equation( expression, data = NULL, prefix = "c", protect_exponents = TRUE )
expression |
Character string with the equation, or a
|
data |
Optional data frame. When supplied, its column names are avoided when naming the new parameters, so that a parameter can never collide with a variable. |
prefix |
Prefix for the generated parameter names. |
protect_exponents |
Logical. When |
The expression is taken apart with R's own parser and rebuilt from its
syntax tree, never by textual substitution: a regular expression cannot tell
the 3 of a coefficient from the 3 inside a variable named
Z3.
A list with
expression |
The parameterised expression, as a character string. |
start |
Named list of starting values, one per generated parameter. |
n_parameters |
Number of parameters generated. |
An expression carrying no numeric literal is returned unchanged, with an
empty start and n_parameters = 0; that is not an error, it
simply means there is nothing to re-estimate.
fit_specified_equation, cross_validate
# A constant already fitted by the search becomes a free parameter parameterize_equation("(x * x) * 3.9305288261198923") # Digits inside a variable name are not constants parameterize_equation("(4.07 * sinX) - Z3") # Nothing to estimate: returned unchanged parameterize_equation("YmX")# A constant already fitted by the search becomes a free parameter parameterize_equation("(x * x) * 3.9305288261198923") # Digits inside a variable name are not constants parameterize_equation("(4.07 * sinX) - Z3") # Nothing to estimate: returned unchanged parameterize_equation("YmX")
Creates a scatter plot with optional nonparametric fit and marginal distributions.
plot_bivariate( data, x_var, y_var, color_var = NULL, show_fit = TRUE, show_marginals = FALSE )plot_bivariate( data, x_var, y_var, color_var = NULL, show_fit = TRUE, show_marginals = FALSE )
data |
Data frame. |
x_var |
X variable name. |
y_var |
Y variable name. |
color_var |
Optional variable for color mapping. |
show_fit |
Add smooth fit? |
show_marginals |
Add marginal histograms? |
A ggplot object.
Visualizes the trade-off between equation complexity and fit quality.
plot_pareto_front(results, highlight_selection = "knee", show_all = FALSE)plot_pareto_front(results, highlight_selection = "knee", show_all = FALSE)
results |
A symbolic_search_result object. |
highlight_selection |
Which equation to highlight ("knee", "BIC", "AIC", or index). |
show_all |
Show all equations or just Pareto front? |
A ggplot object.
Creates a phase diagram plotting dZ vs Z, useful for visualizing autonomous dynamics and identifying fixed points.
plot_phase_1d( data, z_var, dz_var = NULL, show_zero_line = TRUE, show_fit = TRUE, fit_method = "loess" )plot_phase_1d( data, z_var, dz_var = NULL, show_zero_line = TRUE, show_fit = TRUE, fit_method = "loess" )
data |
Data frame. |
z_var |
Name of state variable Z. |
dz_var |
Name of derivative dZ (auto-constructed if starts with "d_"). |
show_zero_line |
Add horizontal line at dZ = 0? |
show_fit |
Add nonparametric fit? |
fit_method |
Method for fit: "loess", "gam", or "spline". |
In the phase diagram:
Points where the curve crosses dZ = 0 are fixed points
Negative slope at crossing indicates stability
Positive slope indicates instability
A ggplot object.
Creates a multi-panel diagnostic plot for residual analysis.
plot_residual_diagnostics_panel(x, ...)plot_residual_diagnostics_panel(x, ...)
x |
Object of class residual_diagnostics |
... |
Additional arguments passed to plotting functions |
Invisibly returns the input object
Creates a 3D surface or contour plot showing how the target variable depends on two predictors.
plot_surface_3d( data, x_var, y_var, z_var, type = c("contour", "filled_contour", "persp"), n_grid = 30, method = "loess" )plot_surface_3d( data, x_var, y_var, z_var, type = c("contour", "filled_contour", "persp"), n_grid = 30, method = "loess" )
data |
Data frame. |
x_var |
First predictor variable. |
y_var |
Second predictor variable. |
z_var |
Response variable (target). |
type |
Plot type: "contour", "filled_contour", or "persp". |
n_grid |
Grid resolution for surface estimation. |
method |
Surface fitting method: "loess", "gam", or "linear". |
A plot (base graphics for persp, ggplot for contour).
Creates a time series plot with optional trend line and change point detection.
plot_timeseries( data, var, time = NULL, show_trend = TRUE, highlight_changes = TRUE )plot_timeseries( data, var, time = NULL, show_trend = TRUE, highlight_changes = TRUE )
data |
Data frame. |
var |
Variable name to plot. |
time |
Time variable name. |
show_trend |
Add trend line? |
highlight_changes |
Highlight potential structural breaks? |
A ggplot object.
Plots the trajectory of a system in the (Z, X) plane, useful for visualizing attractors and limit cycles.
plot_trajectory_2d( data, x_var, y_var, time_var = NULL, show_arrows = TRUE, arrow_spacing = 10, show_start_end = TRUE )plot_trajectory_2d( data, x_var, y_var, time_var = NULL, show_arrows = TRUE, arrow_spacing = 10, show_start_end = TRUE )
data |
Data frame. |
x_var |
First state variable. |
y_var |
Second state variable. |
time_var |
Time variable (for coloring trajectory). |
show_arrows |
Add direction arrows? |
arrow_spacing |
Spacing between arrows (every nth point). |
show_start_end |
Mark start and end points? |
A ggplot object.
Produces a four-panel diagnostic plot showing the original series, estimated derivative, reconstruction comparison, and residuals.
plot_tvr_diagnostic(Z, t = NULL, dZ_tvr)plot_tvr_diagnostic(Z, t = NULL, dZ_tvr)
Z |
Numeric vector of original observations. |
t |
Numeric vector of time points. |
dZ_tvr |
TVR derivative object (from |
Invisibly returns a list of diagnostic values.
Plot Bifurcation Diagram
## S3 method for class 'bifurcation_analysis' plot(x, ...)## S3 method for class 'bifurcation_analysis' plot(x, ...)
x |
Object of class bifurcation_analysis. |
... |
Additional arguments (ignored). |
A ggplot object.
Plot CV Results
## S3 method for class 'cv_result' plot(x, type = c("predictions", "folds", "both"), ...)## S3 method for class 'cv_result' plot(x, type = c("predictions", "folds", "both"), ...)
x |
Object of class cv_result. |
type |
Type of plot: "predictions", "folds", or "both". |
... |
Additional arguments (ignored). |
A ggplot object or a list of ggplot objects.
Plot Simulated Trajectories
## S3 method for class 'trajectory_simulation' plot( x, observed_data = NULL, show_trajectories = TRUE, n_show = 20, alpha_traj = 0.2, ... )## S3 method for class 'trajectory_simulation' plot( x, observed_data = NULL, show_trajectories = TRUE, n_show = 20, alpha_traj = 0.2, ... )
x |
Object of class trajectory_simulation. |
observed_data |
Optional observed data to overlay. |
show_trajectories |
Show individual trajectories? |
n_show |
Number of trajectories to show. |
alpha_traj |
Transparency for trajectories. |
... |
Additional arguments (ignored). |
A ggplot object.
Plot Method for TVR Derivative
## S3 method for class 'tvr_derivative' plot(x, t = NULL, ...)## S3 method for class 'tvr_derivative' plot(x, t = NULL, ...)
x |
A tvr_derivative object. |
t |
Time vector (uses attribute if NULL). |
... |
Additional plot arguments. |
Invisibly returns the input object (called for side effects).
Plot Validation Results
## S3 method for class 'validation_result' plot(x, ...)## S3 method for class 'validation_result' plot(x, ...)
x |
Object of class validation_result. |
... |
Additional arguments (ignored). |
A list of ggplot objects (invisible).
Predict from Variance Model
## S3 method for class 'variance_model' predict(object, newdata, ...)## S3 method for class 'variance_model' predict(object, newdata, ...)
object |
Variance model object |
newdata |
New data for prediction |
... |
Additional arguments |
Numeric vector of predicted standard deviations.
Functions for data preparation, variable specification, and numerical differentiation including Total Variation Regularized (TVR) differentiation for noisy economic data.
Print Analysis Summary
print_summary(results)print_summary(results)
results |
Analysis results list. |
Invisibly returns the input results object (called for side effects).
Print a Bifurcation Analysis
## S3 method for class 'bifurcation_analysis' print(x, ...)## S3 method for class 'bifurcation_analysis' print(x, ...)
x |
Object of class bifurcation_analysis. |
... |
Additional arguments (ignored). |
x, invisibly.
Print CV Results
## S3 method for class 'cv_result' print(x, ...)## S3 method for class 'cv_result' print(x, ...)
x |
Object of class cv_result. |
... |
Additional arguments (ignored). |
Invisibly returns the input object (called for side effects).
Print Qualitative Check Results
## S3 method for class 'qualitative_check' print(x, ...)## S3 method for class 'qualitative_check' print(x, ...)
x |
Object of class qualitative_check. |
... |
Additional arguments (ignored). |
Invisibly returns the input object (called for side effects).
Print Residual Diagnostics
## S3 method for class 'residual_diagnostics' print(x, ...)## S3 method for class 'residual_diagnostics' print(x, ...)
x |
Object of class residual_diagnostics |
... |
Additional arguments (ignored) |
Invisibly returns the input object (called for side effects).
Print Method for TVR Derivative
## S3 method for class 'tvr_derivative' print(x, ...)## S3 method for class 'tvr_derivative' print(x, ...)
x |
A tvr_derivative object. |
... |
Additional arguments (ignored). |
Invisibly returns the input object (called for side effects).
Print Validation Results
## S3 method for class 'validation_result' print(x, ...)## S3 method for class 'validation_result' print(x, ...)
x |
Object of class validation_result. |
... |
Additional arguments (ignored). |
Invisibly returns the input object (called for side effects).
Reads time series or panel data from various file formats and prepares it for use with EmpiricalDynamics functions.
read_empirical_data(file, time_col = NULL, date_format = NULL, ...)read_empirical_data(file, time_col = NULL, date_format = NULL, ...)
file |
Path to the data file (CSV, RDS, or RData). |
time_col |
Name of the time/date column (auto-detected if NULL). |
date_format |
Date format string if time column is character. |
... |
Additional arguments passed to read.csv. |
A data.frame with time column converted to numeric if needed.
Functions for analyzing residual structure, modeling conditional variance, constructing stochastic differential equations (SDEs), and iterative GLS estimation for heteroscedastic systems.
Performs a battery of statistical tests on model residuals to check for autocorrelation, heteroscedasticity, and normality.
residual_diagnostics( residuals, data = NULL, predictors = NULL, max_lag = 10, plot = TRUE )residual_diagnostics( residuals, data = NULL, predictors = NULL, max_lag = 10, plot = TRUE )
residuals |
Numeric vector of residuals (or model object) |
data |
Optional data frame for conditional tests |
predictors |
Variable names for heteroscedasticity tests |
max_lag |
Maximum lag for autocorrelation tests |
plot |
Produce diagnostic plots? |
A list of test results with class "residual_diagnostics"
Saves all diagnostic plots to files.
save_plots( results, output_dir, prefix = "empirical_dynamics", format = c("png", "pdf"), width = 8, height = 6, dpi = 300 )save_plots( results, output_dir, prefix = "empirical_dynamics", format = c("png", "pdf"), width = 8, height = 6, dpi = 300 )
results |
Analysis results list. |
output_dir |
Output directory (required, no default to comply with CRAN policy). |
prefix |
File name prefix. |
format |
Image format: "png", "pdf", "svg". |
width |
Plot width in inches. |
height |
Plot height in inches. |
dpi |
Resolution for raster formats. |
List of paths to created files.
Selects the best equation from the Pareto front using a specified criterion.
select_equation( results, criterion = c("knee", "BIC", "AIC", "min_complexity", "min_error"), n = NULL )select_equation( results, criterion = c("knee", "BIC", "AIC", "min_complexity", "min_error"), n = NULL )
results |
A symbolic_search_result object. |
criterion |
Selection criterion: "knee", "BIC", "AIC", "min_complexity", or "min_error". |
n |
Sample size (required for BIC/AIC if not stored). |
A symbolic_equation object.
Selects the regularization parameter lambda using leave-one-out-like cross-validation. Problems are internally rescaled for numerical stability, and a solver chain is used for each candidate lambda.
select_lambda_cv_tvr( Z, t, A = NULL, D = NULL, solver = "clarabel", lambda_seq = 10^seq(-4, 2, length.out = 30), verbose = FALSE )select_lambda_cv_tvr( Z, t, A = NULL, D = NULL, solver = "clarabel", lambda_seq = 10^seq(-4, 2, length.out = 30), verbose = FALSE )
Z |
Numeric vector of observations. |
t |
Numeric vector of time points. |
A |
Integration matrix. |
D |
Difference matrix. |
solver |
Preferred solver (default "clarabel"). Falls back to other solvers if the preferred one fails. |
lambda_seq |
Sequence of lambda values to evaluate. |
verbose |
Print progress? |
Selected lambda value.
Examines how sensitive the model predictions are to parameter perturbations.
sensitivity_analysis( equation, data, response = NULL, derivative_col = NULL, perturbation_pct = 10, n_bootstrap = 100 )sensitivity_analysis( equation, data, response = NULL, derivative_col = NULL, perturbation_pct = 10, n_bootstrap = 100 )
equation |
Fitted equation object. |
data |
Data for evaluation. |
response |
Name of response column. |
derivative_col |
Alias for response. |
perturbation_pct |
Percentage perturbation (default 10%). |
n_bootstrap |
Number of bootstrap samples for uncertainty. |
Data frame with sensitivity metrics for each parameter.
Checks if Julia and the required SymbolicRegression.jl package are installed.
setup_julia_backend()setup_julia_backend()
Logical indicating if the backend is ready.
Simulates trajectories using the discovered SDE to assess whether the model can reproduce observed dynamics.
simulate_trajectory( sde, initial_conditions, times, n_sims = 100, method = c("euler", "milstein", "rk4"), exogenous_data = NULL, seed = NULL )simulate_trajectory( sde, initial_conditions, times, n_sims = 100, method = c("euler", "milstein", "rk4"), exogenous_data = NULL, seed = NULL )
sde |
SDE object from |
initial_conditions |
Named vector of initial values for all variables. |
times |
Numeric vector of time points. |
n_sims |
Number of Monte Carlo simulations (for stochastic models). |
method |
Integration method: "euler", "milstein", "rk4" (deterministic only). |
exogenous_data |
Data frame with exogenous variable trajectories (if any). |
seed |
Random seed for reproducibility. |
Object of class "trajectory_simulation" containing:
trajectories |
Array of simulated trajectories (time x variable x simulation) |
times |
Time points |
summary |
Summary statistics (mean, quantiles) at each time |
# Toy example: dX = 0.5 * X # Mock SDE object structure sde <- list( drift = list(expression = "0.5 * X"), diffusion = list(expression = "0.1"), # Add noise variable = "X" ) class(sde) <- "sde_model" # Simulation sim <- simulate_trajectory( sde = sde, initial_conditions = c(X = 1), times = seq(0, 1, by = 0.1), n_sims = 10, seed = 123 ) print(sim$summary$mean)# Toy example: dX = 0.5 * X # Mock SDE object structure sde <- list( drift = list(expression = "0.5 * X"), diffusion = list(expression = "0.1"), # Add noise variable = "X" ) class(sde) <- "sde_model" # Simulation sim <- simulate_trajectory( sde = sde, initial_conditions = c(X = 1), times = seq(0, 1, by = 0.1), n_sims = 10, seed = 123 ) print(sim$summary$mean)
Classifies variables in a dataset according to their role in the dynamical system being studied.
specify_variables( data, endogenous = NULL, endogenous_coupled = NULL, coupled = NULL, exogenous = NULL, slow_parameter = NULL, time = NULL, time_col = NULL )specify_variables( data, endogenous = NULL, endogenous_coupled = NULL, coupled = NULL, exogenous = NULL, slow_parameter = NULL, time = NULL, time_col = NULL )
data |
A data.frame containing the time series data. |
endogenous |
Character vector of endogenous state variable names. |
endogenous_coupled |
Character vector of coupled endogenous variables. |
coupled |
Alias for endogenous_coupled (for compatibility). |
exogenous |
Character vector of exogenous forcing variable names. |
slow_parameter |
Character vector of slowly-varying parameter names. |
time |
Name of the time column (auto-detected if NULL). |
time_col |
Alias for time (for compatibility). |
Variable types:
endogenous: Variables whose dynamics are modeled (appear as dZ/dt).
endogenous_coupled: Variables that co-evolve with endogenous vars.
exogenous: Variables that influence the system but are not modeled.
slow_parameter: Variables that change on much longer timescales.
The input data.frame with variable specifications added as the "var_spec" attribute.
data <- data.frame( time = 1:10, profit_rate = runif(10), capital_stock = runif(10), interest_rate = runif(10) ) data <- specify_variables(data, endogenous = "profit_rate", endogenous_coupled = "capital_stock", exogenous = "interest_rate" ) attr(data, "var_spec")data <- data.frame( time = 1:10, profit_rate = runif(10), capital_stock = runif(10), interest_rate = runif(10) ) data <- specify_variables(data, endogenous = "profit_rate", endogenous_coupled = "capital_stock", exogenous = "interest_rate" ) attr(data, "var_spec")
Analyzes the time series to recommend the most appropriate differentiation method based on detected features like trend, periodicity, shocks, and noise.
suggest_differentiation_method(Z, t = NULL)suggest_differentiation_method(Z, t = NULL)
Z |
Numeric vector of observations. |
t |
Numeric vector of time points. |
List with suggested method and diagnostic information.
t <- 1:100 Z <- 0.1 * t + rnorm(100) # Trend with noise result <- suggest_differentiation_method(Z, t) print(result$suggested_method)t <- 1:100 Z <- 0.1 * t + rnorm(100) # Trend with noise result <- suggest_differentiation_method(Z, t) print(result$suggested_method)
Functions for discovering functional forms through symbolic regression using genetic algorithms. Interfaces with Julia's SymbolicRegression.jl for advanced search, with fallback to R-native methods for simpler cases.
Discovers the functional form of a differential equation from data using genetic/evolutionary algorithms. Returns a Pareto front of equations trading off complexity against fit.
symbolic_search( target, predictors, operators = NULL, constraints = NULL, n_runs = 5, complexity_penalty = 0.05, parsimony_pressure = c("adaptive", "constant", "none"), backend = c("r_genetic", "julia", "r_exhaustive"), julia_options = NULL, weights = NULL, verbose = TRUE )symbolic_search( target, predictors, operators = NULL, constraints = NULL, n_runs = 5, complexity_penalty = 0.05, parsimony_pressure = c("adaptive", "constant", "none"), backend = c("r_genetic", "julia", "r_exhaustive"), julia_options = NULL, weights = NULL, verbose = TRUE )
target |
Numeric vector of target values (typically derivatives). |
predictors |
Data frame of predictor variables. |
operators |
List specifying allowed operators:
|
constraints |
List of constraints:
|
n_runs |
Number of independent runs for robustness. |
complexity_penalty |
Penalty per unit complexity. |
parsimony_pressure |
Type of parsimony: "constant", "adaptive", or "none". |
backend |
Computation backend: "julia", "r_genetic", or "r_exhaustive". |
julia_options |
List of options passed to SymbolicRegression.jl. |
weights |
Optional weight vector for weighted regression. Honoured by
every backend: the Julia arms pass it to SymbolicRegression.jl through a
|
verbose |
Print progress messages? |
An object of class "symbolic_search_result" containing:
pareto_front: Data frame of Pareto-optimal equations
all_equations: All discovered equations
best_by_complexity: Best equation at each complexity level
run_diagnostics: Information about each run
# Toy example using R-native exhaustive search (fastest for demo) data <- data.frame( x = seq(1, 10, length.out = 20), y = seq(1, 10, length.out = 20)^2 + rnorm(20, sd = 0.1) ) # Discover y ~ x^2 results <- symbolic_search( target = data$y, predictors = data["x"], backend = "r_exhaustive" ) print(head(results$pareto_front))# Toy example using R-native exhaustive search (fastest for demo) data <- data.frame( x = seq(1, 10, length.out = 20), y = seq(1, 10, length.out = 20)^2 + rnorm(20, sd = 0.1) ) # Discover y ~ x^2 results <- symbolic_search( target = data$y, predictors = data["x"], backend = "r_exhaustive" ) print(head(results$pareto_front))
Performs symbolic search with weighted least squares.
symbolic_search_weighted( target, predictors, weights, operators = NULL, constraints = NULL, n_runs = 3, complexity_penalty = 0.05, verbose = TRUE )symbolic_search_weighted( target, predictors, weights, operators = NULL, constraints = NULL, n_runs = 3, complexity_penalty = 0.05, verbose = TRUE )
target |
Numeric vector of target values (typically derivatives). |
predictors |
Data frame of predictor variables. |
weights |
Optional weight vector for weighted regression. Honoured by
every backend: the Julia arms pass it to SymbolicRegression.jl through a
|
operators |
List specifying allowed operators:
|
constraints |
List of constraints:
|
n_runs |
Number of independent runs for robustness. |
complexity_penalty |
Penalty per unit complexity. |
verbose |
Print progress messages? |
A symbolic_search_result object
Converts a discovered equation to LaTeX format for publication.
to_latex( equation, variable = "Z", precision = 3, scientific_notation = TRUE, include_uncertainty = FALSE, se_values = NULL )to_latex( equation, variable = "Z", precision = 3, scientific_notation = TRUE, include_uncertainty = FALSE, se_values = NULL )
equation |
Fitted equation object. |
variable |
Name of the dependent variable (for dZ/dt notation). |
precision |
Number of decimal places for coefficients. |
scientific_notation |
Use scientific notation for large/small coefficients. |
include_uncertainty |
Include standard errors in parentheses. |
se_values |
Named vector of standard errors (optional). |
Character string with LaTeX equation.
# Toy example using a linear model data <- data.frame(Z = 1:10, dZ = 2 * (1:10) + 3) model <- stats::lm(dZ ~ Z, data = data) # Convert to LaTeX latex_eq <- to_latex(model, variable = "Z") cat(latex_eq)# Toy example using a linear model data <- data.frame(Z = 1:10, dZ = 2 * (1:10) + 3) model <- stats::lm(dZ ~ Z, data = data) # Convert to LaTeX latex_eq <- to_latex(model, variable = "Z") cat(latex_eq)
Runs a battery of validation tests on the discovered equation.
validate_model( equation, sde = NULL, data, response = NULL, derivative_col = NULL, variable, time_col = "time", cv_folds = 5, n_sims = 50, expected_features = list(), verbose = TRUE )validate_model( equation, sde = NULL, data, response = NULL, derivative_col = NULL, variable, time_col = "time", cv_folds = 5, n_sims = 50, expected_features = list(), verbose = TRUE )
equation |
Fitted equation object. |
sde |
SDE object (optional, for trajectory validation). |
data |
Original data frame. |
response |
Name of response column. |
derivative_col |
Alias for response (for compatibility). |
variable |
Main variable name. |
time_col |
Time column name. |
cv_folds |
Number of CV folds. |
n_sims |
Number of trajectory simulations. |
expected_features |
Named list of expected qualitative features, as
|
verbose |
Print progress. |
Object of class "validation_result".
Functions for validating discovered differential equations through cross-validation, trajectory simulation, and qualitative behavior analysis.