inst/doc/L0TFinv-vignette.R

In L0TFinv: A Splicing Approach to the Inverse Problem of L0 Trend Filtering

## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(comment = "#>", warning=FALSE, message=FALSE)

## ----echo = FALSE-------------------------------------------------------------
# Thanks to Yihui Xie for providing this code
library(knitr)
hook_output <- knit_hooks$get("output")
knit_hooks$set(output = function(x, options) {
   lines <- options$output.lines
   if (is.null(lines)) {
     return(hook_output(x, options))  # pass to default hook
   }
   x <- unlist(strsplit(x, "\n"))
   more <- "..."
   if (length(lines)==1) {        # first n lines
     if (length(x) > lines) {
       # truncate the output, but add ....
       x <- c(head(x, lines), more)
     }
   } else {
     x <- c(more, x[lines], more)
   }
   # paste these lines together
   x <- paste(c(x, ""), collapse = "\n")
   hook_output(x, options)
 })

## ----eval=FALSE---------------------------------------------------------------
#  install.packages("L0TFinv_1.0.1.tar.gz", repos = NULL, type = "source")

## ----eval=FALSE---------------------------------------------------------------
#  if(!require(ggplot2)) install.packages("ggplot2")
#  if(!require(Matrix)) install.packages("Matrix")
#  if(!require(stats)) install.packages("stats")

## ----results="hide"-----------------------------------------------------------
library(L0TFinv)

## -----------------------------------------------------------------------------
n = 10
q = 1
D <- DiffMat(n, q)
X <- XMat(n, q)

## -----------------------------------------------------------------------------
print(D)
print(X)
print(D%*%X)

## -----------------------------------------------------------------------------
tau = c(0.2, 0.3, 0.5, 0.65, 0.85)
h = c(-1,3,-2,0,4,-3)
BlocksData <- SimuBlocksInv(n = 500, sigma = 0.1, seed = 50, tau = tau ,h = h)
plot(BlocksData$x, BlocksData$y, xlab="", ylab="") ## The piecewise linear simulated data
lines(BlocksData$x, BlocksData$y0, col = "red") ## The underlying trend
print(BlocksData$setA) ## The set of position indicators of change points
print(BlocksData$tau)

## -----------------------------------------------------------------------------
tau1 = c(0.1, 0.3, 0.4, 0.7, 0.9)
h1 = c(-2, 5, -3, 2, -1, 4)
a0 = -10
WaveData <- SimuWaveInv(n = 500, sigma = 0.1, seed = 50, tau = tau1, h = h1, a0 = a0)
plot(WaveData$x, WaveData$y, xlab="", ylab="")
lines(WaveData$x, WaveData$y0, col = "red")
print(WaveData$setA)
print(WaveData$tau)

## -----------------------------------------------------------------------------
FitBlocks.fix <- L0TFinv.fix(y=BlocksData$y, k=10, q=0, first=0.01, last=1)
FitWave.fix <- L0TFinv.fix(y=WaveData$y, k=8, q=1, first=0, last=0.99)

## -----------------------------------------------------------------------------
FitBlocks.opt <- L0TFinv.opt(y=BlocksData$y, kmax=20, q=0, first=0.01, last=1, penalty="sic")
FitWave.opt <- L0TFinv.opt(y=WaveData$y, kmax=15, q=1, first=0, last=0.99, penalty="bic")

## -----------------------------------------------------------------------------
coef(FitBlocks.fix, k=6)
coef(FitBlocks.opt, k=FitBlocks.opt$kopt)

## ----eval=FALSE---------------------------------------------------------------
#  coef(FitBlocks.fix)
#  coef(FitWave.opt)

## -----------------------------------------------------------------------------
print(FitBlocks.opt)[["mse"]]
print(FitBlocks.opt)[["bic"]]
print(FitBlocks.opt)[["sic"]]

## -----------------------------------------------------------------------------
metrics <- TFmetrics(BlocksData$y0,BlocksData$tau,FitBlocks.opt$yopt,FitBlocks.opt$Aopt/length(BlocksData$y0))
print(metrics)

## -----------------------------------------------------------------------------
plot(FitBlocks.opt,type="yhat")
plot(FitBlocks.opt,type="bic")
plot(FitWave.opt,type="yhat",k=4)
plot(FitWave.opt,type="mse")