In natbprice/tvdiff: Numerical differentiation of noisy, nonsmooth data

library(knitr)

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "man/figures/README-",
  out.width = "100%"
)

local({
  hook_plot = knit_hooks$get('plot')
  knit_hooks$set(plot = function(x, options) {
    x = paste(x, collapse = '.')
    if (!grepl('\\.svg', x)) return(hook_plot(x, options))
    # read the content of the svg image and write it out without <?xml ... ?>
    paste0('<img src= "./', x, '">')
  })
})

tvdiff

The tvdiff package is a solution to the problem of how to estimate the rate of change (i.e., derivative) when available data is noisy and nonsmooth. If the data is noisy, then simple finite difference methods of estimating the derivative are often inaccurate. If the data generating process is nonsmooth, then fitting a smooth model to the data in order to estimate the derivative may yield inaccurate estimates.

This package is based on the Matlab implementation of the Total Variation Regularized Numerical Differentiation algorithm by Rick Chartrand. This package uses a C++ method for the preconditioned conjugate method based on the cPCG package. See references section.

Installation

The tvdiff package is currently only available from Github.

# Install development version from GitHub (install function requires devtools package)
devtools::install_github("natbprice/tvdiff")

Installation and use requirements: - R v. 4.0.1+ - C++ compiler - Mac OS users will likely need to install GFortran, unless previously installed. See here for most recent binares (archived versions are unlikely to be compatible). See here for further discussion regarding R compiler tools for R package Rcpp.

Example

A simple example based on the function $f(x) = \mid x - 0.5 \mid$ with Gaussian noise of standard deviation 0.05. The derivative is estimated from the noisy observations using Total Variation Regularized Differentiation. A prediction of the original function is obtained from the estimated derivative through numerical integration.

# load packages
library(tvdiff)

Load demo data

Example data is included in the package for demonstration.

data("smalldemodata")
str(smalldemodata)

Estimate derivative

The TVRegDiffR function is used to estimate the derivative. The main tuning parameter is the regularization parameter alph. Start by varying alph by orders of magnitude until reasonable results are obtained. The number of iterations iter must also be specified since there is no default stopping criteria.

# Unpack data
x <- smalldemodata$x
obs <- smalldemodata$obs
true <- smalldemodata$true

# True derivative
dydx_true <- rep(-1, length(x))
dydx_true[x > 0.5] <- 1
dx <- x[2] - x[1]

# Extimate derivative
dydx <- TVRegDiffR(
  data = obs, # Vector of data to be differentiated
  iter = 1e3, # Number of iterations
  alph = 0.2, # Regularlization parameter
  dx = dx     # Grid spacing 
  )
dydx <- dydx[-1]

# Prediction
pred <- obs[1] + cumsum(dydx*dx)

Plot results

For this simple example we can compare the estimates to the true values.

# Plot observed vs predicted
plot(x, obs, pch = 20, col = "#57575F", ylab = "y")
lines(x, pred, col = "#CA3542", lty = 2, lwd = 2)
lines(x, true, col = "#27647B", lty = 1, lwd = 2)
legend("bottomleft", 
       legend = c("Observed", "Predicted", "True"),
       col = c("#57575F", "#CA3542", "#27647B"),
       pch = c(20, NA, NA),
       lty = c(NA, 2, 1),
       cex = 0.8,
       lwd = c(NA, 2, 2))

# Plot derivative
plot(x, dydx, ylim = c(-1,1), col = "#CA3542", type = "l", lty = 2, lwd = 2,
     ylab = expression(dy/dx))
lines(x, dydx_true, col = "#27647B", lty = 1, lwd = 2)
legend("topleft", 
       legend = c("Predicted", "True"),
       col = c("#CA3542", "#27647B"),
       lty = c(2, 1),
       lwd = c(2, 2),
       cex = 0.8)

References {#references}

Rick Chartrand, “Numerical Differentiation of Noisy, Nonsmooth Data,” ISRN Applied Mathematics, vol. 2011, Article ID 164564, 11 pages, 2011. https://doi.org/10.5402/2011/164564.

Matlab code: https://sites.google.com/site/dnartrahckcir/home/tvdiff-code

Python translation: https://github.com/stur86/tvregdiff

C++ implementation of preconditioned conjugate gradient method: https://github.com/styvon/cPCG

natbprice/tvdiff documentation built on Jan. 10, 2021, 6:51 a.m.