R Under development (unstable) (2024-02-01 r85851 ucrt) -- "Unsuffered Consequences" 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. > suppressMessages(library(cobs)) > > source(system.file("util.R", package = "cobs")) > (doExtra <- doExtras()) [1] FALSE > source(system.file("test-tools-1.R", package="Matrix", mustWork=TRUE)) Loading required package: tools > showProc.time() # timing here (to be faster by default) Time (user system elapsed): 0 0 0 > > data(DublinWind) > attach(DublinWind)##-> speed & day (instead of "wind.x" & "DUB.") > iday <- sort.list(day) > > if(!dev.interactive(orNone=TRUE)) pdf("wind.pdf", width=10) > > stopifnot(identical(day,c(rep(c(rep(1:365,3),1:366),4), + rep(1:365,2)))) > co50.1 <- cobs(day, speed, constraint= "periodic", tau= .5, lambda= 2.2, + degree = 1) > co50.2 <- cobs(day, speed, constraint= "periodic", tau= .5, lambda= 2.2, + degree = 2) > > showProc.time() Time (user system elapsed): 0.32 0.05 0.36 > > plot(day,speed, pch = ".", col = "gray20") > lines(day[iday], fitted(co50.1)[iday], col="orange", lwd = 2) > lines(day[iday], fitted(co50.2)[iday], col="sky blue", lwd = 2) > rug(knots(co50.1), col=3, lwd=2) > > nknots <- 13 > > > if(doExtra) { + ## Compute the quadratic median smoothing B-spline using SIC + ## lambda selection + co.o50 <- + cobs(day, speed, knots.add = TRUE, constraint="periodic", nknots = nknots, + tau = .5, lambda = -1, method = "uniform") + summary(co.o50) # [does print] + + showProc.time() + + op <- par(mfrow = c(3,1), mgp = c(1.5, 0.6,0), mar=.1 + c(3,3:1)) + with(co.o50, plot(pp.sic ~ pp.lambda, type ="o", + col=2, log = "x", main = "co.o50: periodic")) + with(co.o50, plot(pp.sic ~ pp.lambda, type ="o", ylim = robrng(pp.sic), + col=2, log = "x", main = "co.o50: periodic")) + of <- 0.64430538125795 + with(co.o50, plot(pp.sic - of ~ pp.lambda, type ="o", ylim = c(6e-15, 8e-15), + ylab = paste("sic -",formatC(of, dig=14, small.m = "'")), + col=2, log = "x", main = "co.o50: periodic")) + par(op) + } > > showProc.time() Time (user system elapsed): 0.03 0 0.04 > > ## cobs99: Since SIC chooses a lambda that corresponds to the smoothest > ## possible fit, rerun cobs with a larger lstart value > ## (lstart <- log(.Machine$double.xmax)^3) # 3.57 e9 > ## > co.o50. <- + cobs(day,speed, knots.add = TRUE, constraint = "periodic", nknots = 10, + tau = .5, lambda = -1, method = "quantile") Searching for optimal lambda. This may take a while. While you are waiting, here is something you can consider to speed up the process: (a) Use a smaller number of knots; (b) Set lambda==0 to exclude the penalty term; (c) Use a coarser grid by reducing the argument 'lambda.length' from the default value of 25. The algorithm has converged. You might plot() the returned object (which plots 'sic' against 'lambda') to see if you have found the global minimum of the information criterion so that you can determine if you need to adjust any or all of 'lambda.lo', 'lambda.hi' and 'lambda.length' and refit the model. > summary(co.o50.) COBS smoothing spline (degree = 2) from call: cobs(x = day, y = speed, constraint = "periodic", nknots = 10, method = "quantile", tau = 0.5, lambda = -1, knots.add = TRUE) {tau=0.5}-quantile; dimensionality of fit: 7 from {14,13,11,8,7,30} x$knots[1:10]: 0.999635, 41.000000, 82.000000, ... , 366.000365 lambda = 101002.6, selected via SIC, out of 25 ones. coef[1:12]: 1.121550e+01, 1.139573e+01, 1.089025e+01, 9.954427e+00, 8.148158e+00, ... , 5.373106e-04 R^2 = 8.22% ; empirical tau (over all): 3287/6574 = 0.5 (target tau= 0.5) > summary(pc.5 <- predict(co.o50., interval = "both")) z fit cb.lo cb.up Min. : 0.9996 Min. : 7.212 Min. : 6.351 Min. : 7.951 1st Qu.: 92.2498 1st Qu.: 7.790 1st Qu.: 7.000 1st Qu.: 8.600 Median :183.5000 Median : 9.436 Median : 8.555 Median :10.326 Mean :183.5000 Mean : 9.314 Mean : 8.388 Mean :10.241 3rd Qu.:274.7502 3rd Qu.:10.798 3rd Qu.: 9.716 3rd Qu.:11.787 Max. :366.0004 Max. :11.290 Max. :10.347 Max. :13.416 ci.lo ci.up Min. : 6.782 Min. : 7.598 1st Qu.: 7.370 1st Qu.: 8.213 Median : 8.974 Median : 9.901 Mean : 8.830 Mean : 9.798 3rd Qu.:10.197 3rd Qu.:11.311 Max. :10.797 Max. :12.366 > > showProc.time() Time (user system elapsed): 1.59 0.3 1.89 > > if(doExtra) { ## + repeat.delete.add + co.o50.. <- cobs(day,speed, knots.add = TRUE, repeat.delete.add=TRUE, + constraint = "periodic", nknots = 10, + tau = .5, lambda = -1, method = "quantile") + summary(co.o50..) + showProc.time() + } > > co.o9 <- ## Compute the .9 quantile smoothing B-spline + cobs(day,speed,knots.add = TRUE, constraint = "periodic", nknots = 10, + tau = .9,lambda = -1, method = "uniform") Searching for optimal lambda. This may take a while. While you are waiting, here is something you can consider to speed up the process: (a) Use a smaller number of knots; (b) Set lambda==0 to exclude the penalty term; (c) Use a coarser grid by reducing the argument 'lambda.length' from the default value of 25. WARNING: Some lambdas had problems in rq.fit.sfnc(): lambda icyc ifl fidel sum|res|_s k [1,] 1482516 72 18 6199.497 0.0002848763 4 The algorithm has converged. You might plot() the returned object (which plots 'sic' against 'lambda') to see if you have found the global minimum of the information criterion so that you can determine if you need to adjust any or all of 'lambda.lo', 'lambda.hi' and 'lambda.length' and refit the model. Warning message: In cobs(day, speed, knots.add = TRUE, constraint = "periodic", nknots = 10, : drqssbc2(): Not all flags are normal (== 1), ifl : 11111111111111111111118111 > summary(co.o9) COBS smoothing spline (degree = 2) from call: cobs(x = day, y = speed, constraint = "periodic", nknots = 10, method = "uniform", tau = 0.9, lambda = -1, knots.add = TRUE) * Warning in algorithm: some ifl != 1 {tau=0.9}-quantile; dimensionality of fit: 12 from {13,12,10,9,7,4,14} x$knots[1:10]: 0.999635, 41.555556, 82.111111, ... , 366.000365 lambda = 917.6266, selected via SIC, out of 25 ones. coef[1:12]: 19.08631533, 19.27928101, 17.83662511, 18.71788096, 12.75341051, ... , 0.00416216 empirical tau (over all): 5917/6574 = 0.9000608 (target tau= 0.9) > summary(pc.9 <- predict(co.o9,interval = "both")) z fit cb.lo cb.up Min. : 0.9996 Min. :13.00 Min. :11.25 Min. :14.54 1st Qu.: 92.2498 1st Qu.:13.51 1st Qu.:11.89 1st Qu.:15.20 Median :183.5000 Median :16.84 Median :15.02 Median :18.64 Mean :183.5000 Mean :16.27 Mean :14.47 Mean :18.07 3rd Qu.:274.7502 3rd Qu.:18.33 3rd Qu.:16.44 3rd Qu.:19.98 Max. :366.0004 Max. :19.13 Max. :17.28 Max. :23.38 ci.lo ci.up Min. :12.26 Min. :13.69 1st Qu.:12.81 1st Qu.:14.25 Median :16.09 Median :17.61 Mean :15.50 Mean :17.04 3rd Qu.:17.55 3rd Qu.:18.99 Max. :18.33 Max. :20.92 > > showProc.time() Time (user system elapsed): 1.8 0.17 1.97 > > co.o1 <- ## Compute the .1 quantile smoothing B-spline + cobs(day,speed,knots.add = TRUE, constraint = "periodic",nknots = nknots, + tau = .1,lambda = -1, method = "uniform") Searching for optimal lambda. This may take a while. While you are waiting, here is something you can consider to speed up the process: (a) Use a smaller number of knots; (b) Set lambda==0 to exclude the penalty term; (c) Use a coarser grid by reducing the argument 'lambda.length' from the default value of 25. The algorithm has converged. You might plot() the returned object (which plots 'sic' against 'lambda') to see if you have found the global minimum of the information criterion so that you can determine if you need to adjust any or all of 'lambda.lo', 'lambda.hi' and 'lambda.length' and refit the model. > summary(co.o1) COBS smoothing spline (degree = 2) from call: cobs(x = day, y = speed, constraint = "periodic", nknots = nknots, method = "uniform", tau = 0.1, lambda = -1, knots.add = TRUE) {tau=0.1}-quantile; dimensionality of fit: 7 from {17,16,15,12,10,9,7,6,24} x$knots[1:13]: 0.999635, 31.416667, 61.833333, ... , 366.000365 lambda = 101002.6, selected via SIC, out of 25 ones. coef[1:15]: 4.8378305298, 4.9188085699, 4.8723731947, 4.6175513654, 4.1543428036, ... , 0.0002252407 empirical tau (over all): 657/6574 = 0.09993915 (target tau= 0.1) > summary(pc.1 <- predict(co.o1, interval = "both")) z fit cb.lo cb.up Min. : 0.9996 Min. :3.061 Min. :2.099 Min. :3.841 1st Qu.: 92.2498 1st Qu.:3.309 1st Qu.:2.429 1st Qu.:4.180 Median :183.5000 Median :3.998 Median :3.106 Median :4.849 Mean :183.5000 Mean :3.990 Mean :3.041 Mean :4.938 3rd Qu.:274.7502 3rd Qu.:4.673 3rd Qu.:3.622 3rd Qu.:5.578 Max. :366.0004 Max. :4.901 Max. :4.058 Max. :7.157 ci.lo ci.up Min. :2.561 Min. :3.472 1st Qu.:2.828 1st Qu.:3.768 Median :3.545 Median :4.451 Mean :3.494 Mean :4.486 3rd Qu.:4.086 3rd Qu.:5.150 Max. :4.455 Max. :6.050 > > showProc.time() Time (user system elapsed): 1.7 0.19 1.89 > > op <- par(mfrow = c(1,2), mgp = c(1.5, .6,0), mar = .1 + c(3,3,1,1)) > plot(day,speed, pch = 3, cex=0.6, xlab = "DAYS", ylab = "SPEED (knots)") > lines(pc.5, lwd = 2.5, col = 2) > lines(pc.9, lwd = 2., col = "blue") > lines(pc.1, lwd = 2., col = "blue") > plot(day,speed,type = "n",xlab = "DAYS", ylab = "SPEED (knots)") > lines(pc.5, lwd = 1.5) > lines(pc.9, col = 3) > lines(pc.1, col = 3) > abline(v = co.o50.$knots, lty = 3, col = "gray70") > ## rather rug(co.o5$knots, lty = 2) > par(op) > > showProc.time() Time (user system elapsed): 0.03 0 0.03 > > proc.time() user system elapsed 7.15 0.92 8.07