Fit a trend filtering model (back-end function with more options for experts)

Fit a trend filtering model.

.trendfilter(
  x,
  y,
  weights = NULL,
  lambda,
  edf = NULL,
  k = 2L,
  obj_tol = 1e-10,
  max_iter = length(y),
  ...
)

Arguments

x	Vector of observed values for the input variable.
y	Vector of observed values for the output variable.
weights	Weights for the output measurements. Output weights are defined as the inverse variance of the additive noise that contaminates the output signal. When the noise is expected to have a constant variance $\sigma^2$ over all outputs, a scalar may be passed to `weights`, i.e. `weights =` $1/\sigma^2$. Otherwise, `weights` must be a vector with the same length as `x` and `y`.
lambda	One or more hyperparameter values to fit a trend filtering estimate for.
edf	(Not yet available) Alternative hyperparametrization for the trend filtering model(s). Vector of the desired number of effective degrees of freedom in each model.
k	Degree of the polynomials that make up the piecewise-polynomial trend filtering estimate. Defaults to `k = 2` (i.e. a piecewise quadratic estimate). Must be one of `k = 0,1,2`. Higher order polynomials are disallowed since they yield no statistical benefit over `k = 2` and their use can lead to instability in the convex optimization.
obj_tol	Stopping criterion for the trend filtering convex optimization. If the relative change in the trend filtering objective function between two successive iterations is less than `obj_tol`, the algorithm terminates. Defaults to `obj_tol = 1e-10`.
max_iter	Maximum number of iterations that we will tolerate for the trend filtering convex optimization algorithm. Defaults to `max_iter = length(y)`.
...	Additional named arguments. Currently unused.
scale	A logical indicating whether to scale the inputs and outputs.

Value

An object of class 'trendfilter'. Generic functions such as predict(), fitted(), and residuals() may be called on any object of class (or subclass) 'trendfilter'. A 'trendfilter' object is a list with the following elements:

x: Vector of observed values for the input variable.
y: Vector of observed values for the output variable (if originally present, observations with is.na(y) or weights == 0 are dropped).
weights: Vector of weights for the observed outputs.
k: Degree of the trend filtering estimate.
lambda: Vector of candidate hyperparameter values (always returned in descending order).
edf: Number of effective degrees of freedom in the trend filtering estimator, for every hyperparameter value in lambda.
fitted_values: The fitted values of the trend filtering estimate(s). If length(lambda) == 1, fitted values for the single fit are returned as a numeric vector. Otherwise, fitted values are returned in a matrix with length(lambda) columns, with fitted_values[,i] corresponding to the trend filtering estimate with hyperparameter lambda[i].
admm_params: A list of the parameter values used by the ADMM algorithm used to solve the trend filtering convex optimization.
obj_func: The relative change in the objective function over the ADMM algorithm's final iteration, for every hyperparameter value in lambda.
n_iter: Total number of iterations taken by the ADMM algorithm, for every hyperparameter value in lambda. If an element of n_iter is exactly equal to admm_params$max_iter, then the ADMM algorithm stopped before reaching the objective tolerance admm_params$obj_tol. In these situations, you may need to increase the maximum number of tolerable iterations in order to ensure that the ADMM solution has converged to satisfactory precision. This can be done by passing an extra argument max_iter to the .trendfilter function call and increasing it from its default value max_iter = length(y).
status: For internal use. Output from the C solver.
call: The function call.
scale_xy: For internal use.

References

Politsch et al. (2020a). Trend filtering – I. A modern statistical tool for time-domain astronomy and astronomical spectroscopy. MNRAS, 492(3), p. 4005-4018. [Publisher] [arXiv].
Politsch et al. (2020b). Trend Filtering – II. Denoising astronomical signals with varying degrees of smoothness. MNRAS, 492(3), p. 4019-4032. [Publisher] [arXiv].

Examples

data("eclipsing_binary")
head(eclipsing_binary)
#> # A tibble: 6 × 3
#>    phase  flux std_err
#>    <dbl> <dbl>   <dbl>
#> 1 -0.499 0.938  0.0102
#> 2 -0.498 0.930  0.0102
#> 3 -0.496 0.944  0.0102
#> 4 -0.494 0.963  0.0102
#> 5 -0.494 0.934  0.0102
#> 6 -0.494 0.949  0.0102

x <- eclipsing_binary$phase
y <- eclipsing_binary$flux
weights <- 1 / eclipsing_binary$std_err^2

fit <- .trendfilter(
  x,
  y,
  weights,
  lambda = exp(10),
  obj_tol = 1e-6,
  max_iter = 1e4
)