- using R Under development (unstable) (2024-10-22 r87264)
- using platform: x86_64-pc-linux-gnu
- R was compiled by
clang version 19.1.2
flang-new version 19.1.2
- running under: Fedora Linux 36 (Workstation Edition)
- using session charset: UTF-8
- using option ‘--no-stop-on-test-error’
- checking for file ‘cobs/DESCRIPTION’ ... OK
- this is package ‘cobs’ version ‘1.3-8’
- checking package namespace information ... OK
- checking package dependencies ... OK
- checking if this is a source package ... OK
- checking if there is a namespace ... OK
- checking for executable files ... OK
- checking for hidden files and directories ... OK
- checking for portable file names ... OK
- checking for sufficient/correct file permissions ... OK
- checking whether package ‘cobs’ can be installed ... [18s/46s] OK
See the install log for details.
- used C compiler: ‘clang version 19.1.2’
- checking installed package size ... OK
- checking package directory ... OK
- checking DESCRIPTION meta-information ... OK
- checking top-level files ... OK
- checking for left-over files ... OK
- checking index information ... OK
- checking package subdirectories ... OK
- checking code files for non-ASCII characters ... OK
- checking R files for syntax errors ... OK
- checking whether the package can be loaded ... OK
- checking whether the package can be loaded with stated dependencies ... OK
- checking whether the package can be unloaded cleanly ... OK
- checking whether the namespace can be loaded with stated dependencies ... OK
- checking whether the namespace can be unloaded cleanly ... OK
- checking loading without being on the library search path ... OK
- checking whether startup messages can be suppressed ... OK
- checking use of S3 registration ... OK
- checking dependencies in R code ... OK
- checking S3 generic/method consistency ... OK
- checking replacement functions ... OK
- checking foreign function calls ... OK
- checking R code for possible problems ... [27s/63s] OK
- checking Rd files ... OK
- checking Rd metadata ... OK
- checking Rd line widths ... OK
- checking Rd cross-references ... OK
- checking for missing documentation entries ... OK
- checking for code/documentation mismatches ... OK
- checking Rd \usage sections ... OK
- checking Rd contents ... OK
- checking for unstated dependencies in examples ... OK
- checking contents of ‘data’ directory ... OK
- checking data for non-ASCII characters ... OK
- checking LazyData ... OK
- checking data for ASCII and uncompressed saves ... OK
- checking line endings in C/C++/Fortran sources/headers ... OK
- checking line endings in Makefiles ... OK
- checking compilation flags in Makevars ... OK
- checking for GNU extensions in Makefiles ... OK
- checking for portable use of $(BLAS_LIBS) and $(LAPACK_LIBS) ... OK
- checking use of PKG_*FLAGS in Makefiles ... OK
- checking use of SHLIB_OPENMP_*FLAGS in Makefiles ... OK
- checking pragmas in C/C++ headers and code ... OK
- checking compilation flags used ... OK
- checking compiled code ... OK
- checking examples ... [19s/46s] OK
- checking for unstated dependencies in ‘tests’ ... OK
- checking tests ... [65s/167s] ERROR
Running ‘0_pt-ex.R’
Running ‘ex1.R’ [13s/33s]
Running ‘ex2-long.R’ [8s/21s]
Running ‘ex3.R’
Comparing ‘ex3.Rout’ to ‘ex3.Rout.save’ ... OK
Running ‘multi-constr.R’ [7s/19s]
Running ‘roof.R’ [5s/14s]
Running ‘small-ex.R’ [4s/12s]
Comparing ‘small-ex.Rout’ to ‘small-ex.Rout.save’ ... OK
Running ‘spline-ex.R’
Comparing ‘spline-ex.Rout’ to ‘spline-ex.Rout.save’ ... OK
Running ‘temp.R’ [6s/14s]
Comparing ‘temp.Rout’ to ‘temp.Rout.save’ ... OK
Running ‘wind.R’ [10s/27s]
Running the tests in ‘tests/ex2-long.R’ failed.
Complete output:
> ####
> 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()
Time (user system elapsed): 0.003 0 0.003
>
> options(digits = 5)
> if(!dev.interactive(orNone=TRUE)) pdf("ex2.pdf")
>
> set.seed(821)
> x <- round(sort(rnorm(200)), 3) # rounding -> multiple values
> sum(duplicated(x)) # 9
[1] 3
> y <- (fx <- exp(-x)) + rt(200,4)/4
> summaryCobs(cxy <- cobs(x,y, "decrease"))
qbsks2():
Performing general knot selection ...
Deleting unnecessary knots ...
List of 24
$ call : language cobs(x = x, y = y, constraint = "decrease")
$ tau : num 0.5
$ degree : num 2
$ constraint : chr "decrease"
$ ic : chr "AIC"
$ pointwise : NULL
$ select.knots : logi TRUE
$ select.lambda: logi FALSE
$ x : num [1:200] -2.56 -2.14 -1.91 -1.81 -1.78 ...
$ y : num [1:200] 12.7 8.24 6.67 5.88 6.42 ...
$ resid : num [1:200] 0.72 -0.149 0 -0.195 0.545 ...
$ fitted : num [1:200] 11.98 8.39 6.67 6.07 5.87 ...
$ coef : num [1:5] 11.9769 3.5917 1.0544 0.0295 0.0295
$ knots : num [1:4] -2.557 -0.813 0.418 2.573
$ k0 : num 5
$ k : num 5
$ x.ps :Formal class 'matrix.csr' [package "SparseM"] with 4 slots
$ SSy : num 488
$ lambda : num 0
$ icyc : int 11
$ ifl : int 1
$ pp.lambda : NULL
$ pp.sic : NULL
$ i.mask : NULL
cb.lo ci.lo fit ci.up cb.up
1 11.4448128 11.6875576 11.976923 12.26629 12.50903
2 10.9843366 11.2126114 11.484728 11.75684 11.98512
3 10.5344633 10.7489871 11.004712 11.26044 11.47496
4 10.0951784 10.2966768 10.536874 10.77707 10.97857
5 9.6664684 9.8556730 10.081215 10.30676 10.49596
6 9.2483213 9.4259693 9.637736 9.84950 10.02715
7 8.8407282 9.0075609 9.206435 9.40531 9.57214
8 8.4436848 8.6004453 8.787313 8.97418 9.13094
9 8.0571928 8.2046236 8.380369 8.55612 8.70355
10 7.6812627 7.8201015 7.985605 8.15111 8.28995
11 7.3159159 7.4468904 7.603020 7.75915 7.89012
12 6.9611870 7.0850095 7.232613 7.38022 7.50404
13 6.6171269 6.7344861 6.874385 7.01428 7.13164
14 6.2838041 6.3953578 6.528336 6.66131 6.77287
15 5.9613061 6.0676719 6.194466 6.32126 6.42763
16 5.6497392 5.7514863 5.872775 5.99406 6.09581
17 5.3492272 5.4468683 5.563262 5.67966 5.77730
18 5.0599086 5.1538933 5.265928 5.37796 5.47195
19 4.7819325 4.8726424 4.980774 5.08891 5.17961
20 4.5154542 4.6031999 4.707798 4.81240 4.90014
21 4.2606295 4.3456507 4.447001 4.54835 4.63337
22 4.0176099 4.1000771 4.198383 4.29669 4.37916
23 3.7865383 3.8665567 3.961943 4.05733 4.13735
24 3.5675443 3.6451602 3.737683 3.83021 3.90782
25 3.3607413 3.4359491 3.525601 3.61525 3.69046
26 3.1662231 3.2389744 3.325698 3.41242 3.48517
27 2.9840608 3.0542750 3.137974 3.22167 3.29189
28 2.8142997 2.8818753 2.962429 3.04298 3.11056
29 2.6569546 2.7217833 2.799063 2.87634 2.94117
30 2.5120031 2.5739870 2.647875 2.72176 2.78375
31 2.3793776 2.4384496 2.508867 2.57928 2.63836
32 2.2589520 2.3151025 2.382037 2.44897 2.50512
33 2.1505256 2.2038366 2.267386 2.33094 2.38425
34 2.0538038 2.1044916 2.164914 2.22534 2.27602
35 1.9677723 2.0162522 2.074043 2.13183 2.18031
36 1.8846710 1.9316617 1.987677 2.04369 2.09068
37 1.8024456 1.8486425 1.903712 1.95878 2.00498
38 1.7213655 1.7673410 1.822146 1.87695 1.92293
39 1.6417290 1.6879196 1.742982 1.79804 1.84423
40 1.5638322 1.6105393 1.666217 1.72189 1.76860
41 1.4879462 1.5353474 1.591852 1.64836 1.69576
42 1.4143040 1.4624707 1.519888 1.57731 1.62547
43 1.3430975 1.3920136 1.450324 1.50864 1.55755
44 1.2744792 1.3240589 1.383161 1.44226 1.49184
45 1.2085658 1.2586702 1.318397 1.37812 1.42823
46 1.1454438 1.1958944 1.256034 1.31617 1.36662
47 1.0851730 1.1357641 1.196072 1.25638 1.30697
48 1.0277900 1.0782992 1.138509 1.19872 1.24923
49 0.9733099 1.0235079 1.083347 1.14319 1.19338
50 0.9217268 0.9713870 1.030585 1.08978 1.13944
51 0.8730129 0.9219214 0.980223 1.03852 1.08743
52 0.8271160 0.8750827 0.932262 0.98944 1.03741
53 0.7839554 0.8308269 0.886700 0.94257 0.98945
54 0.7434158 0.7890916 0.843540 0.89799 0.94366
55 0.7053406 0.7497913 0.802779 0.85577 0.90022
56 0.6695233 0.7128138 0.764419 0.81602 0.85931
57 0.6357022 0.6780170 0.728459 0.77890 0.82121
58 0.6035616 0.6452289 0.694899 0.74457 0.78624
59 0.5724566 0.6139693 0.663455 0.71294 0.75445
60 0.5410437 0.5829503 0.632905 0.68286 0.72477
61 0.5094333 0.5521679 0.603110 0.65405 0.69679
62 0.4778879 0.5217649 0.574069 0.62637 0.67025
63 0.4466418 0.4918689 0.545782 0.59970 0.64492
64 0.4158910 0.4625864 0.518250 0.57391 0.62061
65 0.3857918 0.4340022 0.491472 0.54894 0.59715
66 0.3564634 0.4061813 0.465448 0.52471 0.57443
67 0.3279928 0.3791711 0.440179 0.50119 0.55236
68 0.3004403 0.3530042 0.415663 0.47832 0.53089
69 0.2738429 0.3277009 0.391903 0.45610 0.50996
70 0.2482184 0.3032707 0.368896 0.43452 0.48957
71 0.2235676 0.2797141 0.346644 0.41357 0.46972
72 0.1998762 0.2570233 0.325146 0.39327 0.45042
73 0.1771158 0.2351830 0.304402 0.37362 0.43169
74 0.1552452 0.2141706 0.284413 0.35466 0.41358
75 0.1342101 0.1939567 0.265178 0.33640 0.39615
76 0.1139444 0.1745054 0.246697 0.31889 0.37945
77 0.0943704 0.1557743 0.228971 0.30217 0.36357
78 0.0753996 0.1377153 0.211999 0.28628 0.34860
79 0.0569347 0.1202755 0.195781 0.27129 0.33463
80 0.0388708 0.1033980 0.180318 0.25724 0.32177
81 0.0210989 0.0870233 0.165609 0.24419 0.31012
82 0.0035089 0.0710917 0.151654 0.23222 0.29980
83 -0.0140062 0.0555449 0.138454 0.22136 0.29091
84 -0.0315470 0.0403283 0.126008 0.21169 0.28356
85 -0.0492034 0.0253928 0.114316 0.20324 0.27783
86 -0.0670524 0.0106968 0.103378 0.19606 0.27381
87 -0.0851561 -0.0037936 0.093195 0.19018 0.27155
88 -0.1035613 -0.0181039 0.083766 0.18564 0.27109
89 -0.1223000 -0.0322515 0.075091 0.18243 0.27248
90 -0.1413914 -0.0462467 0.067171 0.18059 0.27573
91 -0.1608432 -0.0600938 0.060005 0.18010 0.28085
92 -0.1806546 -0.0737923 0.053594 0.18098 0.28784
93 -0.2008180 -0.0873382 0.047936 0.18321 0.29669
94 -0.2213213 -0.1007247 0.043033 0.18679 0.30739
95 -0.2421494 -0.1139438 0.038884 0.19171 0.31992
96 -0.2632855 -0.1269863 0.035490 0.19797 0.33427
97 -0.2847123 -0.1398427 0.032850 0.20554 0.35041
98 -0.3064126 -0.1525038 0.030964 0.21443 0.36834
99 -0.3283696 -0.1649603 0.029833 0.22463 0.38804
100 -0.3505674 -0.1772037 0.029456 0.23611 0.40948
knots :
[1] -2.557 -0.813 0.418 2.573
coef :
[1] 11.976924 3.591747 1.054378 0.029456 0.029456
> 1 - sum(cxy $ resid ^ 2) / sum((y - mean(y))^2) # R^2 = 97.6%
[1] 0.95969
> showProc.time()
Time (user system elapsed): 0.783 0.089 2.113
>
> if(doExtra) {
+ ## Interpolation
+ cxyI <- cobs(x,y, "decrease", knots = unique(x))
+ ## takes quite long : 63 sec. (Pent. III, 700 MHz) --- this is because
+ ## each knot is added sequentially... {{improve!}}
+
+ summaryCobs(cxyI)# only 7 knots remaining!
+ showProc.time()
+ }
>
> summaryCobs(cxy1 <- cobs(x,y, "decrease", lambda = 0.1))
List of 24
$ call : language cobs(x = x, y = y, constraint = "decrease", lambda = 0.1)
$ tau : num 0.5
$ degree : num 2
$ constraint : chr "decrease"
$ ic : NULL
$ pointwise : NULL
$ select.knots : logi TRUE
$ select.lambda: logi FALSE
$ x : num [1:200] -2.56 -2.14 -1.91 -1.81 -1.78 ...
$ y : num [1:200] 12.7 8.24 6.67 5.88 6.42 ...
$ resid : num [1:200] 0 -0.315 0 -0.161 0.586 ...
$ fitted : num [1:200] 12.7 8.56 6.67 6.04 5.83 ...
$ coef : num [1:22] 12.7 5.78 3.16 2.43 2.11 ...
$ knots : num [1:20] -2.557 -1.34 -1.03 -0.901 -0.772 ...
$ k0 : int 15
$ k : int 15
$ x.ps :Formal class 'matrix.csr' [package "SparseM"] with 4 slots
$ SSy : num 488
$ lambda : num 0.1
$ icyc : int 23
$ ifl : int 1
$ pp.lambda : NULL
$ pp.sic : NULL
$ i.mask : NULL
cb.lo ci.lo fit ci.up cb.up
1 12.0912847 12.4849933 12.6970034 12.90901 13.30272
2 11.5452819 11.9166521 12.1166331 12.31661 12.68798
3 11.0146966 11.3650966 11.5537853 11.74247 12.09287
4 10.4995535 10.8303355 11.0084599 11.18658 11.51737
5 9.9998870 10.3123808 10.4806571 10.64893 10.96143
6 9.5157430 9.8112485 9.9703768 10.12951 10.42501
7 9.0471805 9.3269594 9.4776191 9.62828 9.90806
8 8.5942728 8.8595392 9.0023838 9.14523 9.41049
9 8.1571088 8.4090188 8.5446710 8.68032 8.93223
10 7.7357927 7.9754347 8.1044808 8.23353 8.47317
11 7.3304438 7.5588289 7.6818131 7.80480 8.03318
12 6.9411951 7.1592477 7.2766679 7.39409 7.61214
13 6.5681906 6.7767415 6.8890452 7.00135 7.20990
14 6.2115819 6.4113636 6.5189450 6.62653 6.82631
15 5.8715240 6.0631680 6.1663674 6.26957 6.46121
16 5.5481704 5.7322086 5.8313123 5.93042 6.11445
17 5.2416676 5.4185366 5.5137796 5.60902 5.78589
18 4.9521494 5.1221988 5.2137695 5.30534 5.47539
19 4.6797308 4.8432355 4.9312819 5.01933 5.18283
20 4.4245017 4.5816781 4.6663169 4.75096 4.90813
21 4.1865199 4.3375470 4.4188743 4.50020 4.65123
22 3.9658032 4.1108482 4.1889542 4.26706 4.41211
23 3.7623206 3.9015710 3.9765567 4.05154 4.19079
24 3.5759813 3.7096836 3.7816817 3.85368 3.98738
25 3.4043771 3.5329043 3.6021155 3.67133 3.79985
26 3.2347309 3.3585931 3.4252922 3.49199 3.61585
27 3.0652721 3.1848437 3.2492325 3.31362 3.43319
28 2.8962030 3.0117271 3.0739363 3.13615 3.25167
29 2.7276530 2.8392885 2.8994037 2.95952 3.07115
30 2.5596612 2.6675415 2.7256346 2.78373 2.89161
31 2.3944947 2.4988186 2.5549966 2.61117 2.71550
32 2.2444821 2.3455939 2.4000421 2.45449 2.55560
33 2.1114672 2.2097080 2.2626102 2.31551 2.41375
34 1.9954176 2.0911496 2.1427009 2.19425 2.28998
35 1.8963846 1.9899366 2.0403140 2.09069 2.18424
36 1.8125024 1.9041996 1.9535781 2.00296 2.09465
37 1.7347658 1.8248332 1.8733340 1.92183 2.01190
38 1.6620975 1.7506630 1.7983550 1.84605 1.93461
39 1.5945123 1.6816941 1.7286411 1.77559 1.86277
40 1.5278221 1.6138190 1.6601279 1.70644 1.79243
41 1.4573347 1.5423451 1.5881227 1.63390 1.71891
42 1.3839943 1.4682138 1.5135655 1.55892 1.64314
43 1.3227219 1.4063482 1.4513806 1.49641 1.58004
44 1.2787473 1.3619265 1.4067181 1.45151 1.53469
45 1.2488624 1.3317463 1.3763789 1.42101 1.50390
46 1.2168724 1.2994789 1.3439621 1.38845 1.47105
47 1.1806389 1.2628708 1.3071522 1.35143 1.43367
48 1.1401892 1.2219316 1.2659495 1.30997 1.39171
49 1.0941843 1.1754044 1.2191410 1.26288 1.34410
50 1.0326549 1.1134412 1.1569442 1.20045 1.28123
51 0.9535058 1.0339215 1.0772249 1.12053 1.20094
52 0.8632281 0.9433870 0.9865521 1.02972 1.10988
53 0.7875624 0.8676441 0.9107678 0.95389 1.03397
54 0.7267897 0.8069673 0.8501425 0.89332 0.97350
55 0.6673925 0.7477244 0.7909827 0.83424 0.91457
56 0.6072642 0.6877460 0.7310850 0.77442 0.85491
57 0.5471548 0.6278279 0.6712700 0.71471 0.79539
58 0.4995140 0.5804770 0.6240752 0.66767 0.74864
59 0.4686435 0.5499607 0.5937495 0.63754 0.71886
60 0.4531016 0.5348803 0.5789177 0.62296 0.70473
61 0.4381911 0.5206110 0.5649937 0.60938 0.69180
62 0.4199957 0.5032331 0.5480561 0.59288 0.67612
63 0.4036491 0.4879280 0.5333117 0.57870 0.66297
64 0.3952493 0.4807890 0.5268517 0.57291 0.65845
65 0.3926229 0.4796600 0.5265291 0.57340 0.66044
66 0.3900185 0.4787485 0.5265291 0.57431 0.66304
67 0.3870480 0.4776752 0.5264774 0.57528 0.66591
68 0.3738545 0.4665585 0.5164792 0.56640 0.65910
69 0.3432056 0.4380737 0.4891596 0.54025 0.63511
70 0.2950830 0.3922142 0.4445189 0.49682 0.59395
71 0.2295290 0.3291123 0.3827373 0.43636 0.53595
72 0.1670195 0.2693294 0.3244228 0.37952 0.48183
73 0.1216565 0.2269375 0.2836308 0.34032 0.44561
74 0.0934100 0.2019260 0.2603613 0.31880 0.42731
75 0.0787462 0.1907702 0.2510947 0.31142 0.42344
76 0.0658428 0.1813823 0.2435998 0.30582 0.42136
77 0.0538230 0.1727768 0.2368329 0.30089 0.41984
78 0.0427388 0.1649719 0.2307938 0.29662 0.41885
79 0.0325663 0.1579592 0.2254827 0.29301 0.41840
80 0.0232151 0.1517072 0.2208995 0.29009 0.41858
81 0.0145359 0.1461634 0.2170442 0.28792 0.41955
82 0.0063272 0.1412575 0.2139168 0.28658 0.42151
83 -0.0016568 0.1369034 0.2115173 0.28613 0.42469
84 -0.0096967 0.1330028 0.2098457 0.28669 0.42939
85 -0.0180957 0.1294496 0.2089021 0.28835 0.43590
86 -0.0272134 0.1260791 0.2086264 0.29117 0.44447
87 -0.0387972 0.1210358 0.2071052 0.29317 0.45301
88 -0.0534279 0.1135207 0.2034217 0.29332 0.46027
89 -0.0709531 0.1035871 0.1975762 0.29157 0.46611
90 -0.0912981 0.0912612 0.1895684 0.28788 0.47043
91 -0.1144525 0.0765465 0.1793985 0.28225 0.47325
92 -0.1404576 0.0594287 0.1670665 0.27470 0.47459
93 -0.1693951 0.0398791 0.1525723 0.26527 0.47454
94 -0.2013769 0.0178586 0.1359159 0.25397 0.47321
95 -0.2365365 -0.0066795 0.1170974 0.24087 0.47073
96 -0.2750210 -0.0337868 0.0961167 0.22602 0.46725
97 -0.3169840 -0.0635170 0.0729738 0.20946 0.46293
98 -0.3625797 -0.0959240 0.0476688 0.19126 0.45792
99 -0.4119579 -0.1310604 0.0202016 0.17146 0.45236
100 -0.4652595 -0.1689754 -0.0094278 0.15012 0.44640
knots :
[1] -2.557 -1.340 -1.030 -0.901 -0.772 -0.586 -0.448 -0.305 -0.092 0.054
[11] 0.163 0.329 0.481 0.606 0.722 0.859 1.065 1.244 1.837 2.573
coef :
[1] 12.6970048 5.7788265 3.1620633 2.4291174 2.1069607 1.8462166
[7] 1.6371062 1.4304905 1.3348346 1.1758220 0.9413974 0.7863913
[13] 0.5998958 0.5697029 0.5265291 0.5265291 0.5265291 0.2707227
[19] 0.2086712 0.2086712 -0.0094278 6.5257497
> 1 - sum(cxy1 $ resid ^ 2) / sum((y - mean(y))^2) # R^2 = 98.2%
[1] 0.96169
>
> summaryCobs(cxy2 <- cobs(x,y, "decrease", lambda = 1e-2))
List of 24
$ call : language cobs(x = x, y = y, constraint = "decrease", lambda = 0.01)
$ tau : num 0.5
$ degree : num 2
$ constraint : chr "decrease"
$ ic : NULL
$ pointwise : NULL
$ select.knots : logi TRUE
$ select.lambda: logi FALSE
$ x : num [1:200] -2.56 -2.14 -1.91 -1.81 -1.78 ...
$ y : num [1:200] 12.7 8.24 6.67 5.88 6.42 ...
$ resid : num [1:200] 0 -0.146 0.1468 -0.0463 0.6868 ...
$ fitted : num [1:200] 12.7 8.39 6.52 5.92 5.73 ...
$ coef : num [1:22] 12.7 5.34 3.59 2.19 2.13 ...
$ knots : num [1:20] -2.557 -1.34 -1.03 -0.901 -0.772 ...
$ k0 : int 21
$ k : int 21
$ x.ps :Formal class 'matrix.csr' [package "SparseM"] with 4 slots
$ SSy : num 488
$ lambda : num 0.01
$ icyc : int 35
$ ifl : int 1
$ pp.lambda : NULL
$ pp.sic : NULL
$ i.mask : NULL
cb.lo ci.lo fit ci.up cb.up
1 12.0477594 12.4997491 12.6970071 12.89427 13.34625
2 11.4687308 11.8950752 12.0811411 12.26721 12.69355
3 10.9090823 11.3113523 11.4869116 11.66247 12.06474
4 10.3688404 10.7485883 10.9143185 11.08005 11.45980
5 9.8480420 10.2067945 10.3633618 10.51993 10.87868
6 9.3467363 9.6859859 9.8340417 9.98210 10.32135
7 8.8649866 9.1861815 9.3263579 9.46653 9.78773
8 8.4028715 8.7074055 8.8403106 8.97322 9.27775
9 7.9604861 8.2496865 8.3758998 8.50211 8.79131
10 7.5379421 7.8130586 7.9331254 8.05319 8.32831
11 7.1353676 7.3975607 7.5119874 7.62641 7.88861
12 6.7529050 7.0032361 7.1124859 7.22174 7.47207
13 6.3907086 6.6301316 6.7346209 6.83911 7.07853
14 6.0489410 6.2782966 6.3783923 6.47849 6.70784
15 5.7277684 5.9477816 6.0438001 6.13982 6.35983
16 5.4273551 5.6386366 5.7308444 5.82305 6.03433
17 5.1478583 5.3509094 5.4395252 5.52814 5.73119
18 4.8894214 5.0846433 5.1698424 5.25504 5.45026
19 4.6521676 4.8398760 4.9217960 5.00372 5.19142
20 4.4361933 4.6166367 4.6953861 4.77414 4.95458
21 4.2415605 4.4149443 4.4906127 4.56628 4.73966
22 4.0682883 4.2348044 4.3074756 4.38015 4.54666
23 3.9163432 4.0762071 4.1459751 4.21574 4.37561
24 3.7856282 3.9391227 4.0061110 4.07310 4.22659
25 3.6683774 3.8159306 3.8803259 3.94472 4.09227
26 3.5214653 3.6636629 3.7257209 3.78778 3.92998
27 3.3383583 3.4756303 3.5355387 3.59545 3.73272
28 3.1192735 3.2518988 3.3097793 3.36766 3.50028
29 2.8643493 2.9925103 3.0484425 3.10437 3.23254
30 2.5736278 2.6974778 2.7515286 2.80558 2.92943
31 2.2696062 2.3893733 2.4416422 2.49391 2.61368
32 2.0718959 2.1879754 2.2386350 2.28929 2.40537
33 1.9979346 2.1107181 2.1599392 2.20916 2.32194
34 1.9710324 2.0809358 2.1288999 2.17686 2.28677
35 1.9261503 2.0335510 2.0804229 2.12729 2.23470
36 1.8645775 1.9698487 2.0157914 2.06173 2.16701
37 1.7927585 1.8961587 1.9412848 1.98641 2.08981
38 1.7116948 1.8133707 1.8577443 1.90212 2.00379
39 1.6214021 1.7214896 1.7651699 1.80885 1.90894
40 1.5242004 1.6229275 1.6660141 1.70910 1.80783
41 1.4229217 1.5205162 1.5631086 1.60570 1.70330
42 1.3194940 1.4161806 1.4583766 1.50057 1.59726
43 1.2442053 1.3402109 1.3821098 1.42401 1.52001
44 1.2075941 1.3030864 1.3447613 1.38644 1.48193
45 1.2023778 1.2975311 1.3390581 1.38059 1.47574
46 1.1914924 1.2863272 1.3277152 1.36910 1.46394
47 1.1698641 1.2642688 1.3054691 1.34667 1.44107
48 1.1375221 1.2313649 1.2723199 1.31327 1.40712
49 1.0934278 1.1866710 1.2273643 1.26806 1.36130
50 1.0300956 1.1228408 1.1633168 1.20379 1.29654
51 0.9459780 1.0382977 1.0785880 1.11888 1.21120
52 0.8492712 0.9412961 0.9814577 1.02162 1.11364
53 0.7724392 0.8643755 0.9044985 0.94462 1.03656
54 0.7154255 0.8074718 0.8476428 0.88781 0.97986
55 0.6587891 0.7510125 0.7912608 0.83151 0.92373
56 0.5994755 0.6918710 0.7321944 0.77252 0.86491
57 0.5383570 0.6309722 0.6713915 0.71181 0.80443
58 0.4898228 0.5827709 0.6233354 0.66390 0.75685
59 0.4588380 0.5521926 0.5929345 0.63368 0.72703
60 0.4438719 0.5377564 0.5787296 0.61970 0.71359
61 0.4293281 0.5239487 0.5652432 0.60654 0.70116
62 0.4110511 0.5066103 0.5483143 0.59002 0.68558
63 0.3944126 0.4911673 0.5333932 0.57562 0.67237
64 0.3857958 0.4839980 0.5268556 0.56971 0.66792
65 0.3830000 0.4829213 0.5265291 0.57014 0.67006
66 0.3802084 0.4820731 0.5265291 0.57099 0.67285
67 0.3770181 0.4810608 0.5264673 0.57187 0.67592
68 0.3616408 0.4680678 0.5145149 0.56096 0.66739
69 0.3254129 0.4343244 0.4818557 0.52939 0.63830
70 0.2683149 0.3798245 0.4284897 0.47715 0.58866
71 0.1904294 0.3047541 0.3546478 0.40454 0.51887
72 0.1179556 0.2354105 0.2866704 0.33793 0.45539
73 0.0689088 0.1897746 0.2425231 0.29527 0.41614
74 0.0432569 0.1678366 0.2222059 0.27658 0.40115
75 0.0359906 0.1645977 0.2207246 0.27685 0.40546
76 0.0301934 0.1628364 0.2207246 0.27861 0.41126
77 0.0245630 0.1611257 0.2207246 0.28032 0.41689
78 0.0191553 0.1594827 0.2207246 0.28197 0.42229
79 0.0139446 0.1578996 0.2207246 0.28355 0.42750
80 0.0088340 0.1563468 0.2207246 0.28510 0.43262
81 0.0036634 0.1547759 0.2207246 0.28667 0.43779
82 -0.0017830 0.1531211 0.2207246 0.28833 0.44323
83 -0.0077688 0.1513025 0.2207246 0.29015 0.44922
84 -0.0145948 0.1492286 0.2207246 0.29222 0.45604
85 -0.0225859 0.1468007 0.2207246 0.29465 0.46404
86 -0.0321107 0.1438739 0.2206774 0.29748 0.47347
87 -0.0445016 0.1389916 0.2190720 0.29915 0.48265
88 -0.0601227 0.1315395 0.2151851 0.29883 0.49049
89 -0.0788103 0.1215673 0.2090164 0.29647 0.49684
90 -0.1004844 0.1090993 0.2005661 0.29203 0.50162
91 -0.1251339 0.0941388 0.1898342 0.28553 0.50480
92 -0.1528032 0.0766725 0.1768206 0.27697 0.50644
93 -0.1835797 0.0566736 0.1615253 0.26638 0.50663
94 -0.2175834 0.0341058 0.1439484 0.25379 0.50548
95 -0.2549574 0.0089256 0.1240898 0.23925 0.50314
96 -0.2958592 -0.0189149 0.1019496 0.22281 0.49976
97 -0.3404537 -0.0494657 0.0775277 0.20452 0.49551
98 -0.3889062 -0.0827771 0.0508241 0.18443 0.49055
99 -0.4413769 -0.1188979 0.0218389 0.16258 0.48505
100 -0.4980173 -0.1578738 -0.0094279 0.13902 0.47916
knots :
[1] -2.557 -1.340 -1.030 -0.901 -0.772 -0.586 -0.448 -0.305 -0.092 0.054
[11] 0.163 0.329 0.481 0.606 0.722 0.859 1.065 1.244 1.837 2.573
coef :
[1] 12.697009 5.337850 3.591398 2.187733 2.133993 1.936435 1.631856
[8] 1.340650 1.340650 1.185401 0.931750 0.789326 0.598245 0.570221
[15] 0.526529 0.526529 0.526529 0.220725 0.220725 0.220725 -0.009428
[22] 46.342964
> 1 - sum(cxy2 $ resid ^ 2) / sum((y - mean(y))^2) # R^2 = 98.2% (tiny bit better)
[1] 0.96257
>
> summaryCobs(cxy3 <- cobs(x,y, "decrease", lambda = 1e-6, nknots = 60))
List of 24
$ call : language cobs(x = x, y = y, constraint = "decrease", nknots = 60, lambda = 1e-06)
$ tau : num 0.5
$ degree : num 2
$ constraint : chr "decrease"
$ ic : NULL
$ pointwise : NULL
$ select.knots : logi TRUE
$ select.lambda: logi FALSE
$ x : num [1:200] -2.56 -2.14 -1.91 -1.81 -1.78 ...
$ y : num [1:200] 12.7 8.24 6.67 5.88 6.42 ...
$ resid : num [1:200] 0 0 0 -0.382 0.309 ...
$ fitted : num [1:200] 12.7 8.24 6.67 6.26 6.11 ...
$ coef : num [1:62] 12.7 7.69 6.09 4.35 3.73 3.73 2.74 2.57 2.57 2.25 ...
$ knots : num [1:60] -2.56 -1.81 -1.73 -1.38 -1.23 ...
$ k0 : int 61
$ k : int 61
$ x.ps :Formal class 'matrix.csr' [package "SparseM"] with 4 slots
$ SSy : num 488
$ lambda : num 1e-06
$ icyc : int 46
$ ifl : int 1
$ pp.lambda : NULL
$ pp.sic : NULL
$ i.mask : NULL
cb.lo ci.lo fit ci.up cb.up
1 12.0247124 12.56890432 12.6970139 12.825123 13.36932
2 11.3797843 11.89599414 12.0175164 12.139039 12.65525
3 10.7668218 11.25721357 11.3726579 11.488102 11.97849
4 10.1860204 10.65259986 10.7624385 10.872277 11.33886
5 9.6375946 10.08219388 10.1868581 10.291522 10.73612
6 9.1217734 9.54603927 9.6459167 9.745794 10.17006
7 8.6387946 9.04418136 9.1396144 9.235048 9.64043
8 8.1888978 8.57666578 8.6679512 8.759237 9.14700
9 7.7723156 8.14353686 8.2309270 8.318317 8.68954
10 7.3892646 7.74483589 7.8285418 7.912248 8.26782
11 7.0399352 7.38059913 7.4607957 7.540992 7.88166
12 6.7244802 7.05085572 7.1276886 7.204521 7.53090
13 6.4430029 6.75562533 6.8292205 6.902816 7.21544
14 6.1955428 6.49491547 6.5653915 6.635868 6.93524
15 5.9820595 6.26871848 6.3362016 6.403685 6.69034
16 5.7696526 6.04428975 6.1089428 6.173596 6.44823
17 5.4339991 5.69759119 5.7596440 5.821697 6.08529
18 5.0454361 5.29908138 5.3587927 5.418504 5.67215
19 4.6993977 4.94405130 5.0016458 5.059240 5.30389
20 4.3963458 4.63268699 4.6883247 4.743962 4.98030
21 4.1365583 4.36504142 4.4188292 4.472617 4.70110
22 3.9202312 4.14115193 4.1931594 4.245167 4.46609
23 3.7474595 3.96103662 4.0113153 4.061594 4.27517
24 3.6182953 3.82478434 3.8733944 3.922005 4.12849
25 3.5335861 3.73343196 3.7804782 3.827524 4.02737
26 3.4937186 3.68729597 3.7328665 3.778437 3.97201
27 3.4752667 3.66292175 3.7070981 3.751274 3.93893
28 3.3043525 3.48641351 3.5292729 3.572132 3.75419
29 2.9458452 3.12249549 3.1640812 3.205667 3.38232
30 2.4899112 2.66132542 2.7016785 2.742031 2.91345
31 2.3652956 2.53186083 2.5710724 2.610284 2.77685
32 2.2382402 2.40029503 2.4384448 2.476594 2.63865
33 2.0486975 2.20653724 2.2436947 2.280852 2.43869
34 2.0511798 2.20522276 2.2414864 2.277750 2.43179
35 2.0553528 2.20601792 2.2414864 2.276955 2.42762
36 2.0385642 2.18623332 2.2209965 2.255760 2.40343
37 1.8391470 1.98414706 2.0182819 2.052417 2.19742
38 1.6312788 1.77395114 1.8075380 1.841125 1.98380
39 1.5314449 1.67192652 1.7049976 1.738069 1.87855
40 1.5208780 1.65927041 1.6918497 1.724429 1.86282
41 1.4986364 1.63513027 1.6672626 1.699395 1.83589
42 1.4498027 1.58470514 1.6164629 1.648221 1.78312
43 1.2247043 1.35830771 1.3897596 1.421211 1.55481
44 1.1772885 1.30980813 1.3410049 1.372202 1.50472
45 1.1781750 1.30997706 1.3410049 1.372033 1.50383
46 1.1786125 1.31005757 1.3410014 1.371945 1.50339
47 1.1644262 1.29555858 1.3264288 1.357299 1.48843
48 1.1223208 1.25286982 1.2836027 1.314336 1.44488
49 1.0583227 1.18805529 1.2185960 1.249137 1.37887
50 1.0360396 1.16504088 1.1954094 1.225778 1.35478
51 1.0366880 1.16516444 1.1954094 1.225654 1.35413
52 0.9728290 1.10089058 1.1310379 1.161185 1.28925
53 0.6458992 0.77387319 0.8039998 0.834127 0.96210
54 0.6278378 0.75589463 0.7860408 0.816187 0.94424
55 0.6233664 0.75144260 0.7815933 0.811744 0.93982
56 0.6203139 0.74853170 0.7787158 0.808900 0.93712
57 0.4831205 0.61171664 0.6419898 0.672263 0.80086
58 0.4152141 0.54435194 0.5747526 0.605153 0.73429
59 0.4143942 0.54419570 0.5747526 0.605309 0.73511
60 0.4133407 0.54399495 0.5747526 0.605510 0.73616
61 0.3912541 0.52305164 0.5540784 0.585105 0.71690
62 0.3615872 0.49479624 0.5261553 0.557514 0.69072
63 0.3595156 0.49440150 0.5261553 0.557909 0.69279
64 0.3572502 0.49396981 0.5261553 0.558341 0.69506
65 0.3545874 0.49346241 0.5261553 0.558848 0.69772
66 0.3515435 0.49288238 0.5261553 0.559428 0.70077
67 0.3482098 0.49224713 0.5261553 0.560063 0.70410
68 0.3447026 0.49157882 0.5261553 0.560732 0.70761
69 0.3265062 0.47651151 0.5118246 0.547138 0.69714
70 0.2579257 0.41132297 0.4474346 0.483546 0.63694
71 0.2081857 0.36515737 0.4021105 0.439064 0.59604
72 0.1349572 0.29569526 0.3335350 0.371375 0.53211
73 0.0020438 0.16674762 0.2055209 0.244294 0.40900
74 -0.0243664 0.14460810 0.1843868 0.224166 0.39314
75 -0.0362635 0.13720915 0.1780468 0.218884 0.39236
76 -0.0421115 0.13609478 0.1780468 0.219999 0.39820
77 -0.0482083 0.13493301 0.1780468 0.221161 0.40430
78 -0.0546034 0.13371440 0.1780468 0.222379 0.41070
79 -0.0610386 0.13248816 0.1780468 0.223605 0.41713
80 -0.0674722 0.13126221 0.1780468 0.224831 0.42357
81 -0.0740291 0.13001276 0.1780468 0.226081 0.43012
82 -0.0809567 0.12869267 0.1780468 0.227401 0.43705
83 -0.0885308 0.12724941 0.1780468 0.228844 0.44462
84 -0.0966886 0.12569491 0.1780468 0.230399 0.45278
85 -0.1053882 0.12403716 0.1780468 0.232056 0.46148
86 -0.1147206 0.12225885 0.1780468 0.233835 0.47081
87 -0.1248842 0.12032213 0.1780468 0.235771 0.48098
88 -0.1360096 0.11820215 0.1780468 0.237891 0.49210
89 -0.1480747 0.11590310 0.1780468 0.240190 0.50417
90 -0.1611528 0.11337745 0.1780053 0.242633 0.51716
91 -0.1772967 0.10838384 0.1756366 0.242889 0.52857
92 -0.1976403 0.09964452 0.1696291 0.239614 0.53690
93 -0.2221958 0.08715720 0.1599828 0.232808 0.54216
94 -0.2510614 0.07090314 0.1466976 0.222492 0.54446
95 -0.2844042 0.05085051 0.1297736 0.208697 0.54395
96 -0.3224450 0.02695723 0.1092109 0.191465 0.54087
97 -0.3654434 -0.00082617 0.0850093 0.170845 0.53546
98 -0.4136843 -0.03255395 0.0571689 0.146892 0.52802
99 -0.4674640 -0.06828261 0.0256897 0.119662 0.51884
100 -0.5270786 -0.10806856 -0.0094284 0.089212 0.50822
knots :
[1] -2.557 -1.812 -1.726 -1.384 -1.233 -1.082 -1.046 -1.009 -0.932 -0.902
[11] -0.877 -0.838 -0.813 -0.765 -0.707 -0.665 -0.568 -0.498 -0.460 -0.413
[21] -0.347 -0.333 -0.299 -0.274 -0.226 -0.089 -0.024 -0.011 0.063 0.094
[31] 0.118 0.136 0.231 0.285 0.328 0.392 0.460 0.473 0.517 0.551
[41] 0.602 0.623 0.692 0.715 0.742 0.787 0.812 0.892 0.934 0.988
[51] 1.070 1.162 1.178 1.276 1.402 1.655 1.877 1.988 2.047 2.573
coef :
[1] 12.6970155 7.6878537 6.0937652 4.3540061 3.7259911 3.7259911
[7] 2.7408131 2.5727608 2.5727608 2.2478639 2.2414864 2.2414864
[13] 2.2414864 2.2414864 2.2414864 1.9875889 1.6964374 1.6964374
[19] 1.6623718 1.6623718 1.3410049 1.3410049 1.3410049 1.3410049
[25] 1.3410049 1.3410049 1.1954094 1.1954094 1.1954094 1.1954094
[31] 0.9829296 0.8091342 0.7815933 0.7815933 0.7815933 0.5747526
[37] 0.5747526 0.5747526 0.5747526 0.5747526 0.5261553 0.5261553
[43] 0.5261553 0.5261553 0.5261553 0.5261553 0.5261553 0.5261553
[49] 0.5261553 0.5261553 0.4273578 0.3741431 0.2060752 0.1780468
[55] 0.1780468 0.1780468 0.1780468 0.1780468 0.1780468 0.1780468
[61] -0.0094285 432.6957871
> 1 - sum(cxy3 $ resid ^ 2) / sum((y - mean(y))^2) # R^2 = 98.36%
[1] 0.96502
> showProc.time()
Time (user system elapsed): 0.282 0.006 0.649
>
> cpuTime(cxy4 <- cobs(x,y, "decrease", lambda = 1e-6, nknots = 100))# ~ 3 sec.
Time elapsed: 0.12
> 1 - sum(cxy4 $ resid ^ 2) / sum((y - mean(y))^2) # R^2 = 98.443%
[1] 0.96603
>
> cpuTime(cxy5 <- cobs(x,y, "decrease", lambda = 1e-6, nknots = 150))# ~ 8.7 sec.
Time elapsed: 0.152
> 1 - sum(cxy5 $ resid ^ 2) / sum((y - mean(y))^2) # R^2 = 98.4396%
[1] 0.96835
> showProc.time()
Time (user system elapsed): 0.604 0.016 1.651
>
>
> ## regularly spaced x :
> X <- seq(-1,1, len = 201)
> xx <- c(seq(-1.1, -1, len = 11), X,
+ seq( 1, 1.1, len = 11))
> y <- (fx <- exp(-X)) + rt(201,4)/4
> summaryCobs(cXy <- cobs(X,y, "decrease"))
qbsks2():
Performing general knot selection ...
Deleting unnecessary knots ...
List of 24
$ call : language cobs(x = X, y = y, constraint = "decrease")
$ tau : num 0.5
$ degree : num 2
$ constraint : chr "decrease"
$ ic : chr "AIC"
$ pointwise : NULL
$ select.knots : logi TRUE
$ select.lambda: logi FALSE
$ x : num [1:201] -1 -0.99 -0.98 -0.97 -0.96 -0.95 -0.94 -0.93 -0.92 -0.91 ...
$ y : num [1:201] 2.67 2.77 3.46 3.14 1.79 ...
$ resid : num [1:201] 0 0.125 0.84 0.555 -0.77 ...
$ fitted : num [1:201] 2.67 2.64 2.62 2.59 2.56 ...
$ coef : num [1:4] 2.672 1.556 0.7 0.356
$ knots : num [1:3] -1 -0.2 1
$ k0 : num 4
$ k : num 4
$ x.ps :Formal class 'matrix.csr' [package "SparseM"] with 4 slots
$ SSy : num 100
$ lambda : num 0
$ icyc : int 9
$ ifl : int 1
$ pp.lambda : NULL
$ pp.sic : NULL
$ i.mask : NULL
cb.lo ci.lo fit ci.up cb.up
1 2.46750 2.55064 2.67153 2.79242 2.87556
2 2.42251 2.50122 2.61568 2.73013 2.80884
3 2.37783 2.45240 2.56081 2.66923 2.74379
4 2.33345 2.40414 2.50694 2.60973 2.68043
5 2.28933 2.35645 2.45404 2.55164 2.61876
6 2.24548 2.30932 2.40214 2.49496 2.55879
7 2.20189 2.26274 2.35122 2.43970 2.50055
8 2.15855 2.21672 2.30129 2.38586 2.44402
9 2.11547 2.17124 2.25234 2.33344 2.38922
10 2.07265 2.12633 2.20438 2.28244 2.33611
11 2.03013 2.08199 2.15741 2.23283 2.28470
12 1.98791 2.03824 2.11142 2.18461 2.23494
13 1.94605 1.99510 2.06642 2.13775 2.18680
14 1.90459 1.95260 2.02241 2.09222 2.14023
15 1.86359 1.91078 1.97938 2.04799 2.09517
16 1.82311 1.86966 1.93734 2.00502 2.05157
17 1.78322 1.82929 1.89629 1.96328 2.00936
18 1.74397 1.78971 1.85622 1.92273 1.96847
19 1.70544 1.75096 1.81714 1.88332 1.92883
20 1.66769 1.71307 1.77904 1.84502 1.89039
21 1.63079 1.67608 1.74193 1.80779 1.85308
22 1.59478 1.64002 1.70581 1.77160 1.81684
23 1.55972 1.60493 1.67067 1.73642 1.78163
24 1.52564 1.57083 1.63653 1.70222 1.74741
25 1.49260 1.53773 1.60336 1.66899 1.71412
26 1.46062 1.50567 1.57118 1.63670 1.68175
27 1.42972 1.47466 1.53999 1.60533 1.65026
28 1.39994 1.44470 1.50979 1.57488 1.61964
29 1.37128 1.41581 1.48057 1.54533 1.58987
30 1.34375 1.38800 1.45234 1.51668 1.56093
31 1.31736 1.36126 1.42510 1.48893 1.53283
32 1.29211 1.33560 1.39884 1.46207 1.50556
33 1.26800 1.31101 1.37357 1.43612 1.47914
34 1.24500 1.28749 1.34928 1.41107 1.45356
35 1.22310 1.26502 1.32598 1.38694 1.42886
36 1.20228 1.24360 1.30367 1.36374 1.40505
37 1.18250 1.22319 1.28234 1.34150 1.38218
38 1.16372 1.20377 1.26200 1.32023 1.36028
39 1.14589 1.18532 1.24265 1.29998 1.33941
40 1.12894 1.16779 1.22428 1.28077 1.31962
41 1.11271 1.15106 1.20683 1.26259 1.30094
42 1.09639 1.13439 1.18963 1.24488 1.28287
43 1.07982 1.11760 1.17253 1.22747 1.26525
44 1.06303 1.10072 1.15553 1.21034 1.24803
45 1.04607 1.08378 1.13862 1.19346 1.23117
46 1.02898 1.06681 1.12181 1.17681 1.21463
47 1.01180 1.04982 1.10509 1.16037 1.19838
48 0.99458 1.03284 1.08847 1.14411 1.18237
49 0.97734 1.01589 1.07195 1.12801 1.16656
50 0.96011 0.99899 1.05552 1.11205 1.15092
51 0.94294 0.98216 1.03919 1.09621 1.13543
52 0.92585 0.96541 1.02295 1.08049 1.12005
53 0.90885 0.94877 1.00681 1.06485 1.10477
54 0.89197 0.93223 0.99076 1.04930 1.08956
55 0.87523 0.91581 0.97482 1.03382 1.07440
56 0.85865 0.89952 0.95896 1.01840 1.05928
57 0.84223 0.88337 0.94321 1.00304 1.04419
58 0.82598 0.86736 0.92755 0.98773 1.02911
59 0.80991 0.85150 0.91198 0.97246 1.01405
60 0.79403 0.83579 0.89651 0.95723 0.99899
61 0.77834 0.82023 0.88114 0.94205 0.98394
62 0.76284 0.80482 0.86586 0.92690 0.96888
63 0.74753 0.78956 0.85068 0.91180 0.95383
64 0.73241 0.77446 0.83559 0.89673 0.93878
65 0.71747 0.75950 0.82060 0.88171 0.92374
66 0.70271 0.74468 0.80571 0.86674 0.90871
67 0.68812 0.73001 0.79091 0.85182 0.89371
68 0.67368 0.71546 0.77621 0.83696 0.87874
69 0.65939 0.70104 0.76161 0.82217 0.86382
70 0.64523 0.68674 0.74710 0.80745 0.84896
71 0.63118 0.67254 0.73268 0.79282 0.83419
72 0.61722 0.65844 0.71836 0.77829 0.81951
73 0.60333 0.64441 0.70414 0.76388 0.80495
74 0.58948 0.63045 0.69002 0.74958 0.79055
75 0.57565 0.61654 0.67599 0.73544 0.77632
76 0.56181 0.60266 0.66205 0.72145 0.76230
77 0.54792 0.58879 0.64821 0.70764 0.74851
78 0.53395 0.57491 0.63447 0.69403 0.73500
79 0.51986 0.56100 0.62083 0.68065 0.72179
80 0.50563 0.54705 0.60728 0.66750 0.70892
81 0.49121 0.53302 0.59382 0.65462 0.69643
82 0.47657 0.51891 0.58046 0.64202 0.68435
83 0.46169 0.50468 0.56720 0.62972 0.67271
84 0.44652 0.49033 0.55403 0.61774 0.66155
85 0.43105 0.47584 0.54096 0.60609 0.65087
86 0.41526 0.46119 0.52799 0.59478 0.64072
87 0.39912 0.44638 0.51511 0.58383 0.63109
88 0.38264 0.43141 0.50233 0.57324 0.62202
89 0.36579 0.41626 0.48964 0.56302 0.61349
90 0.34858 0.40093 0.47705 0.55317 0.60552
91 0.33101 0.38542 0.46455 0.54368 0.59810
92 0.31307 0.36975 0.45215 0.53456 0.59123
93 0.29478 0.35390 0.43985 0.52580 0.58492
94 0.27615 0.33788 0.42764 0.51741 0.57914
95 0.25717 0.32170 0.41553 0.50936 0.57389
96 0.23787 0.30536 0.40352 0.50167 0.56917
97 0.21824 0.28888 0.39160 0.49431 0.56495
98 0.19830 0.27225 0.37977 0.48730 0.56125
99 0.17806 0.25547 0.36804 0.48062 0.55803
100 0.15752 0.23857 0.35641 0.47426 0.55531
knots :
[1] -1.0 -0.2 1.0
coef :
[1] 2.67153 1.55592 0.70045 0.35641
> 1 - sum(cXy $ resid ^ 2) / sum((y - mean(y))^2) # R^2 = 77.2%
[1] 0.77644
> showProc.time()
Time (user system elapsed): 0.177 0.011 0.51
>
> (cXy.9 <- cobs(X,y, "decrease", tau = 0.9))
qbsks2():
Performing general knot selection ...
Deleting unnecessary knots ...
COBS regression spline (degree = 2) from call:
cobs(x = X, y = y, constraint = "decrease", tau = 0.9)
{tau=0.9}-quantile; dimensionality of fit: 6 from {6}
x$knots[1:5]: -1.0, -0.6, -0.2, 0.2, 1.0
> (cXy.1 <- cobs(X,y, "decrease", tau = 0.1))
qbsks2():
Performing general knot selection ...
WARNING! Since the number of 6 knots selected by AIC reached the
upper bound during general knot selection, you might want to rerun
cobs with a larger number of knots.
Deleting unnecessary knots ...
WARNING! Since the number of 6 knots selected by AIC reached the
upper bound during general knot selection, you might want to rerun
cobs with a larger number of knots.
COBS regression spline (degree = 2) from call:
cobs(x = X, y = y, constraint = "decrease", tau = 0.1)
{tau=0.1}-quantile; dimensionality of fit: 4 from {4}
x$knots[1:3]: -1.0, 0.6, 1.0
> (cXy.99<- cobs(X,y, "decrease", tau = 0.99))
qbsks2():
Performing general knot selection ...
Deleting unnecessary knots ...
COBS regression spline (degree = 2) from call:
cobs(x = X, y = y, constraint = "decrease", tau = 0.99)
{tau=0.99}-quantile; dimensionality of fit: 4 from {4}
x$knots[1:3]: -1.0, -0.2, 1.0
> (cXy.01<- cobs(X,y, "decrease", tau = 0.01))
qbsks2():
Performing general knot selection ...
Deleting unnecessary knots ...
COBS regression spline (degree = 2) from call:
cobs(x = X, y = y, constraint = "decrease", tau = 0.01)
{tau=0.01}-quantile; dimensionality of fit: 6 from {6}
x$knots[1:5]: -1.0, -0.6, -0.2, 0.2, 1.0
> plot(X,y, xlim = range(xx),
+ main = "cobs(*, \"decrease\"), N=201, tau = 50% (Med.), 1,10, 90,99%")
> lines(predict(cXy, xx), col = 2)
> lines(predict(cXy.1, xx), col = 3)
> lines(predict(cXy.9, xx), col = 3)
> lines(predict(cXy.01, xx), col = 4)
> lines(predict(cXy.99, xx), col = 4)
>
> showProc.time()
Time (user system elapsed): 0.873 0.007 2.057
>
> ## Interpolation
> cpuTime(cXyI <- cobs(X,y, "decrease", knots = unique(X)))
qbsks2():
Performing general knot selection ...
Error in x %*% coefficients : NA/NaN/Inf in foreign function call (arg 2)
Calls: cpuTime ... cobs -> qbsks2 -> drqssbc2 -> rq.fit.sfnc -> %*% -> %*%
In addition: Warning message:
In cobs(X, y, "decrease", knots = unique(X)) :
The number of knots can't be equal to the number of unique x for degree = 2.
'cobs' has automatically deleted the middle knot.
Timing stopped at: 1.287 0.016 3.571
Execution halted
Running the tests in ‘tests/multi-constr.R’ failed.
Complete output:
> #### 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(), ...
> Lnx <- Sys.info()[["sysname"]] == "Linux"
> isMac <- Sys.info()[["sysname"]] == "Darwin"
> x86 <- (arch <- Sys.info()[["machine"]]) == "x86_64"
> noLdbl <- (.Machine$sizeof.longdouble <= 8) ## TRUE when --disable-long-double
> ## IGNORE_RDIFF_BEGIN
> Sys.info()
sysname
"Linux"
release
"6.2.15-100.fc36.x86_64"
version
"#1 SMP PREEMPT_DYNAMIC Thu May 11 16:51:53 UTC 2023"
nodename
"gannet.stats.ox.ac.uk"
machine
"x86_64"
login
"ripley"
user
"ripley"
effective_user
"ripley"
> noLdbl
[1] FALSE
> ## IGNORE_RDIFF_END
>
>
> Rsq <- function(obj) {
+ stopifnot(inherits(obj, "cobs"), is.numeric(res <- obj$resid))
+ 1 - sum(res^2)/obj$SSy
+ }
> list_ <- function (...) `names<-`(list(...), vapply(sys.call()[-1L], as.character, ""))
> is.cobs <- function(x) inherits(x, "cobs")
>
> 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
> ## IGNORE_RDIFF_BEGIN
> r2i$warning
<simpleWarning in cobs(x, y, constraint = "increase", pointwise = con): drqssbc2(): Not all flags are normal (== 1), ifl : 21>
> ## IGNORE_RDIFF_END
> ## 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%)
loo.design2(): -> Xeq 50 x 3 (nz = 150 =^= 1%)
Xieq 3 x 3 (nz = 9 =^= 1%)
>
> ## now gives warning (not error): but no warning on M1 mac -> IGNORE
> ## IGNORE_RDIFF_BEGIN
> 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
> ## IGNORE_RDIFF_END
> 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%)
Warning message:
In cobs(x, y, degree = 1, pointwise = con, trace = 3) :
drqssbc2(): Not all flags are normal (== 1), ifl : 19
> 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'
>
> ## Named list of all cobs() results above -- sort() collation order matters for ls() !
> (curLC <- Sys.getlocale("LC_COLLATE"))
[1] "C"
> Sys.setlocale("LC_COLLATE", "C")
[1] "C"
> cobsL <- mget(Filter(\(nm) is.cobs(.GlobalEnv[[nm]]), ls(patt="c[12p]")),
+ envir = .GlobalEnv)
> Sys.setlocale("LC_COLLATE", curLC) # reverting
[1] "C"
>
> knL <- lapply(cobsL, `[[`, "knots")
> str(knL[order(lengths(knL))])
List of 16
$ c2IC : num [1:2] -1 2
$ cp2IC: num [1:2] -1 2
$ cp2c : num [1:2] -1 2
$ c1IC : num [1:3] -1 0.776 2
$ c1c : num [1:3] -1 0.776 2
$ c2 : num [1:3] -1 -0.449 2
$ c2c : num [1:3] -1 0.163 2
$ cp1IC: num [1:3] -1 0.776 2
$ cp1c : num [1:3] -1 0.776 2
$ c2i : num [1:4] -1 0.163 0.776 2
$ cp2 : num [1:4] -1 -0.449 0.776 2
$ cp2i : num [1:4] -1 0.163 0.776 2
$ c1 : num [1:5] -1 -0.449 0.163 0.776 2
$ c1i : num [1:5] -1 -0.449 0.163 0.776 2
$ cp1 : num [1:5] -1 -0.449 0.163 0.776 2
$ cp1i : num [1:5] -1 -0.449 0.163 0.776 2
>
> gotRsqrs <- sapply(cobsL, Rsq)
> Rsqrs <- c(c1 = 0.95079126, c1IC = 0.92974549, c1c = 0.92974549, c1i = 0.95079126,
+ c2 = 0.94637437, c2IC = 0.91375404, c2c = 0.92505977, c2i = 0.95022829,
+ cp1 = 0.9426453, cp1IC = 0.92223149, cp1c = 0.92223149, cp1i = 0.9426453,
+ cp2 = 0.94988863, cp2IC= 0.90051964, cp2c = 0.91375409, cp2i = 0.93611487)
> ## M1 mac " = " , cp2IC= 0.91704726, " = " , cp2i = 0.94620178
> ## noLD " = " , cp2IC=-0.08244284, " = " , cp2i = 0.94636815
> ## ATLAS " = " , cp2IC= 0.91471729, " = " , cp2i = 0.94506339
> ## openBLAS " = " , cp2IC= 0.91738019, " = " , cp2i = 0.93589404
> ## MKL " = " , cp2IC= 0.91765403, " = " , cp2i = 0.94501205
> ## Intel " = " , cp2IC= 0.91765403, " = " , cp2i = 0.94501205
> ## ^^^^^^^^^^ ^^^^^^^^^^
> ## remove these two from testing, notably for the M1 Mac & noLD .. :
> ##iR2 <- if(!x86 || noLdbl) setdiff(names(cobsL), c("cp2IC", "cp2i")) else TRUE
> ## actually everywhere, because of ATLAS, openBLAS, MKL, Intel... :
> iR2 <- setdiff(names(cobsL), nR2 <- c("cp2IC", "cp2i"))
> ## IGNORE_RDIFF_BEGIN
> dput(signif(gotRsqrs, digits=8))
c(c1 = 0.95079126, c1IC = 0.92974549, c1c = 0.92974549, c1i = 0.95079126,
c2 = 0.94637437, c2IC = 0.91375404, c2c = 0.92505977, c2i = 0.95022829,
cp1 = 0.95318954, cp1IC = 0.92223149, cp1c = 0.92223149, cp1i = 0.9426453,
cp2 = 0.94988863, cp2IC = 0.90051964, cp2c = 0.91375409, cp2i = 0.94404129
)
> all.equal(Rsqrs[iR2], gotRsqrs[iR2], tolerance=0)# 2.6277e-9 (Lnx F 38); 2.6898e-9 (M1 mac)
[1] "Mean relative difference: 0.000805528"
> all.equal(Rsqrs[nR2], gotRsqrs[nR2], tolerance=0)# differ; drastically only for 'noLD'
[1] "Mean relative difference: 0.00431573"
> ## IGNORE_RDIFF_END
> stopifnot(exprs = {
+ all.equal(Rsqrs[iR2], gotRsqrs[iR2])
+ identical(c(5L, 3L, 3L, 5L,
+ 3L, 2L, 3L, 4L,
+ 5L, 3L, 3L, 5L,
+ 4L, 2L, 2L, 4L), unname(lengths(knL)))
+ })
Error: Rsqrs[iR2] and gotRsqrs[iR2] are not equal:
Mean relative difference: 0.000805528
Execution halted
Running the tests in ‘tests/roof.R’ failed.
Complete output:
> suppressMessages(library(cobs))
>
> data(USArmyRoofs)
> attach(USArmyRoofs)#-> "age" and "fci"
>
> if(!dev.interactive(orNone=TRUE)) pdf("roof.pdf", width=10)
>
> ## Compute the quadratic median smoothing B-spline with SIC
> ## chosen lambda
> a50 <- cobs(age,fci,constraint = "decrease",lambda = -1,nknots = 10,
+ degree = 2,pointwise = rbind(c(0,0,100)),
+ trace = 2)# trace > 1 : more tracing
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.
fieq=TRUE -> Tnobs = 184, n0 = 29, |ptConstr| = 2
Error in drqssbc2(x, y, w, pw = pw, knots = knots, degree = degree, Tlambda = if (select.lambda) lambdaSet else lambda, :
The problem is degenerate for the range of lambda specified.
Calls: cobs -> drqssbc2
In addition: Warning message:
In min(sol1["k", i.keep]) : no non-missing arguments to min; returning Inf
Execution halted
Running the tests in ‘tests/wind.R’ failed.
Complete output:
> 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.002 0 0.007
>
> 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.692 0.033 2.284
>
> 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.05 0.003 0.136
>
> ## 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): 2.921 0.031 8.021
>
> 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.
Error in x %*% coefficients : NA/NaN/Inf in foreign function call (arg 2)
Calls: cobs -> drqssbc2 -> rq.fit.sfnc -> %*% -> %*%
Execution halted
- checking PDF version of manual ... [8s/23s] OK
- checking HTML version of manual ... OK
- checking for non-standard things in the check directory ... OK
- checking for detritus in the temp directory ... OK
- DONE
Status: 1 ERROR