R version 4.4.0 RC (2024-04-16 r86441 ucrt) -- "Puppy Cup" Copyright (C) 2024 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > require(lokern) Loading required package: lokern > data(xSim) > n <- length(xSim) > stopifnot(n == 75) > tt <- ((1:n) - 1/2)/n # equidistant x > > (gk <- glkerns(tt, xSim)) Call: glkerns.default(x = tt, y = xSim) List of 23 $ x : num [1:75] 0.006666667 0.02 0.033333333 0.046666667 0.06 ... $ y : num [1:75] 1.9666 1.9 1.6449 1.4275 2.4 1.8487 1.3383 1.9514 1.5722 1.1978 ... $ x.out : num [1:300] 0.006666667 0.011032448 0.01539823 0.019764012 0.02 ... $ est : num [1:300] 1.867502 1.860786 1.85407 1.847354 1.846991 ... $ nobs : int 75 $ n.out : int 300 $ deriv : int 0 $ korder : int 2 $ hetero : int 0 $ is.rand : int 1 $ inputb : int 0 $ iter : int 15 $ xl : num 0.067 $ xu : num 0.933 $ s : num [1:76] 0.004117607 0.013333333 0.026666667 0.04 0.053333333 ... $ sig : num 0.1451065 $ bandwidth : num 0.04845241 $ ind.x : int [1:75] 1 5 9 13 17 21 25 29 33 37 ... $ seqXmethod: chr "aim" $ m1 : int 400 $ isOrd : logi TRUE $ ord : NULL $ x.inOut : logi TRUE > summary(gk$est) Min. 1st Qu. Median Mean 3rd Qu. Max. -3.0683 -1.6676 0.7053 0.4094 1.7819 4.2125 > gk$bandwidth [1] 0.04845241 > glkerns(tt,xSim, deriv = 1)$bandwidth [1] 0.08137163 > glkerns(tt,xSim, deriv = 2)$bandwidth [1] 0.1190404 > > if(!grDevices::dev.interactive(orNone=TRUE)) pdf("glk-derivs-etc.pdf") > > demo("glk-derivs", ask = FALSE) demo(glk-derivs) ---- ~~~~~~~~~~ > ## Q: Can we have *same* kernel, *same* bandwidth with different 'deriv' > ## similarly to smooth.spline() ? > ## > ## Answer: not really, mainly because don't have enough choices > ## (nu, k_{ord}), i.e., because currently, nu - k must be even ... > > ## "dput(.) of simple integer vector to character : > myF <- function(d, A="(", O=")", B = length(d) > 1) + paste(if(B) A, paste(d, collapse=", "), if(B) O, sep="") > library(lokern) > p.3glks <- function(x.dat, y.dat, korder, derivs = 0:2, + is.rand=FALSE, useBandwidth, bw.factor = 1.8, + col = 2, lwd = 1.5) + { + ## Purpose: Plot glkerns(*, deriv = {0, 1, 2}) + ## ---------------------------------------------------------------------- + ## Arguments: (x.dat, y.dat): the numeric data vectors + ## korder : the kernel order -- automatically is diminuished by one + ## if needed to keep 'korder - deriv' an even number + ## derivs : integer vectors of derivatives to compute + ## useBandwidth: possibly a user specified bandwidth + ## ---------------------------------------------------------------------- + ## Author: Martin Maechler, Date: 2 Jul 2009, 09:24 + + if(!missing(useBandwidth) && is.numeric(useBandwidth) && useBandwidth > 0) + bw <- useBandwidth + else { + ## Determine the fixed bandwidth : + bw0 <- glkerns(x.dat, y.dat, korder=korder, is.rand=is.rand)$bandwidth + bw <- bw0 * bw.factor # more smoothing for the derivatives + } + + stopifnot(derivs == (d <- as.integer(derivs)), length(derivs) >= 1) + derivs <- d + stopifnot(0 <= derivs, derivs <= 4) + glist <- as.list(derivs) + ## Estimates for g, g' , g'' ... {well, depending on derivs} : + for(i.d in seq_along(derivs)) { + nu <- derivs[i.d] + k0 <- korder + if((korder - nu) %% 2 != 0) { ## 'korder - nu' must be even {theory; Fortran code} + k0 <- korder - 1 + message(gettextf("deriv = %d: modifying korder from %d to %d", + nu, korder, k0)) + } + glist[[i.d]] <- + glkerns(x.dat, y.dat, korder=k0, is.rand=is.rand, bandw = bw, deriv = nu) + } + + names(glist) <- paste("nu=", derivs, sep="") + + ## --------- Plots --------------- + + op <- par(mfrow= c(length(derivs), 1), mgp = c(1.25, 0.6, 0), + mar = c(3,3,2.5,1) + .1, oma = c(0,0, 2, 0)) + on.exit(par(op)) + + for(i.d in seq_along(derivs)) { + nu <- derivs[i.d] + tit <- + switch(nu + 1, + expression(widehat(g)(.)), # 0 + expression(widehat(g * minute)(.)), # 1: g' + expression(widehat(g * second)(.)), # 2: g'' + expression(widehat(g * minute*second)(.)),# 3: g''' + expression(widehat(g ^ {(4)})), + expression(widehat(g ^ {(5)})), + expression(widehat(g ^ {(6)}))) + + with(glist[[i.d]], { + plot(est ~ x.out, type = "l", main = tit, col=col, lwd=lwd) + if(nu == 0) ## data + points(y ~ x, cex = 0.5) + else ## y = 0 line (to see zeros): + abline(h = 0, col = "gray", lty=3) + + mtext(substitute(list(bw == B,k[ord] == K), + list(B = formatC(bandwidth), K = korder)), + adj = 1, line = .25) }) + } + mtext(sprintf("glkerns(*, deriv = %s, bandwidth = , korder = %d)", + myF(derivs, "{","}"), korder), + line = 0.5, outer = TRUE, cex = par("cex.main"), font = par("font.main")) + + invisible(glist) + } > data(xSim) > n <- length(xSim) > tt <- ((1:n) - 1/2)/n # equidistant x > p.3glks(tt, xSim, kord = 4) deriv = 1: modifying korder from 4 to 3 > ## Chose bandwidth yourself; see all available derivatives: > ## Store results > r <- p.3glks(tt, xSim, kord = 6, + derivs = 0:4, useBand = 0.15) deriv = 1: modifying korder from 6 to 5 deriv = 3: modifying korder from 6 to 5 > ## and inspect them > if(interactive()) print(r) > ## ------------ > stopifnot(lengths(r.est <- lapply(r, `[[`, "est")) == 300) > rq <- sapply(r.est, quantile) > stopifnot( + all.equal( + rq, + array(c(-3.05424, -1.741756, 0.6749366, 1.758692, 4.387521, + -39.49606, -7.86697, -5.754756, -2.398544, 34.82773, + -770.4308, -57.48376, 2.249302, 197.8807, 449.9687, + -9002.222, -1537.481, 267.3216, 1870.717, 8486.839, + -198315.3, -24259.81, 18925.91, 87331.48, 216498.5), + dim = c(5,5), dimnames = list(names(quantile(1)), names(r.est))), + tol = 7e-7) # 64b : 9.9e-8 + ) > > p.3glks(tt, xSim, kord = 3) deriv = 0: modifying korder from 3 to 2 deriv = 2: modifying korder from 3 to 2 Warning messages: 1: In .glkerns(x = x, y = y, x.out = x.out, nobs = n, n.out = n.out, : 'korder' reset from 3 to 2, internally 2: In .glkerns(x = x, y = y, x.out = x.out, nobs = n, n.out = n.out, : 'korder' reset from 2 to 4, internally > > p.3glks(tt, xSim, kord = 4, useB = 0.15) deriv = 1: modifying korder from 4 to 3 > > ## Some summary output, but not to too high precision ( platform diffs ) : > > g3k5 <- p.3glks(tt, xSim, kord = 5, useB = 0.12) # k.ord = (4,5,4) => less sensical? deriv = 0: modifying korder from 5 to 4 deriv = 2: modifying korder from 5 to 4 > stopifnot(lengths(gk5.est <- lapply(g3k5, `[[`, "est")) == 300) > gk5.q <- sapply(gk5.est, quantile, prob = (1:9)/10) > print(gk5.q, digits = 5) nu=0 nu=1 nu=2 10% -2.39274 -28.4630 -351.95788 20% -1.89222 -10.4884 -38.50880 30% -1.39110 -7.8650 -34.30356 40% 0.40576 -6.2751 -19.39833 50% 0.71256 -5.2738 -0.57926 60% 1.16704 -4.0419 46.68432 70% 1.61158 -3.1877 97.17764 80% 1.86995 -1.2868 167.94221 90% 3.34403 20.1072 242.05082 > > g3k6 <- p.3glks(tt, xSim, kord = 6, useB = 0.2, derivs = 0:3) # k.ord = (6,5,6, 5) deriv = 1: modifying korder from 6 to 5 deriv = 3: modifying korder from 6 to 5 > stopifnot(lengths(gk6.est <- lapply(g3k6, `[[`, "est")) == 300) > gk6.q <- sapply(gk6.est, quantile, prob = (1:9)/10) > print(gk6.q, digits = 4) nu=0 nu=1 nu=2 nu=3 10% -2.3847 -27.9617 -386.195 -3421.7 20% -1.9068 -11.0193 -54.248 -1595.0 30% -1.4120 -7.6688 -43.218 -1118.2 40% 0.3670 -6.8423 -24.507 -519.6 50% 0.7081 -5.5052 1.395 129.7 60% 1.1887 -4.1167 44.646 775.4 70% 1.6358 -3.2823 114.243 1404.2 80% 1.8669 -0.1737 200.441 1932.6 90% 3.3561 19.7495 291.180 3116.0 > > > ## "FIXME" visually compare with numerical derivatives (e.g. from splines). > > proc.time() user system elapsed 0.40 0.03 0.42