R Under development (unstable) (2024-01-28 r85838 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. > #### Examples which use the new feature of more than one 'constraint'. > > suppressMessages(library(cobs)) > > ## do *not* show platform info here (as have *.Rout.save), but in 0_pt-ex.R > options(digits = 6) > > if(!dev.interactive(orNone=TRUE)) pdf("multi-constr.pdf") > > source(system.file("util.R", package = "cobs")) > source(system.file(package="Matrix", "test-tools-1.R", mustWork=TRUE)) Loading required package: tools > ##--> tryCatch.W.E(), showProc.time(), assertError(), relErrV(), ... > > > set.seed(908) > x <- seq(-1,2, len = 50) > f.true <- pnorm(2*x) > y <- f.true + rnorm(50)/10 > plot(x,y); lines(x, f.true, col="gray", lwd=2, lty=3) > > ## constraint on derivative at right end: > (con <- rbind(c(2 , max(x), 0))) # f'(x_n) == 0 [,1] [,2] [,3] [1,] 2 2 0 > > ## Using 'trace = 3' --> 'trace = 2' inside drqssbc2() > > ## Regression splines (lambda = 0) > c2 <- cobs(x,y, trace = 3) qbsks2(): Performing general knot selection ... loo.design2(): -> Xeq 50 x 3 (nz = 150 =^= 1%) loo.design2(): -> Xeq 50 x 4 (nz = 150 =^= 0.75%) loo.design2(): -> Xeq 50 x 5 (nz = 150 =^= 0.6%) loo.design2(): -> Xeq 50 x 6 (nz = 150 =^= 0.5%) loo.design2(): -> Xeq 50 x 7 (nz = 150 =^= 0.43%) Deleting unnecessary knots ... loo.design2(): -> Xeq 50 x 4 (nz = 150 =^= 0.75%) loo.design2(): -> Xeq 50 x 4 (nz = 150 =^= 0.75%) loo.design2(): -> Xeq 50 x 3 (nz = 150 =^= 1%) loo.design2(): -> Xeq 50 x 4 (nz = 150 =^= 0.75%) > c2i <- cobs(x,y, constraint = c("increase"), trace = 3) qbsks2(): Performing general knot selection ... loo.design2(): -> Xeq 50 x 3 (nz = 150 =^= 1%) Xieq 2 x 3 (nz = 6 =^= 1%) loo.design2(): -> Xeq 50 x 4 (nz = 150 =^= 0.75%) Xieq 3 x 4 (nz = 9 =^= 0.75%) loo.design2(): -> Xeq 50 x 5 (nz = 150 =^= 0.6%) Xieq 4 x 5 (nz = 12 =^= 0.6%) loo.design2(): -> Xeq 50 x 6 (nz = 150 =^= 0.5%) Xieq 5 x 6 (nz = 15 =^= 0.5%) loo.design2(): -> Xeq 50 x 7 (nz = 150 =^= 0.43%) Xieq 6 x 7 (nz = 18 =^= 0.43%) Deleting unnecessary knots ... loo.design2(): -> Xeq 50 x 5 (nz = 150 =^= 0.6%) Xieq 4 x 5 (nz = 12 =^= 0.6%) loo.design2(): -> Xeq 50 x 5 (nz = 150 =^= 0.6%) Xieq 4 x 5 (nz = 12 =^= 0.6%) loo.design2(): -> Xeq 50 x 5 (nz = 150 =^= 0.6%) Xieq 4 x 5 (nz = 12 =^= 0.6%) loo.design2(): -> Xeq 50 x 4 (nz = 150 =^= 0.75%) Xieq 3 x 4 (nz = 9 =^= 0.75%) loo.design2(): -> Xeq 50 x 4 (nz = 150 =^= 0.75%) Xieq 3 x 4 (nz = 9 =^= 0.75%) loo.design2(): -> Xeq 50 x 5 (nz = 150 =^= 0.6%) Xieq 4 x 5 (nz = 12 =^= 0.6%) > c2c <- cobs(x,y, constraint = c("concave"), trace = 3) qbsks2(): Performing general knot selection ... loo.design2(): -> Xeq 50 x 3 (nz = 150 =^= 1%) Xieq 1 x 3 (nz = 3 =^= 1%) loo.design2(): -> Xeq 50 x 4 (nz = 150 =^= 0.75%) Xieq 2 x 4 (nz = 6 =^= 0.75%) loo.design2(): -> Xeq 50 x 5 (nz = 150 =^= 0.6%) Xieq 3 x 5 (nz = 9 =^= 0.6%) loo.design2(): -> Xeq 50 x 6 (nz = 150 =^= 0.5%) Xieq 4 x 6 (nz = 12 =^= 0.5%) loo.design2(): -> Xeq 50 x 7 (nz = 150 =^= 0.43%) Xieq 5 x 7 (nz = 15 =^= 0.43%) Deleting unnecessary knots ... loo.design2(): -> Xeq 50 x 3 (nz = 150 =^= 1%) Xieq 1 x 3 (nz = 3 =^= 1%) loo.design2(): -> Xeq 50 x 4 (nz = 150 =^= 0.75%) Xieq 2 x 4 (nz = 6 =^= 0.75%) > > c2IC <- cobs(x,y, constraint = c("inc", "concave"), trace = 3) qbsks2(): Performing general knot selection ... loo.design2(): -> Xeq 50 x 3 (nz = 150 =^= 1%) Xieq 3 x 3 (nz = 9 =^= 1%) loo.design2(): -> Xeq 50 x 4 (nz = 150 =^= 0.75%) Xieq 5 x 4 (nz = 15 =^= 0.75%) loo.design2(): -> Xeq 50 x 5 (nz = 150 =^= 0.6%) Xieq 7 x 5 (nz = 21 =^= 0.6%) loo.design2(): -> Xeq 50 x 6 (nz = 150 =^= 0.5%) Xieq 9 x 6 (nz = 27 =^= 0.5%) loo.design2(): -> Xeq 50 x 7 (nz = 150 =^= 0.43%) Xieq 11 x 7 (nz = 33 =^= 0.43%) Deleting unnecessary knots ... loo.design2(): -> Xeq 50 x 3 (nz = 150 =^= 1%) Xieq 3 x 3 (nz = 9 =^= 1%) > ## here, it's the same as just "i": > all.equal(fitted(c2i), fitted(c2IC)) [1] "Mean relative difference: 0.0808156" > > c1 <- cobs(x,y, degree = 1, trace = 3) qbsks2(): Performing general knot selection ... l1.design2(): -> Xeq 50 x 2 (nz = 100 =^= 1%) l1.design2(): -> Xeq 50 x 3 (nz = 100 =^= 0.67%) l1.design2(): -> Xeq 50 x 4 (nz = 100 =^= 0.5%) l1.design2(): -> Xeq 50 x 5 (nz = 100 =^= 0.4%) l1.design2(): -> Xeq 50 x 6 (nz = 100 =^= 0.33%) Deleting unnecessary knots ... l1.design2(): -> Xeq 50 x 4 (nz = 100 =^= 0.5%) l1.design2(): -> Xeq 50 x 4 (nz = 100 =^= 0.5%) l1.design2(): -> Xeq 50 x 4 (nz = 100 =^= 0.5%) l1.design2(): -> Xeq 50 x 5 (nz = 100 =^= 0.4%) > c1i <- cobs(x,y, degree = 1, constraint = c("increase"), trace = 3) qbsks2(): Performing general knot selection ... l1.design2(): -> Xeq 50 x 2 (nz = 100 =^= 1%) Xieq 1 x 2 (nz = 2 =^= 1%) l1.design2(): -> Xeq 50 x 3 (nz = 100 =^= 0.67%) Xieq 2 x 3 (nz = 4 =^= 0.67%) l1.design2(): -> Xeq 50 x 4 (nz = 100 =^= 0.5%) Xieq 3 x 4 (nz = 6 =^= 0.5%) l1.design2(): -> Xeq 50 x 5 (nz = 100 =^= 0.4%) Xieq 4 x 5 (nz = 8 =^= 0.4%) l1.design2(): -> Xeq 50 x 6 (nz = 100 =^= 0.33%) Xieq 5 x 6 (nz = 10 =^= 0.33%) Deleting unnecessary knots ... l1.design2(): -> Xeq 50 x 4 (nz = 100 =^= 0.5%) Xieq 3 x 4 (nz = 6 =^= 0.5%) l1.design2(): -> Xeq 50 x 4 (nz = 100 =^= 0.5%) Xieq 3 x 4 (nz = 6 =^= 0.5%) l1.design2(): -> Xeq 50 x 4 (nz = 100 =^= 0.5%) Xieq 3 x 4 (nz = 6 =^= 0.5%) l1.design2(): -> Xeq 50 x 5 (nz = 100 =^= 0.4%) Xieq 4 x 5 (nz = 8 =^= 0.4%) > c1c <- cobs(x,y, degree = 1, constraint = c("concave"), trace = 3) qbsks2(): Performing general knot selection ... l1.design2(): -> Xeq 50 x 2 (nz = 100 =^= 1%) l1.design2(): -> Xeq 50 x 3 (nz = 100 =^= 0.67%) Xieq 1 x 3 (nz = 3 =^= 1%) l1.design2(): -> Xeq 50 x 4 (nz = 100 =^= 0.5%) Xieq 2 x 4 (nz = 6 =^= 0.75%) l1.design2(): -> Xeq 50 x 5 (nz = 100 =^= 0.4%) Xieq 3 x 5 (nz = 9 =^= 0.6%) l1.design2(): -> Xeq 50 x 6 (nz = 100 =^= 0.33%) Xieq 4 x 6 (nz = 12 =^= 0.5%) Deleting unnecessary knots ... l1.design2(): -> Xeq 50 x 3 (nz = 100 =^= 0.67%) Xieq 1 x 3 (nz = 3 =^= 1%) l1.design2(): -> Xeq 50 x 3 (nz = 100 =^= 0.67%) Xieq 1 x 3 (nz = 3 =^= 1%) l1.design2(): -> Xeq 50 x 2 (nz = 100 =^= 1%) l1.design2(): -> Xeq 50 x 3 (nz = 100 =^= 0.67%) Xieq 1 x 3 (nz = 3 =^= 1%) > > plot(c1) > lines(predict(c1i), col="forest green") > all.equal(fitted(c1), fitted(c1i), tol = 1e-9)# but not 1e-10 [1] TRUE > > ## now gives warning (not error): > c1IC <- cobs(x,y, degree = 1, constraint = c("inc", "concave"), trace = 3) qbsks2(): Performing general knot selection ... l1.design2(): -> Xeq 50 x 2 (nz = 100 =^= 1%) Xieq 1 x 2 (nz = 2 =^= 1%) l1.design2(): -> Xeq 50 x 3 (nz = 100 =^= 0.67%) Xieq 3 x 3 (nz = 7 =^= 0.78%) l1.design2(): -> Xeq 50 x 4 (nz = 100 =^= 0.5%) Xieq 5 x 4 (nz = 12 =^= 0.6%) l1.design2(): -> Xeq 50 x 5 (nz = 100 =^= 0.4%) Xieq 7 x 5 (nz = 17 =^= 0.49%) l1.design2(): -> Xeq 50 x 6 (nz = 100 =^= 0.33%) Xieq 9 x 6 (nz = 22 =^= 0.41%) Deleting unnecessary knots ... l1.design2(): -> Xeq 50 x 3 (nz = 100 =^= 0.67%) Xieq 3 x 3 (nz = 7 =^= 0.78%) l1.design2(): -> Xeq 50 x 3 (nz = 100 =^= 0.67%) Xieq 3 x 3 (nz = 7 =^= 0.78%) l1.design2(): -> Xeq 50 x 2 (nz = 100 =^= 1%) Xieq 1 x 2 (nz = 2 =^= 1%) l1.design2(): -> Xeq 50 x 3 (nz = 100 =^= 0.67%) Xieq 3 x 3 (nz = 7 =^= 0.78%) Warning messages: 1: In l1.design2(x, w, constraint, ptConstr, knots, pw, nrq = n, nl1, : too few knots ==> nk <= 4; could not add constraint 'concave' 2: In l1.design2(x, w, constraint, ptConstr, knots, pw, nrq = n, nl1, : too few knots ==> nk <= 4; could not add constraint 'concave' > > cp2 <- cobs(x,y, pointwise = con, trace = 3) qbsks2(): Performing general knot selection ... loo.design2(): -> Xeq 50 x 3 (nz = 150 =^= 1%) Xieq 2 x 3 (nz = 6 =^= 1%) loo.design2(): -> Xeq 50 x 4 (nz = 150 =^= 0.75%) Xieq 2 x 4 (nz = 6 =^= 0.75%) loo.design2(): -> Xeq 50 x 5 (nz = 150 =^= 0.6%) Xieq 2 x 5 (nz = 6 =^= 0.6%) loo.design2(): -> Xeq 50 x 6 (nz = 150 =^= 0.5%) Xieq 2 x 6 (nz = 6 =^= 0.5%) loo.design2(): -> Xeq 50 x 7 (nz = 150 =^= 0.43%) Xieq 2 x 7 (nz = 6 =^= 0.43%) Deleting unnecessary knots ... loo.design2(): -> Xeq 50 x 4 (nz = 150 =^= 0.75%) Xieq 2 x 4 (nz = 6 =^= 0.75%) loo.design2(): -> Xeq 50 x 4 (nz = 150 =^= 0.75%) Xieq 2 x 4 (nz = 6 =^= 0.75%) loo.design2(): -> Xeq 50 x 5 (nz = 150 =^= 0.6%) Xieq 2 x 5 (nz = 6 =^= 0.6%) > > ## Here, warning ".. 'ifl'.. " on *some* platforms (e.g. Windows 32bit) : > r2i <- tryCatch.W.E( cobs(x,y, constraint = "increase", pointwise = con) ) qbsks2(): Performing general knot selection ... Deleting unnecessary knots ... > cp2i <- r2i$value > if(doExtras()) print(r2i$warning) # not by default as long as have multi-constr.Rout.save > ## when plotting it, we see that it gave a trivial constant!! > cp2c <- cobs(x,y, constraint = "concave", pointwise = con, trace = 3) qbsks2(): Performing general knot selection ... loo.design2(): -> Xeq 50 x 3 (nz = 150 =^= 1%) Xieq 3 x 3 (nz = 9 =^= 1%) loo.design2(): -> Xeq 50 x 4 (nz = 150 =^= 0.75%) Xieq 4 x 4 (nz = 12 =^= 0.75%) loo.design2(): -> Xeq 50 x 5 (nz = 150 =^= 0.6%) Xieq 5 x 5 (nz = 15 =^= 0.6%) loo.design2(): -> Xeq 50 x 6 (nz = 150 =^= 0.5%) Xieq 6 x 6 (nz = 18 =^= 0.5%) loo.design2(): -> Xeq 50 x 7 (nz = 150 =^= 0.43%) Xieq 7 x 7 (nz = 21 =^= 0.43%) Deleting unnecessary knots ... loo.design2(): -> Xeq 50 x 3 (nz = 150 =^= 1%) Xieq 3 x 3 (nz = 9 =^= 1%) > > ## now gives warning (not error): > cp2IC <- cobs(x,y, constraint = c("inc", "concave"), pointwise = con, trace = 3) qbsks2(): Performing general knot selection ... loo.design2(): -> Xeq 50 x 3 (nz = 150 =^= 1%) Xieq 5 x 3 (nz = 15 =^= 1%) loo.design2(): -> Xeq 50 x 4 (nz = 150 =^= 0.75%) Xieq 7 x 4 (nz = 21 =^= 0.75%) loo.design2(): -> Xeq 50 x 5 (nz = 150 =^= 0.6%) Xieq 9 x 5 (nz = 27 =^= 0.6%) loo.design2(): -> Xeq 50 x 6 (nz = 150 =^= 0.5%) Xieq 11 x 6 (nz = 33 =^= 0.5%) loo.design2(): -> Xeq 50 x 7 (nz = 150 =^= 0.43%) Xieq 13 x 7 (nz = 39 =^= 0.43%) Deleting unnecessary knots ... loo.design2(): -> Xeq 50 x 3 (nz = 150 =^= 1%) Xieq 5 x 3 (nz = 15 =^= 1%) Warning message: In cobs(x, y, constraint = c("inc", "concave"), pointwise = con, : drqssbc2(): Not all flags are normal (== 1), ifl : 18 > > cp1 <- cobs(x,y, degree = 1, pointwise = con, trace = 3) qbsks2(): Performing general knot selection ... l1.design2(): -> Xeq 50 x 2 (nz = 100 =^= 1%) Xieq 2 x 2 (nz = 4 =^= 1%) l1.design2(): -> Xeq 50 x 3 (nz = 100 =^= 0.67%) Xieq 2 x 3 (nz = 4 =^= 0.67%) l1.design2(): -> Xeq 50 x 4 (nz = 100 =^= 0.5%) Xieq 2 x 4 (nz = 4 =^= 0.5%) l1.design2(): -> Xeq 50 x 5 (nz = 100 =^= 0.4%) Xieq 2 x 5 (nz = 4 =^= 0.4%) l1.design2(): -> Xeq 50 x 6 (nz = 100 =^= 0.33%) Xieq 2 x 6 (nz = 4 =^= 0.33%) Deleting unnecessary knots ... l1.design2(): -> Xeq 50 x 4 (nz = 100 =^= 0.5%) Xieq 2 x 4 (nz = 4 =^= 0.5%) l1.design2(): -> Xeq 50 x 4 (nz = 100 =^= 0.5%) Xieq 2 x 4 (nz = 4 =^= 0.5%) l1.design2(): -> Xeq 50 x 4 (nz = 100 =^= 0.5%) Xieq 2 x 4 (nz = 4 =^= 0.5%) l1.design2(): -> Xeq 50 x 5 (nz = 100 =^= 0.4%) Xieq 2 x 5 (nz = 4 =^= 0.4%) > cp1i <- cobs(x,y, degree = 1, constraint = "increase", pointwise = con, trace = 3) qbsks2(): Performing general knot selection ... l1.design2(): -> Xeq 50 x 2 (nz = 100 =^= 1%) Xieq 3 x 2 (nz = 6 =^= 1%) l1.design2(): -> Xeq 50 x 3 (nz = 100 =^= 0.67%) Xieq 4 x 3 (nz = 8 =^= 0.67%) l1.design2(): -> Xeq 50 x 4 (nz = 100 =^= 0.5%) Xieq 5 x 4 (nz = 10 =^= 0.5%) l1.design2(): -> Xeq 50 x 5 (nz = 100 =^= 0.4%) Xieq 6 x 5 (nz = 12 =^= 0.4%) l1.design2(): -> Xeq 50 x 6 (nz = 100 =^= 0.33%) Xieq 7 x 6 (nz = 14 =^= 0.33%) Deleting unnecessary knots ... l1.design2(): -> Xeq 50 x 4 (nz = 100 =^= 0.5%) Xieq 5 x 4 (nz = 10 =^= 0.5%) l1.design2(): -> Xeq 50 x 4 (nz = 100 =^= 0.5%) Xieq 5 x 4 (nz = 10 =^= 0.5%) l1.design2(): -> Xeq 50 x 4 (nz = 100 =^= 0.5%) Xieq 5 x 4 (nz = 10 =^= 0.5%) l1.design2(): -> Xeq 50 x 5 (nz = 100 =^= 0.4%) Xieq 6 x 5 (nz = 12 =^= 0.4%) > cp1c <- cobs(x,y, degree = 1, constraint = "concave", pointwise = con, trace = 3) qbsks2(): Performing general knot selection ... l1.design2(): -> Xeq 50 x 2 (nz = 100 =^= 1%) Xieq 2 x 2 (nz = 4 =^= 1%) l1.design2(): -> Xeq 50 x 3 (nz = 100 =^= 0.67%) Xieq 3 x 3 (nz = 7 =^= 0.78%) l1.design2(): -> Xeq 50 x 4 (nz = 100 =^= 0.5%) Xieq 4 x 4 (nz = 10 =^= 0.62%) l1.design2(): -> Xeq 50 x 5 (nz = 100 =^= 0.4%) Xieq 5 x 5 (nz = 13 =^= 0.52%) l1.design2(): -> Xeq 50 x 6 (nz = 100 =^= 0.33%) Xieq 6 x 6 (nz = 16 =^= 0.44%) Deleting unnecessary knots ... l1.design2(): -> Xeq 50 x 3 (nz = 100 =^= 0.67%) Xieq 3 x 3 (nz = 7 =^= 0.78%) l1.design2(): -> Xeq 50 x 3 (nz = 100 =^= 0.67%) Xieq 3 x 3 (nz = 7 =^= 0.78%) l1.design2(): -> Xeq 50 x 2 (nz = 100 =^= 1%) Xieq 2 x 2 (nz = 4 =^= 1%) l1.design2(): -> Xeq 50 x 3 (nz = 100 =^= 0.67%) Xieq 3 x 3 (nz = 7 =^= 0.78%) > > cp1IC <- cobs(x,y, degree = 1, constraint = c("inc", "concave"), pointwise = con, trace = 3) qbsks2(): Performing general knot selection ... l1.design2(): -> Xeq 50 x 2 (nz = 100 =^= 1%) Xieq 3 x 2 (nz = 6 =^= 1%) l1.design2(): -> Xeq 50 x 3 (nz = 100 =^= 0.67%) Xieq 5 x 3 (nz = 11 =^= 0.73%) l1.design2(): -> Xeq 50 x 4 (nz = 100 =^= 0.5%) Xieq 7 x 4 (nz = 16 =^= 0.57%) l1.design2(): -> Xeq 50 x 5 (nz = 100 =^= 0.4%) Xieq 9 x 5 (nz = 21 =^= 0.47%) l1.design2(): -> Xeq 50 x 6 (nz = 100 =^= 0.33%) Xieq 11 x 6 (nz = 26 =^= 0.39%) Deleting unnecessary knots ... l1.design2(): -> Xeq 50 x 3 (nz = 100 =^= 0.67%) Xieq 5 x 3 (nz = 11 =^= 0.73%) l1.design2(): -> Xeq 50 x 3 (nz = 100 =^= 0.67%) Xieq 5 x 3 (nz = 11 =^= 0.73%) l1.design2(): -> Xeq 50 x 2 (nz = 100 =^= 1%) Xieq 3 x 2 (nz = 6 =^= 1%) l1.design2(): -> Xeq 50 x 3 (nz = 100 =^= 0.67%) Xieq 5 x 3 (nz = 11 =^= 0.73%) Warning messages: 1: In l1.design2(x, w, constraint, ptConstr, knots, pw, nrq = n, nl1, : too few knots ==> nk <= 4; could not add constraint 'concave' 2: In l1.design2(x, w, constraint, ptConstr, knots, pw, nrq = n, nl1, : too few knots ==> nk <= 4; could not add constraint 'concave' > > > plot(x,y, main = "cobs(*, degree= 1, constraint = *, pointwise= *)") > matlines(x,cbind(fitted(c1), + fitted(c1i), + fitted(c1c), + fitted(cp1), + fitted(cp1i), + fitted(cp1c)), + col = 1:6, lty=1) > legend("bottomright", inset = .02, col = 1:6, lty=1, + legend = c("none", "increase","concave", + "pt", "pt + incr.", "pt + conc.")) > > if(dev.interactive()) x11() # cheap way to get plot in new window, when testing > > plot(x,y, main = "cobs(*, degree= 2, constraint = *, pointwise= *)") > matlines(x,cbind(fitted(c2), + fitted(c2i), + fitted(c2c), + fitted(cp2), + fitted(cp2i), + fitted(cp2c)), + col = 1:6, lty=1) > legend("bottomright", inset = .02, col = 1:6, lty=1, + legend = c("none", "increase","concave", + "pt", "pt + incr.", "pt + conc.")) > > ##--> "increase + pointwise" gives constant which seems plain wrong <<<< BUG ??? > > proc.time() user system elapsed 2.98 0.31 3.28