Fits the non-parametric Gaussian regression model
y = mu + e, where the mean mu is modelled as mu =
Xb, X is a matrix with columns containing an appropriate basis,
and b is vector with a (sparse) SuSiE prior. In particular, when
order = 0, the jth column of X is a vector with the first j
elements equal to zero, and the remaining elements equal to 1, so
that b_j corresponds to the change in the mean of y between
indices j and j+1. For background on trend filtering, see
Tibshirani (2014). See also the "Trend filtering" vignette,
An n-vector of observations ordered in time or space (assumed to be equally spaced).
An integer specifying the order of trend filtering.
Logical indicating whether to standardize the X
variables ("basis functions");
Logical indicating whether to use the "median
absolute deviation" (MAD) method to the estimate residual
Other arguments passed to
This implementation exploits the special structure of X,
which means that the matrix-vector product X^Ty is fast to
compute; in particular, the computation time is O(n) rather
than O(n^2) if
X were formed explicitly. For
implementation details, see the "Implementation of SuSiE trend
filtering" vignette by running
A "susie" fit; see
susie for details.
R. J. Tibshirani (2014). Adaptive piecewise polynomial estimation via trend filtering. Annals of Statistics 42, 285-323.
1 2 3 4 5 6 7 8 9 10 11 12 13 14
set.seed(1) mu = c(rep(0,50),rep(1,50),rep(3,50),rep(-2,50),rep(0,200)) y = mu + rnorm(400) s = susie_trendfilter(y) plot(y) lines(mu,col = 1,lwd = 3) lines(predict(s),col = 2,lwd = 2) # Calculate credible sets (indices of y that occur just before # changepoints). susie_get_cs(s) # Plot with credible sets for changepoints. susie_plot_changepoint(s,y)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.