## ----setup, include=FALSE----------------------------------------------------- knitr::opts_chunk$set( echo = FALSE, message = FALSE, warning = FALSE, fig.width = 7, fig.height = 4.5, fig.align = "center", dev = "png", dpi = 96 ) library(dplR) # caps() library(Matrix) # the legacy ffcsaps() used sparse matrices ## A colour-blind-safe palette (Okabe and Ito) COLOR_SIM <- "#0072B2" # blue COLOR_ALT <- "#D55E00" # vermillion COLOR_LINE <- "#009E73" # bluish green COLOR_REF <- "#999999" # grey ## ----response-init------------------------------------------------------------ ## Cook, E. R. and Kairiukstis, L. A. (1990) Methods of ## Dendrochronology: Applications in the Environmental Sciences. ## Cook, E. R. and Peters, K. (1981) The smoothing spline: a new ## approach to standardizing forest interior tree-ring width series ## for dendroclimatic studies. ## Smoothing parameter p as a function of the period (nyrs) at which a ## frequency response of f is desired. This is equation (2) below, and ## it is what the Fortran behind caps() computes. pCook <- function(nyrs, f = 0.5) { 6 * f * (cos(2 * pi / nyrs) - 1)^2 / ((1 - f) * (cos(2 * pi / nyrs) + 2)) } ## Frequency response according to Cook and Kairiukstis (citing Cook ## and Peters). This is equation (1) below. respCook <- function(f, p) { pif2 <- 2 * pi * f 1 - 1 / (1 + (p * (cos(pif2) + 2)) / (6 * (cos(pif2) - 1)^2)) } ## ----response-comp------------------------------------------------------------ N <- 1536 K <- 500 NYRS <- c(4, 16, 64) nFreq <- N / 2 + 1 halfseq <- seq_len(nFreq) ratio1 <- array(NA_real_, c(nFreq, K, length(NYRS))) if (!exists(".Random.seed", globalenv(), mode = "numeric")) { foo <- sample(TRUE) } seed <- get(".Random.seed", globalenv()) set.seed(123) for (k in seq_len(K)) { x <- rnorm(N) fftx <- abs(fft(x))[halfseq] for (j in seq_along(NYRS)) { fft1 <- abs(fft(caps(x, nyrs = NYRS[j], f = 0.5)))[halfseq] ratio1[, k, j] <- fft1 / fftx } } assign(".Random.seed", seed, globalenv()) response1 <- matrix(NA_real_, nFreq, length(NYRS)) colnames(response1) <- NYRS for (j in seq_along(NYRS)) { response1[, j] <- rowMeans(ratio1[, , j]) } fftFreq <- seq(from = 0, to = 0.5, length.out = nFreq) ## Simulated response at the nominal cutoff frequency atCutoff <- vapply(seq_along(NYRS), function(j) approx(fftFreq, response1[, j], xout = 1 / NYRS[j])$y, numeric(1)) ## ----response, fig.height=8, fig.cap="**Figure 1.** Theoretical frequency response of the spline filter (equations 1 and 2, green line) versus the response measured with i.i.d. normal series of 1536 samples, mean of 500 repeats, using `caps` (blue circles). The legend on the bottom panel applies to all panels."---- op <- par(mfcol = c(3, 1), mgp = c(2, 0.75, 0), mar = c(4, 4, 2.5, 1)) LWD <- 3 PCH_1 <- 1 ## The simulation has one value per Fourier frequency, which is far too ## dense to plot as points. Show every SUBth one so that the theoretical ## curve underneath stays visible. SUB <- seq(from = 1, to = nFreq, by = 10) for (j in seq_along(NYRS)) { plot(fftFreq, response1[, j], type = "n", ylim = c(0, 1), xlab = "Frequency (1 / year)", ylab = "Amplitude response", main = sprintf("nyrs = %d, f = 0.5", NYRS[j])) lines(fftFreq, respCook(fftFreq, pCook(NYRS[j])), col = COLOR_LINE, lwd = LWD) points(fftFreq[SUB], response1[SUB, j], pch = PCH_1, col = COLOR_SIM, cex = 0.8, lwd = 1.5) abline(h = 0.5, lty = "dashed") abline(v = 1 / NYRS[j], lty = "dashed") text(0.35, 0.5, "50% response", pos = 1, offset = 1) text(1 / NYRS[j], 0.6, sprintf("%d yr period", NYRS[j]), pos = 4, srt = 90, offset = 1) } legend("topright", bg = "white", legend = c("Simulation (caps)", "Theoretical (Cook and Peters)"), col = c(COLOR_SIM, COLOR_LINE), lty = c(NA, "solid"), pch = c(PCH_1, NA), lwd = c(1, LWD)) par(op) ## ----ffcsaps-legacy----------------------------------------------------------- ## Helper used by ffcsaps.legacy() inc <- function(from, to) { if (is.numeric(to) && is.numeric(from) && to >= from) { seq(from = from, to = to) } else { integer(length = 0) } } ## The pure R implementation of ffcsaps() as it stood in dplR 1.6.9, ## before caps() replaced it. Reproduced here only so that the two ## implementations can be compared; use caps() for real work. ffcsaps.legacy <- function(y, x = seq_along(y), nyrs = length(y)/2, f = 0.5) { ffppual <- function(breaks, c1, c2, c3, c4, x, left) { if (left) { ix <- order(x) x2 <- x[ix] } else { x2 <- x } n.breaks <- length(breaks) if (left) { index <- pmax(ffsorted(breaks[-n.breaks], x2), 1) } else { index <- ffsorted2(breaks[-1], x2) } x2 <- x2 - breaks[index] v <- x2 * (x2 * (x2 * c1[index] + c2[index]) + c3[index]) + c4[index] if (left) v[ix] <- v v } ffsorted <- function(meshsites, sites) { index <- order(c(meshsites, sites)) which(index > length(meshsites)) - seq_along(sites) } ffsorted2 <- function(meshsites, sites) { index <- order(c(sites, meshsites)) which(index <= length(sites)) - seq(from = 0, to = length(sites) - 1) } ## Similar in function to spdiags(B, d, n, n) in MATLAB spdiags <- function(B, d, n) { n.d <- length(d) A <- matrix(0, n.d * n, 3) count <- 0 for (k in seq_len(n.d)) { this.diag <- d[k] i <- inc(max(1, 1 - this.diag), min(n, n - this.diag)) n.i <- length(i) if (n.i > 0) { j <- i + this.diag row.idx <- seq(from = count + 1, by = 1, length.out = n.i) A[row.idx, 1] <- i A[row.idx, 2] <- j A[row.idx, 3] <- B[j, k] count <- count + n.i } } A <- A[A[, 3] != 0, , drop = FALSE] A[order(A[, 2], A[, 1]), , drop = FALSE] } y2 <- as.numeric(y) x2 <- as.numeric(x) n <- length(x2) if (n < 3) stop("there must be at least 3 data points") ix <- order(x2) zz1 <- n - 1 xi <- x2[ix] zz2 <- n - 2 diff.xi <- diff(xi) if (any(diff.xi == 0)) stop("the data abscissae must be distinct") if (n != length(y2)) stop("abscissa and ordinate vector must be of the same length") arg2 <- -1:1 odx <- 1 / diff.xi R <- spdiags(cbind(c(diff.xi[-c(1, zz1)], 0), 2 * (diff.xi[-1] + diff.xi[-zz1]), c(0, diff.xi[-c(1, 2)])), arg2, zz2) R2 <- spdiags(cbind(c(odx[-zz1], 0, 0), c(0, -(odx[-1] + odx[-zz1]), 0), c(0, 0, odx[-1])), arg2, n) R2[, 1] <- R2[, 1] - 1 forR <- Matrix(0, zz2, zz2, sparse = TRUE) forR2 <- Matrix(0, zz2, n, sparse = TRUE) forR[R[, 1:2, drop = FALSE]] <- R[, 3] forR2[R2[, 1:2, drop = FALSE]] <- R2[, 3] ## This is equation (4): the ffcsaps parameterization p.inv <- (1 - f) * (cos(2 * pi / nyrs) + 2) / (12 * f * (cos(2 * pi / nyrs) - 1)^2) + 1 yi <- y2[ix] p <- 1 / p.inv mplier <- 6 - 6 / p.inv u <- as.numeric(solve(mplier * tcrossprod(forR2) + forR * p, diff(diff(yi) / diff.xi))) yi <- yi - mplier * diff(c(0, diff(c(0, u, 0)) / diff.xi, 0)) test0 <- xi[-c(1, n)] c3 <- c(0, u / p.inv, 0) x3 <- c(test0, seq(from = xi[1], to = xi[n], length = 101)) cc.1 <- diff(c3) / diff.xi cc.2 <- 3 * c3[-n] cc.3 <- diff(yi) / diff.xi - diff.xi * (2 * c3[-n] + c3[-1]) cc.4 <- yi[-n] to.sort <- c(test0, x3) ix.final <- order(to.sort) tmp <- unique(data.frame( to.sort[ix.final], c(ffppual(xi, cc.1, cc.2, cc.3, cc.4, test0, FALSE), ffppual(xi, cc.1, cc.2, cc.3, cc.4, x3, TRUE))[ix.final])) tmp2 <- tmp tmp2[[1]] <- round(tmp2[[1]], 5) res <- tmp2[[2]][tmp2[[1]] %in% x2] if (length(res) != n) res <- approx(x = tmp[[1]], y = tmp[[2]], xout = x2, ties = "ordered")$y res } ## ----equiv-comp--------------------------------------------------------------- if (!exists(".Random.seed", globalenv(), mode = "numeric")) { foo <- sample(TRUE) } seed <- get(".Random.seed", globalenv()) set.seed(42) data(ca533) cam011 <- as.numeric(na.omit(ca533[, "CAM011"])) cases <- list( "i.i.d. normal, n = 200" = rnorm(200, 100, 20), "noisy sine wave, n = 100" = 5 * sin(seq(from = 0, to = 6 * pi, length.out = 101)[-101]) + rnorm(100) + 20, "AR(1), phi = 0.7, n = 500" = as.numeric(arima.sim(list(ar = 0.7), 500)) + 50, "ca533 series CAM011" = cam011) assign(".Random.seed", seed, globalenv()) equivGrid <- expand.grid(case = names(cases), nyrs = c(10, 32), f = c(0.5, 0.9), stringsAsFactors = FALSE) equivGrid$maxdiff <- vapply(seq_len(nrow(equivGrid)), function(i) { y <- cases[[equivGrid$case[i]]] max(abs(ffcsaps.legacy(y, nyrs = equivGrid$nyrs[i], f = equivGrid$f[i]) - caps(y, nyrs = equivGrid$nyrs[i], f = equivGrid$f[i]))) }, numeric(1)) worstInteger <- max(equivGrid$maxdiff) ## The one place the two differ: caps() coerces nyrs to an integer fracNyrs <- 2 * length(cam011) / 3 diffFrac <- max(abs(ffcsaps.legacy(cam011, nyrs = fracNyrs, f = 0.5) - caps(cam011, nyrs = fracNyrs, f = 0.5))) diffTrunc <- max(abs(ffcsaps.legacy(cam011, nyrs = trunc(fracNyrs), f = 0.5) - caps(cam011, nyrs = fracNyrs, f = 0.5))) camRange <- diff(range(cam011)) ## ----equiv-table-------------------------------------------------------------- ord <- order(match(equivGrid$case, names(cases)), equivGrid$nyrs, equivGrid$f) tab <- equivGrid[ord, c("case", "nyrs", "f", "maxdiff")] tab$maxdiff <- sprintf("%.2e", tab$maxdiff) names(tab) <- c("Series", "nyrs", "f", "max. abs. difference") knitr::kable(tab, row.names = FALSE, align = "lrrr", caption = paste("**Table 1.** Largest absolute difference", "between `ffcsaps.legacy` and `caps` over all", "fitted values, for integer `nyrs`.")) ## ----equiv-fig, fig.height=6, fig.cap="**Figure 2.** Top: series CAM011 of the `ca533` data set (grey) with a 32-year spline fitted by `caps` (blue) and by the legacy `ffcsaps` (vermillion, dashed); the two curves are indistinguishable. Bottom: the difference between the two fitted curves, in ring-width units."---- capsFit <- caps(cam011, nyrs = 32, f = 0.5) ffFit <- ffcsaps.legacy(cam011, nyrs = 32, f = 0.5) op <- par(mfcol = c(2, 1), mgp = c(2, 0.75, 0), mar = c(4, 4, 2.5, 1)) plot(cam011, type = "l", col = COLOR_REF, xlab = "Index", ylab = "Ring width (mm)", main = "CAM011 with a 32-year spline") lines(capsFit, col = COLOR_SIM, lwd = 3) lines(ffFit, col = COLOR_ALT, lwd = 2, lty = "dashed") legend("topright", bty = "n", cex = 0.85, legend = c("data", "caps", "ffcsaps (legacy)"), col = c(COLOR_REF, COLOR_SIM, COLOR_ALT), lty = c("solid", "solid", "dashed"), lwd = c(1, 3, 2)) plot(capsFit - ffFit, type = "l", col = COLOR_SIM, xlab = "Index", ylab = "caps - ffcsaps", main = "Difference between the two fits") abline(h = 0, lty = "dashed", col = COLOR_REF) par(op) ## ----lorenz, fig.width=6, fig.height=6, fig.cap="**Figure 3.** Graphical representation of the Gini coefficient based on areas defined by the Lorenz curve (n = 6). See equations (16), (17), (18) and (19)."---- xg <- sort(c(0.2, 0.4, 0.75, 0.95, 1.2, 2.5)) ng <- length(xg) Xg <- cumsum(xg) px <- c(0, seq_len(ng) / ng) # cumulative portion of population py <- c(0, Xg / Xg[ng]) # cumulative sum of values / total COL_A <- "#FFC0CB" # pink COL_B1 <- "#008080" # teal COL_B2 <- "#00FFFF" # cyan op <- par(mar = c(5, 5, 1, 3), pty = "s") plot(NA, xlim = c(0, 1), ylim = c(0, 1), asp = 1, axes = FALSE, xlab = "Cumulative portion of population\n(ordered from lowest to highest value)", ylab = "Cumulative sum of values\ndivided by total") ## A: between the line of equality and the Lorenz curve polygon(c(px, rev(px)), c(py, rev(px)), col = COL_A, border = NA) ## B2: the staircase of bars under the curve, height = left endpoint for (i in seq_len(ng)) { rect(px[i], 0, px[i + 1], py[i], col = COL_B2, border = NA) } ## B1: the triangles capping each bar for (i in seq_len(ng)) { polygon(c(px[i], px[i + 1], px[i + 1]), c(py[i], py[i], py[i + 1]), col = COL_B1, border = NA) } ## Outlines polygon(c(0, 1, 1), c(0, 0, 1)) # the triangle lines(px, py, lwd = 1.5) # Lorenz curve points(px[-1], py[-1], pch = 21, bg = "white", cex = 0.9) for (i in seq_len(ng)) { # staircase outline lines(c(px[i], px[i + 1], px[i + 1]), c(py[i], py[i], py[i + 1]), col = "grey30", lwd = 0.6) } axis(1, at = px, labels = c("0", paste0(seq_len(ng), "/", ng)), cex.axis = 0.85) axis(4, at = c(0, 1), labels = c("0", "1"), las = 1, cex.axis = 0.85) text(0.42, 0.60, "Line of equality", srt = 45, cex = 0.9) text(0.55, 0.33, "Lorenz curve", srt = 44, cex = 0.9) legend("topleft", bty = "n", cex = 0.9, legend = c("A", expression(B[1]), expression(B[2])), fill = c(COL_A, COL_B1, COL_B2), border = NA) par(op) ## ----gini-rmd, echo=TRUE------------------------------------------------------ ## Gini index is one half of relative mean difference. ## x should not have NA values. gini.rmd <- function(x) { mean(abs(outer(x, x, "-"))) / mean(x) * 0.5 } ## ----gini-check, echo=TRUE---------------------------------------------------- giniMax <- max(abs(vapply(ca533, function(x) { x <- x[!is.na(x)] gini.rmd(x) - gini.coef(x) }, numeric(1)))) giniMax