Optimize a trend filtering model by minimizing Stein's unbiased risk estimate

For every candidate hyperparameter value, compute an unbiased estimate of the trend filtering model's predictive mean-squared error. See the Details section for guidelines on when sure_trendfilter() should be used versus cv_trendfilter().

sure_trendfilter(x, y, weights, nlambda = 250L, ...)

Arguments

x	Vector of observed values for the input variable.
y	Vector of observed values for the output variable.
weights	Weights for the observed outputs, defined as the reciprocal variance of the additive noise that contaminates the output signal. When the noise is expected to have an equal variance $\sigma^2$ for all observations, 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`.
nlambda	Number of hyperparameter values to evaluate during cross validation. Defaults to `nlambda = 250`. The hyperparameter grid is internally constructed to span the full trend filtering model space (which is bookended by a global polynomial solution and an interpolating solution), with `nlambda` controlling the granularity of the hyperparameter grid.
...	Additional named arguments to pass to `.trendfilter()`.

Value

An object of class 'sure_trendfilter'. The object has subclass 'trendfilter' and may therefore be passed to generic stats functions such as predict(), fitted(), and residuals(). More precisely, a sure_trendfilter' object is a list with the elements below, as well as all elements from the trendfilter' call.

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.
error: Vector of mean-squared prediction errors estimated by SURE, for every hyperparameter value in lambda.
se_error: Vector of estimated standard errors for the error.
training_error: The "in-sample" MSE between the observed outputs y and the trend filtering estimate, for every hyperparameter value in lambda.
optimism: SURE-estimated optimisms, i.e. optimism = error - training_error.
lambda_min: Hyperparameter value in lambda that minimizes the SURE validation error curve.
lambda_1se: The largest hyperparameter value in lambda that has a SURE error within one standard error of min(error). We call this the "1-standard-error rule" hyperparameter, and it serves as an Occam's razor-esque heuristic. More precisely, given two models with approximately equal performance (here, in terms of predictive MSE), it may be wise to opt for the simpler model, i.e. the model with the larger hyperparameter value / fewer effective degrees of freedom.
edf_min: Number of effective degrees of freedom in the trend filtering estimator with hyperparameter lambda_min.
edf_1se: Number of effective degrees of freedom in the trend filtering estimator with hyperparameter lambda_1se.
i_min: Index of lambda that gives lambda_min.
i_1se: Index of lambda that gives lambda_1se.
fitted_values: The fitted values of all trend filtering estimates, return as 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 candidate hyperparameter in lambda.
n_iter: Total number of iterations taken by the ADMM algorithm, for every candidate hyperparameter 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 by passing a max_iter argument to sure_trendfilter() in order to ensure that the ADMM solution has converged to satisfactory precision.
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 estimates.
status: For internal use. Output from the C solver.
call: The function call.
scale: For internal use.

Details

Our recommendations for when to use sure_trendfilter() versus cv_trendfilter() are summarized in the table below. See Section 3.5 of Politsch et al. (2020a) for more details.

Scenario	Hyperparameter optimization
`x` is unevenly sampled	`cv_trendfilter()`
`x` is evenly sampled and measurement variances for `y` are not available	`cv_trendfilter()`
`x` is evenly sampled and measurement variances for `y` are available	`sure_trendfilter()`

For our purposes, an evenly sampled data set with some discarded pixels (either sporadically or in wide consecutive chunks) is still considered to be evenly sampled. When x is evenly sampled on a transformed scale, we recommend transforming to that scale and carrying out the full trend filtering analysis on that scale. See the sure_trendfilter() examples for a case when the inputs are evenly sampled on the log10(x) scale.

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("quasar_spectrum")
head(quasar_spectrum)
#> # A tibble: 6 × 3
#>   log10_wavelength   flux weights
#>              <dbl>  <dbl>   <dbl>
#> 1             3.55  0.424  0.0417
#> 2             3.55 -2.11   0.125
#> 3             3.55 -3.78   0.128
#> 4             3.55  3.60   0.121
#> 5             3.55 -4.89   0.134
#> 6             3.55  1.76   0.132 

x <- quasar_spectrum$log10_wavelength
y <- quasar_spectrum$flux
weights <- quasar_spectrum$weights

sure_tf <- sure_trendfilter(x, y, weights)