Bayesian trend filtering estiamtes the function \(f = f(x_1), \ldots, f(x_n)\) from data \(y_i, i=1, \ldots, n\), hypothesized to be generated from the model \[y_i = f(x_i) + \epsilon_i\] by fitting a fully Bayesian hierarchical model that is equivalent to the minimization problem \[argmin_f ||y - f||_2^2 + \lambda||D^k f||_1\] as detailed by R.J. Tibshirani [-@Tibshirani:2014]. The penalty matrix \(D^k\) estimates $k$th derivatives of \(f\) at various inputs \(x_i\). This is the Bayesian analogue to the function genlasso::trendfilter.
The main function of this package, btf::btf, assumes the inputs \(x_i, i=1, \ldots, n\), are equally spaced and in ascending order. The user needs to specify the order of the derivatives \(k\). I recommend \(k<=3\) since calculations with \(D^k\) can otherwise cause numerical instability. Note that unlike genlasso::trendfilter Bayesian trend filtering will not set any elements of the penalty vector \(D^k f\) exactly equal to zero and thus this method should be used as a smoother.
We'll use the time-series observations datasets::nhtemp. First, fit a model by specifying the observations \(y\), the order of the derivative to use, and possibly the number of iterations.
library(btf)
fit <- btf(nhtemp, k=2, iter=2000)
The variable fit is nothing special, but I recommend retrieving the posterior samples with the functions btf::getPost. You can specify which parameter you want to investigate, by default we assume you are interested in the estimate of \(f\). Output from this function is a coda::mcmc object, suitable for any of the methods specified within the coda package.
s2_sample <- getPost(fit, 's2')
plot(s2_sample)
Alternatively, you can ask for a summary of the parameters requested. The function btf::getPostEst will give the output from coda::summary.mcmc applied to the parameters specified, or \(f\) is no paramters are specified. You can also specify a function, via the named argument est, to be applied to each of the paramters requested
getPostEst(fit, 's2', est=mean)
## var1
## 1.054152
A simple plot of the estimate of \(f\) and highest posterior density intervals is provided by the S3 method btf::plot.btf
plot(fit)
