Last updated on 2024-11-24 15:49:22 CET.
Flavor | Version | Tinstall | Tcheck | Ttotal | Status | Flags |
---|---|---|---|---|---|---|
r-devel-linux-x86_64-debian-clang | 1.3-8 | 10.99 | 117.85 | 128.84 | ERROR | |
r-devel-linux-x86_64-debian-gcc | 1.3-8 | 7.71 | 86.23 | 93.94 | ERROR | |
r-devel-linux-x86_64-fedora-clang | 1.3-8 | 219.92 | ERROR | |||
r-devel-linux-x86_64-fedora-gcc | 1.3-8 | 215.39 | ERROR | |||
r-devel-windows-x86_64 | 1.3-8 | 15.00 | 129.00 | 144.00 | ERROR | |
r-patched-linux-x86_64 | 1.3-8 | 10.64 | 115.10 | 125.74 | OK | |
r-release-linux-x86_64 | 1.3-8 | 10.40 | 114.40 | 124.80 | OK | |
r-release-macos-arm64 | 1.3-8 | 69.00 | OK | |||
r-release-macos-x86_64 | 1.3-8 | 96.00 | OK | |||
r-release-windows-x86_64 | 1.3-8 | 14.00 | 135.00 | 149.00 | OK | |
r-oldrel-macos-arm64 | 1.3-8 | 74.00 | OK | |||
r-oldrel-macos-x86_64 | 1.3-8 | 160.00 | OK | |||
r-oldrel-windows-x86_64 | 1.3-8 | 17.00 | 166.00 | 183.00 | OK |
Version: 1.3-8
Check: tests
Result: ERROR
Running ‘0_pt-ex.R’ [2s/3s]
Running ‘ex1.R’ [4s/5s]
Running ‘ex2-long.R’ [10s/12s]
Running ‘ex3.R’ [2s/3s]
Comparing ‘ex3.Rout’ to ‘ex3.Rout.save’ ... OK
Running ‘multi-constr.R’ [4s/5s]
Running ‘roof.R’ [3s/4s]
Running ‘small-ex.R’ [3s/3s]
Comparing ‘small-ex.Rout’ to ‘small-ex.Rout.save’ ...24,25d23
< Warning message:
< In cobs(x, y) : drqssbc2(): Not all flags are normal (== 1), ifl : 22
30,33d27
< 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.
<
35,38d28
<
< 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.
41,51c31,34
<
< **** ERROR in algorithm: ifl = 23
<
<
< {tau=0.5}-quantile; dimensionality of fit: 7 from {7}
< x$knots[1:6]: 0.999989, 2.000000, 3.000000, ... , 12.000011
< coef[1:7]: 6.99997, 8.58325, 23.41667, 28.16667, 45.91667, 67.50003, 96.00010
< R^2 = 100% ; empirical tau (over all): 5/7 = 0.714286 (target tau= 0.5)
< Warning message:
< In cobs(x, y, constraint = "convex") :
< drqssbc2(): Not all flags are normal (== 1), ifl : 23
---
> {tau=0.5}-quantile; dimensionality of fit: 3 from {3}
> x$knots[1:2]: 0.999989, 12.000011
> coef[1:3]: 7.99991, 52.00000, 96.00009
> R^2 = 99.88% ; empirical tau (over all): 3/7 = 0.428571 (target tau= 0.5)
56,59d38
< 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.
<
61,64d39
<
< 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.
67,77c42,45
<
< **** ERROR in algorithm: ifl = 23
<
<
< {tau=0.5}-quantile; dimensionality of fit: 7 from {7}
< x$knots[1:6]: 0.999989, 2.000000, 3.000000, ... , 12.000011
< coef[1:7]: 6.99997, 8.58325, 23.41667, 28.16667, 45.91667, 67.50003, 96.00010
< R^2 = 100% ; empirical tau (over all): 5/7 = 0.714286 (target tau= 0.5)
< Warning message:
< In cobs(x, y, constraint = "concave") :
< drqssbc2(): Not all flags are normal (== 1), ifl : 23
---
> {tau=0.5}-quantile; dimensionality of fit: 3 from {3}
> x$knots[1:2]: 0.999989, 12.000011
> coef[1:3]: 7.54157, 54.29167, 96.00008
> R^2 = 99.84% ; empirical tau (over all): 3/7 = 0.428571 (target tau= 0.5)
80c48
< [1] 6.31089e-28
---
> [1] 7
82c50
< [1] 6.31089e-28
---
> [1] 9.71528
91,192c59,158
< [1,] 0.999989 6.99997 6.99997 6.99997 NaN NaN
< [2,] 1.111100 7.42382 7.42382 7.42382 NaN NaN
< [3,] 1.222212 7.99171 7.99171 7.99171 NaN NaN
< [4,] 1.333323 8.70363 8.70363 8.70363 NaN NaN
< [5,] 1.444434 9.55959 9.55959 9.55959 NaN NaN
< [6,] 1.555546 10.55958 10.55958 10.55958 NaN NaN
< [7,] 1.666657 11.70360 11.70360 11.70360 NaN NaN
< [8,] 1.777768 12.99166 12.99166 12.99166 NaN NaN
< [9,] 1.888880 14.42374 14.42374 14.42374 NaN NaN
< [10,] 1.999991 15.99987 15.99987 15.99987 NaN NaN
< [11,] 2.111102 17.57601 17.57601 17.57601 NaN NaN
< [12,] 2.222214 19.00813 19.00813 19.00813 NaN NaN
< [13,] 2.333325 20.29621 20.29621 20.29621 NaN NaN
< [14,] 2.444436 21.44025 21.44025 21.44025 NaN NaN
< [15,] 2.555548 22.44026 22.44026 22.44026 NaN NaN
< [16,] 2.666659 23.29624 23.29624 23.29624 NaN NaN
< [17,] 2.777770 24.00819 24.00819 24.00819 NaN NaN
< [18,] 2.888882 24.57610 24.57610 24.57610 NaN NaN
< [19,] 2.999993 24.99998 24.99998 24.99998 NaN NaN
< [20,] 3.111104 25.37858 25.37858 25.37858 NaN NaN
< [21,] 3.222216 25.81067 25.81067 25.81067 NaN NaN
< [22,] 3.333327 26.29627 26.29627 26.29627 NaN NaN
< [23,] 3.444438 26.83536 26.83536 26.83536 NaN NaN
< [24,] 3.555550 27.42795 27.42795 27.42795 NaN NaN
< [25,] 3.666661 28.07404 28.07404 28.07404 NaN NaN
< [26,] 3.777772 28.77363 28.77363 28.77363 NaN NaN
< [27,] 3.888884 29.52671 29.52671 29.52671 NaN NaN
< [28,] 3.999995 30.33330 30.33330 30.33330 NaN NaN
< [29,] 4.111106 31.19338 31.19338 31.19338 NaN NaN
< [30,] 4.222218 32.10696 32.10696 32.10696 NaN NaN
< [31,] 4.333329 33.07404 33.07404 33.07404 NaN NaN
< [32,] 4.444440 34.09461 34.09461 34.09461 NaN NaN
< [33,] 4.555552 35.16869 35.16869 35.16869 NaN NaN
< [34,] 4.666663 36.29626 36.29626 36.29626 NaN NaN
< [35,] 4.777774 37.47733 37.47733 37.47733 NaN NaN
< [36,] 4.888886 38.71190 38.71190 38.71190 NaN NaN
< [37,] 4.999997 39.99996 39.99996 39.99996 NaN NaN
< [38,] 5.111108 41.27980 41.27980 41.27980 NaN NaN
< [39,] 5.222220 42.48969 42.48969 42.48969 NaN NaN
< [40,] 5.333331 43.62961 43.62961 43.62961 NaN NaN
< [41,] 5.444442 44.69957 44.69957 44.69957 NaN NaN
< [42,] 5.555554 45.69957 45.69957 45.69957 NaN NaN
< [43,] 5.666665 46.62962 46.62962 46.62962 NaN NaN
< [44,] 5.777776 47.48970 47.48970 47.48970 NaN NaN
< [45,] 5.888888 48.27983 48.27983 48.27983 NaN NaN
< [46,] 5.999999 48.99999 48.99999 48.99999 NaN NaN
< [47,] 6.111110 49.68861 49.68861 49.68861 NaN NaN
< [48,] 6.222222 50.38408 50.38408 50.38408 NaN NaN
< [49,] 6.333333 51.08642 51.08642 51.08642 NaN NaN
< [50,] 6.444444 51.79561 51.79561 51.79561 NaN NaN
< [51,] 6.555556 52.51166 52.51166 52.51166 NaN NaN
< [52,] 6.666667 53.23457 53.23457 53.23457 NaN NaN
< [53,] 6.777778 53.96434 53.96434 53.96434 NaN NaN
< [54,] 6.888890 54.70097 54.70097 54.70097 NaN NaN
< [55,] 7.000001 55.44445 55.44445 55.44445 NaN NaN
< [56,] 7.111112 56.19480 56.19480 56.19480 NaN NaN
< [57,] 7.222224 56.95200 56.95200 56.95200 NaN NaN
< [58,] 7.333335 57.71606 57.71606 57.71606 NaN NaN
< [59,] 7.444446 58.48698 58.48698 58.48698 NaN NaN
< [60,] 7.555558 59.26476 59.26476 59.26476 NaN NaN
< [61,] 7.666669 60.04940 60.04940 60.04940 NaN NaN
< [62,] 7.777780 60.84090 60.84090 60.84090 NaN NaN
< [63,] 7.888892 61.63925 61.63925 61.63925 NaN NaN
< [64,] 8.000003 62.44447 62.44447 62.44447 NaN NaN
< [65,] 8.111114 63.25654 63.25654 63.25654 NaN NaN
< [66,] 8.222226 64.07547 64.07547 64.07547 NaN NaN
< [67,] 8.333337 64.90126 64.90126 64.90126 NaN NaN
< [68,] 8.444448 65.73391 65.73391 65.73391 NaN NaN
< [69,] 8.555560 66.57342 66.57342 66.57342 NaN NaN
< [70,] 8.666671 67.41979 67.41979 67.41979 NaN NaN
< [71,] 8.777782 68.27301 68.27301 68.27301 NaN NaN
< [72,] 8.888894 69.13310 69.13310 69.13310 NaN NaN
< [73,] 9.000005 70.00004 70.00004 70.00004 NaN NaN
< [74,] 9.111116 70.87384 70.87384 70.87384 NaN NaN
< [75,] 9.222228 71.75450 71.75450 71.75450 NaN NaN
< [76,] 9.333339 72.64202 72.64202 72.64202 NaN NaN
< [77,] 9.444450 73.53640 73.53640 73.53640 NaN NaN
< [78,] 9.555562 74.43763 74.43763 74.43763 NaN NaN
< [79,] 9.666673 75.34573 75.34573 75.34573 NaN NaN
< [80,] 9.777784 76.26068 76.26068 76.26068 NaN NaN
< [81,] 9.888896 77.18250 77.18250 77.18250 NaN NaN
< [82,] 10.000007 78.11117 78.11117 78.11117 NaN NaN
< [83,] 10.111118 79.04670 79.04670 79.04670 NaN NaN
< [84,] 10.222230 79.98909 79.98909 79.98909 NaN NaN
< [85,] 10.333341 80.93834 80.93834 80.93834 NaN NaN
< [86,] 10.444452 81.89444 81.89444 81.89444 NaN NaN
< [87,] 10.555564 82.85741 82.85741 82.85741 NaN NaN
< [88,] 10.666675 83.82723 83.82723 83.82723 NaN NaN
< [89,] 10.777786 84.80392 84.80392 84.80392 NaN NaN
< [90,] 10.888898 85.78746 85.78746 85.78746 NaN NaN
< [91,] 11.000009 86.77786 86.77786 86.77786 NaN NaN
< [92,] 11.111120 87.77512 87.77512 87.77512 NaN NaN
< [93,] 11.222232 88.77923 88.77923 88.77923 NaN NaN
< [94,] 11.333343 89.79021 89.79021 89.79021 NaN NaN
< [95,] 11.444454 90.80805 90.80805 90.80805 NaN NaN
< [96,] 11.555566 91.83274 91.83274 91.83274 NaN NaN
< [97,] 11.666677 92.86429 92.86429 92.86429 NaN NaN
< [98,] 11.777788 93.90270 93.90270 93.90270 NaN NaN
< [99,] 11.888900 94.94797 94.94797 94.94797 NaN NaN
< [100,] 12.000011 96.00010 96.00010 96.00010 NaN NaN
< Warning message:
< In qt((1 + level)/2, n - object$k) : NaNs produced
---
> [1,] 0.999989 7.99991 4.82789 11.1719 4.84949 11.1503
> [2,] 1.111100 8.88880 5.84444 11.9332 5.86518 11.9124
> [3,] 1.222212 9.77769 6.85519 12.7002 6.87510 12.6803
> [4,] 1.333323 10.66658 7.85994 13.4732 7.87905 13.4541
> [5,] 1.444434 11.55548 8.85846 14.2525 8.87683 14.2341
> [6,] 1.555546 12.44437 9.85054 15.0382 9.86821 15.0205
> [7,] 1.666657 13.33326 10.83596 15.8305 10.85297 15.8135
> [8,] 1.777768 14.22215 11.81451 16.6298 11.83091 16.6134
> [9,] 1.888880 15.11104 12.78598 17.4361 12.80181 17.4203
> [10,] 1.999991 15.99993 13.75019 18.2497 13.76551 18.2343
> [11,] 2.111102 16.88882 14.70699 19.0706 14.72185 19.0558
> [12,] 2.222214 17.77771 15.65628 19.8991 15.67072 19.8847
> [13,] 2.333325 18.66660 16.59799 20.7352 16.61208 20.7211
> [14,] 2.444436 19.55549 17.53215 21.5788 17.54593 21.5651
> [15,] 2.555548 20.44438 18.45883 22.4299 18.47235 22.4164
> [16,] 2.666659 21.33327 19.37817 23.2884 19.39149 23.2751
> [17,] 2.777770 22.22216 20.29043 24.1539 20.30359 24.1407
> [18,] 2.888882 23.11105 21.19590 25.0262 21.20894 25.0132
> [19,] 2.999993 23.99994 22.09496 25.9049 22.10793 25.8920
> [20,] 3.111104 24.88884 22.98804 26.7896 23.00099 26.7767
> [21,] 3.222216 25.77773 23.87563 27.6798 23.88858 27.6669
> [22,] 3.333327 26.66662 24.75824 28.5750 24.77123 28.5620
> [23,] 3.444438 27.55551 25.63639 29.4746 25.64946 29.4616
> [24,] 3.555550 28.44440 26.51062 30.3782 26.52379 30.3650
> [25,] 3.666661 29.33329 27.38147 31.2851 27.39476 31.2718
> [26,] 3.777772 30.22218 28.24944 32.1949 28.26287 32.1815
> [27,] 3.888884 31.11107 29.11501 33.1071 29.12861 33.0935
> [28,] 3.999995 31.99996 29.97866 34.0213 29.99243 34.0075
> [29,] 4.111106 32.88885 30.84080 34.9369 30.85475 34.9230
> [30,] 4.222218 33.77774 31.70184 35.8536 31.71597 35.8395
> [31,] 4.333329 34.66663 32.56212 36.7711 32.57645 36.7568
> [32,] 4.444440 35.55552 33.42197 37.6891 33.43650 37.6745
> [33,] 4.555552 36.44441 34.28170 38.6071 34.29643 38.5924
> [34,] 4.666663 37.33330 35.14155 39.5251 35.15648 39.5101
> [35,] 4.777774 38.22219 36.00178 40.4426 36.01690 40.4275
> [36,] 4.888886 39.11109 36.86257 41.3596 36.87789 41.3443
> [37,] 4.999997 39.99998 37.72413 42.2758 37.73963 42.2603
> [38,] 5.111108 40.88887 38.58661 43.1911 38.60229 43.1754
> [39,] 5.222220 41.77776 39.45015 44.1054 39.46601 44.0895
> [40,] 5.333331 42.66665 40.31488 45.0184 40.33090 45.0024
> [41,] 5.444442 43.55554 41.18091 45.9302 41.19708 45.9140
> [42,] 5.555554 44.44443 42.04831 46.8405 42.06463 46.8242
> [43,] 5.666665 45.33332 42.91718 47.7495 42.93363 47.7330
> [44,] 5.777776 46.22221 43.78757 48.6569 43.80415 48.6403
> [45,] 5.888888 47.11110 44.65953 49.5627 44.67622 49.5460
> [46,] 5.999999 47.99999 45.53310 50.4669 45.54990 50.4501
> [47,] 6.111110 48.88888 46.40831 51.3695 46.42520 51.3526
> [48,] 6.222222 49.77777 47.28517 52.2704 47.30215 52.2534
> [49,] 6.333333 50.66666 48.16370 53.1696 48.18074 53.1526
> [50,] 6.444444 51.55555 49.04388 54.0672 49.06099 54.0501
> [51,] 6.555556 52.44445 49.92572 54.9632 49.94287 54.9460
> [52,] 6.666667 53.33334 50.80918 55.8575 50.82637 55.8403
> [53,] 6.777778 54.22223 51.69424 56.7502 51.71145 56.7330
> [54,] 6.888890 55.11112 52.58085 57.6414 52.59808 57.6242
> [55,] 7.000001 56.00001 53.46896 58.5311 53.48620 58.5138
> [56,] 7.111112 56.88890 54.35852 59.4193 54.37575 59.4020
> [57,] 7.222224 57.77779 55.24944 60.3061 55.26666 60.2889
> [58,] 7.333335 58.66668 56.14166 61.1917 56.15886 61.1745
> [59,] 7.444446 59.55557 57.03507 62.0761 57.05224 62.0589
> [60,] 7.555558 60.44446 57.92957 62.9593 57.94670 62.9422
> [61,] 7.666669 61.33335 58.82505 63.8417 58.84213 63.8246
> [62,] 7.777780 62.22224 59.72137 64.7231 59.73840 64.7061
> [63,] 7.888892 63.11113 60.61839 65.6039 60.63536 65.5869
> [64,] 8.000003 64.00002 61.51595 66.4841 61.53286 66.4672
> [65,] 8.111114 64.88891 62.41387 67.3640 62.43073 67.3471
> [66,] 8.222226 65.77781 63.31198 68.2436 63.32877 68.2268
> [67,] 8.333337 66.66670 64.21005 69.1233 64.22678 69.1066
> [68,] 8.444448 67.55559 65.10788 70.0033 65.12455 69.9866
> [69,] 8.555560 68.44448 66.00522 70.8837 66.02184 70.8671
> [70,] 8.666671 69.33337 66.90183 71.7649 66.91839 71.7484
> [71,] 8.777782 70.22226 67.79742 72.6471 67.81393 72.6306
> [72,] 8.888894 71.11115 68.69171 73.5306 68.70819 73.5141
> [73,] 9.000005 72.00004 69.58441 74.4157 69.60086 74.3992
> [74,] 9.111116 72.88893 70.47520 75.3027 70.49164 75.2862
> [75,] 9.222228 73.77782 71.36376 76.1919 71.38020 76.1754
> [76,] 9.333339 74.66671 72.24976 77.0837 72.26622 77.0672
> [77,] 9.444450 75.55560 73.13285 77.9783 73.14935 77.9618
> [78,] 9.555562 76.44449 74.01272 78.8763 74.02928 78.8597
> [79,] 9.666673 77.33338 74.88902 79.7777 74.90567 79.7611
> [80,] 9.777784 78.22227 75.76143 80.6831 75.77819 80.6664
> [81,] 9.888896 79.11116 76.62965 81.5927 76.64655 81.5758
> [82,] 10.000007 80.00006 77.49337 82.5067 77.51044 82.4897
> [83,] 10.111118 80.88895 78.35232 83.4256 78.36960 83.4083
> [84,] 10.222230 81.77784 79.20626 84.3494 79.22378 84.3319
> [85,] 10.333341 82.66673 80.05498 85.2785 80.07276 85.2607
> [86,] 10.444452 83.55562 80.89827 86.2130 80.91637 86.1949
> [87,] 10.555564 84.44451 81.73600 87.1530 81.75445 87.1346
> [88,] 10.666675 85.33340 82.56803 88.0988 82.58687 88.0799
> [89,] 10.777786 86.22229 83.39429 89.0503 83.41355 89.0310
> [90,] 10.888898 87.11118 84.21470 90.0077 84.23443 89.9879
> [91,] 11.000009 88.00007 85.02924 90.9709 85.04947 90.9507
> [92,] 11.111120 88.88896 85.83791 91.9400 85.85868 91.9192
> [93,] 11.222232 89.77785 86.64072 92.9150 86.66208 92.8936
> [94,] 11.333343 90.66674 87.43771 93.8958 87.45970 93.8738
> [95,] 11.444454 91.55563 88.22894 94.8823 88.25160 94.8597
> [96,] 11.555566 92.44452 89.01448 95.8746 89.03784 95.8512
> [97,] 11.666677 93.33342 89.79440 96.8724 89.81850 96.8483
> [98,] 11.777788 94.22231 90.56880 97.8758 90.59368 97.8509
> [99,] 11.888900 95.11120 91.33776 98.8846 91.36346 98.8589
> [100,] 12.000011 96.00009 92.10139 99.8988 92.12794 99.8722
195,296c161,260
< [1,] 0.999989 6.99997 6.99997 6.99997 NaN NaN
< [2,] 1.111100 7.42382 7.42382 7.42382 NaN NaN
< [3,] 1.222212 7.99171 7.99171 7.99171 NaN NaN
< [4,] 1.333323 8.70363 8.70363 8.70363 NaN NaN
< [5,] 1.444434 9.55959 9.55959 9.55959 NaN NaN
< [6,] 1.555546 10.55958 10.55958 10.55958 NaN NaN
< [7,] 1.666657 11.70360 11.70360 11.70360 NaN NaN
< [8,] 1.777768 12.99166 12.99166 12.99166 NaN NaN
< [9,] 1.888880 14.42374 14.42374 14.42374 NaN NaN
< [10,] 1.999991 15.99987 15.99987 15.99987 NaN NaN
< [11,] 2.111102 17.57601 17.57601 17.57601 NaN NaN
< [12,] 2.222214 19.00813 19.00813 19.00813 NaN NaN
< [13,] 2.333325 20.29621 20.29621 20.29621 NaN NaN
< [14,] 2.444436 21.44025 21.44025 21.44025 NaN NaN
< [15,] 2.555548 22.44026 22.44026 22.44026 NaN NaN
< [16,] 2.666659 23.29624 23.29624 23.29624 NaN NaN
< [17,] 2.777770 24.00819 24.00819 24.00819 NaN NaN
< [18,] 2.888882 24.57610 24.57610 24.57610 NaN NaN
< [19,] 2.999993 24.99998 24.99998 24.99998 NaN NaN
< [20,] 3.111104 25.37858 25.37858 25.37858 NaN NaN
< [21,] 3.222216 25.81067 25.81067 25.81067 NaN NaN
< [22,] 3.333327 26.29627 26.29627 26.29627 NaN NaN
< [23,] 3.444438 26.83536 26.83536 26.83536 NaN NaN
< [24,] 3.555550 27.42795 27.42795 27.42795 NaN NaN
< [25,] 3.666661 28.07404 28.07404 28.07404 NaN NaN
< [26,] 3.777772 28.77363 28.77363 28.77363 NaN NaN
< [27,] 3.888884 29.52671 29.52671 29.52671 NaN NaN
< [28,] 3.999995 30.33330 30.33330 30.33330 NaN NaN
< [29,] 4.111106 31.19338 31.19338 31.19338 NaN NaN
< [30,] 4.222218 32.10696 32.10696 32.10696 NaN NaN
< [31,] 4.333329 33.07404 33.07404 33.07404 NaN NaN
< [32,] 4.444440 34.09461 34.09461 34.09461 NaN NaN
< [33,] 4.555552 35.16869 35.16869 35.16869 NaN NaN
< [34,] 4.666663 36.29626 36.29626 36.29626 NaN NaN
< [35,] 4.777774 37.47733 37.47733 37.47733 NaN NaN
< [36,] 4.888886 38.71190 38.71190 38.71190 NaN NaN
< [37,] 4.999997 39.99996 39.99996 39.99996 NaN NaN
< [38,] 5.111108 41.27980 41.27980 41.27980 NaN NaN
< [39,] 5.222220 42.48969 42.48969 42.48969 NaN NaN
< [40,] 5.333331 43.62961 43.62961 43.62961 NaN NaN
< [41,] 5.444442 44.69957 44.69957 44.69957 NaN NaN
< [42,] 5.555554 45.69957 45.69957 45.69957 NaN NaN
< [43,] 5.666665 46.62962 46.62962 46.62962 NaN NaN
< [44,] 5.777776 47.48970 47.48970 47.48970 NaN NaN
< [45,] 5.888888 48.27983 48.27983 48.27983 NaN NaN
< [46,] 5.999999 48.99999 48.99999 48.99999 NaN NaN
< [47,] 6.111110 49.68861 49.68861 49.68861 NaN NaN
< [48,] 6.222222 50.38408 50.38408 50.38408 NaN NaN
< [49,] 6.333333 51.08642 51.08642 51.08642 NaN NaN
< [50,] 6.444444 51.79561 51.79561 51.79561 NaN NaN
< [51,] 6.555556 52.51166 52.51166 52.51166 NaN NaN
< [52,] 6.666667 53.23457 53.23457 53.23457 NaN NaN
< [53,] 6.777778 53.96434 53.96434 53.96434 NaN NaN
< [54,] 6.888890 54.70097 54.70097 54.70097 NaN NaN
< [55,] 7.000001 55.44445 55.44445 55.44445 NaN NaN
< [56,] 7.111112 56.19480 56.19480 56.19480 NaN NaN
< [57,] 7.222224 56.95200 56.95200 56.95200 NaN NaN
< [58,] 7.333335 57.71606 57.71606 57.71606 NaN NaN
< [59,] 7.444446 58.48698 58.48698 58.48698 NaN NaN
< [60,] 7.555558 59.26476 59.26476 59.26476 NaN NaN
< [61,] 7.666669 60.04940 60.04940 60.04940 NaN NaN
< [62,] 7.777780 60.84090 60.84090 60.84090 NaN NaN
< [63,] 7.888892 61.63925 61.63925 61.63925 NaN NaN
< [64,] 8.000003 62.44447 62.44447 62.44447 NaN NaN
< [65,] 8.111114 63.25654 63.25654 63.25654 NaN NaN
< [66,] 8.222226 64.07547 64.07547 64.07547 NaN NaN
< [67,] 8.333337 64.90126 64.90126 64.90126 NaN NaN
< [68,] 8.444448 65.73391 65.73391 65.73391 NaN NaN
< [69,] 8.555560 66.57342 66.57342 66.57342 NaN NaN
< [70,] 8.666671 67.41979 67.41979 67.41979 NaN NaN
< [71,] 8.777782 68.27301 68.27301 68.27301 NaN NaN
< [72,] 8.888894 69.13310 69.13310 69.13310 NaN NaN
< [73,] 9.000005 70.00004 70.00004 70.00004 NaN NaN
< [74,] 9.111116 70.87384 70.87384 70.87384 NaN NaN
< [75,] 9.222228 71.75450 71.75450 71.75450 NaN NaN
< [76,] 9.333339 72.64202 72.64202 72.64202 NaN NaN
< [77,] 9.444450 73.53640 73.53640 73.53640 NaN NaN
< [78,] 9.555562 74.43763 74.43763 74.43763 NaN NaN
< [79,] 9.666673 75.34573 75.34573 75.34573 NaN NaN
< [80,] 9.777784 76.26068 76.26068 76.26068 NaN NaN
< [81,] 9.888896 77.18250 77.18250 77.18250 NaN NaN
< [82,] 10.000007 78.11117 78.11117 78.11117 NaN NaN
< [83,] 10.111118 79.04670 79.04670 79.04670 NaN NaN
< [84,] 10.222230 79.98909 79.98909 79.98909 NaN NaN
< [85,] 10.333341 80.93834 80.93834 80.93834 NaN NaN
< [86,] 10.444452 81.89444 81.89444 81.89444 NaN NaN
< [87,] 10.555564 82.85741 82.85741 82.85741 NaN NaN
< [88,] 10.666675 83.82723 83.82723 83.82723 NaN NaN
< [89,] 10.777786 84.80392 84.80392 84.80392 NaN NaN
< [90,] 10.888898 85.78746 85.78746 85.78746 NaN NaN
< [91,] 11.000009 86.77786 86.77786 86.77786 NaN NaN
< [92,] 11.111120 87.77512 87.77512 87.77512 NaN NaN
< [93,] 11.222232 88.77923 88.77923 88.77923 NaN NaN
< [94,] 11.333343 89.79021 89.79021 89.79021 NaN NaN
< [95,] 11.444454 90.80805 90.80805 90.80805 NaN NaN
< [96,] 11.555566 91.83274 91.83274 91.83274 NaN NaN
< [97,] 11.666677 92.86429 92.86429 92.86429 NaN NaN
< [98,] 11.777788 93.90270 93.90270 93.90270 NaN NaN
< [99,] 11.888900 94.94797 94.94797 94.94797 NaN NaN
< [100,] 12.000011 96.00010 96.00010 96.00010 NaN NaN
< Warning message:
< In qt((1 + level)/2, n - object$k) : NaNs produced
---
> [1,] 0.999989 7.54157 3.84088 11.2423 3.86608 11.2171
> [2,] 1.111100 8.48551 4.93375 12.0373 4.95794 12.0131
> [3,] 1.222212 9.42841 6.01883 12.8380 6.04205 12.8148
> [4,] 1.333323 10.37028 7.09586 13.6447 7.11816 13.6224
> [5,] 1.444434 11.31113 8.16461 14.4576 8.18604 14.4362
> [6,] 1.555546 12.25095 9.22482 15.2771 9.24543 15.2565
> [7,] 1.666657 13.18973 10.27623 16.1032 10.29607 16.0834
> [8,] 1.777768 14.12749 11.31858 16.9364 11.33771 16.9173
> [9,] 1.888880 15.06422 12.35165 17.7768 12.37013 17.7583
> [10,] 1.999991 15.99993 13.37523 18.6246 13.39310 18.6067
> [11,] 2.111102 16.93460 14.38913 19.4801 14.40646 19.4627
> [12,] 2.222214 17.86824 15.39324 20.3432 15.41009 20.3264
> [13,] 2.333325 18.80086 16.38748 21.2142 16.40392 21.1978
> [14,] 2.444436 19.73244 17.37188 22.0930 17.38795 22.0769
> [15,] 2.555548 20.66300 18.34652 22.9795 18.36229 22.9637
> [16,] 2.666659 21.59253 19.31158 23.8735 19.32712 23.8579
> [17,] 2.777770 22.52103 20.26734 24.7747 20.28269 24.7594
> [18,] 2.888882 23.44850 21.21415 25.6828 21.22937 25.6676
> [19,] 2.999993 24.37494 22.15246 26.5974 22.16759 26.5823
> [20,] 3.111104 25.30036 23.08276 27.5179 23.09787 27.5028
> [21,] 3.222216 26.22474 24.00563 28.4439 24.02074 28.4287
> [22,] 3.333327 27.14810 24.92165 29.3745 24.93682 29.3594
> [23,] 3.444438 28.07042 25.83145 30.3094 25.84670 30.2941
> [24,] 3.555550 28.99172 26.73565 31.2478 26.75102 31.2324
> [25,] 3.666661 29.91199 27.63487 32.1891 27.65038 32.1736
> [26,] 3.777772 30.83123 28.52970 33.1328 28.54537 33.1171
> [27,] 3.888884 31.74944 29.42071 34.0782 29.43657 34.0623
> [28,] 3.999995 32.66663 30.30844 35.0248 30.32450 35.0087
> [29,] 4.111106 33.58278 31.19339 35.9722 31.20966 35.9559
> [30,] 4.222218 34.49791 32.07601 36.9198 32.09251 36.9033
> [31,] 4.333329 35.41200 32.95673 37.8673 32.97345 37.8505
> [32,] 4.444440 36.32507 33.83593 38.8142 33.85288 38.7973
> [33,] 4.555552 37.23711 34.71394 39.7603 34.73112 39.7431
> [34,] 4.666663 38.14812 35.59108 40.7052 35.60849 40.6877
> [35,] 4.777774 39.05810 36.46761 41.6486 36.48525 41.6309
> [36,] 4.888886 39.96705 37.34379 42.5903 37.36165 42.5725
> [37,] 4.999997 40.87498 38.21982 43.5301 38.23791 43.5120
> [38,] 5.111108 41.78187 39.09591 44.4678 39.11420 44.4495
> [39,] 5.222220 42.68774 39.97220 45.4033 39.99069 45.3848
> [40,] 5.333331 43.59257 40.84885 46.3363 40.86754 46.3176
> [41,] 5.444442 44.49638 41.72598 47.2668 41.74485 47.2479
> [42,] 5.555554 45.39916 42.60369 48.1946 42.62273 48.1756
> [43,] 5.666665 46.30091 43.48208 49.1197 43.50128 49.1005
> [44,] 5.777776 47.20163 44.36122 50.0421 44.38056 50.0227
> [45,] 5.888888 48.10133 45.24116 50.9615 45.26064 50.9420
> [46,] 5.999999 48.99999 46.12195 51.8780 46.14155 51.8584
> [47,] 6.111110 49.89763 47.00362 52.7916 47.02333 52.7719
> [48,] 6.222222 50.79423 47.88620 53.7023 47.90600 53.6825
> [49,] 6.333333 51.68981 48.76968 54.6099 48.78957 54.5901
> [50,] 6.444444 52.58436 49.65408 55.5146 49.67404 55.4947
> [51,] 6.555556 53.47788 50.53937 56.4164 50.55938 56.3964
> [52,] 6.666667 54.37037 51.42553 57.3152 51.44558 57.2952
> [53,] 6.777778 55.26184 52.31252 58.2112 52.33260 58.1911
> [54,] 6.888890 56.15227 53.20029 59.1043 53.22039 59.0841
> [55,] 7.000001 57.04167 54.08879 59.9946 54.10890 59.9745
> [56,] 7.111112 57.93005 54.97794 60.8822 54.99804 60.8621
> [57,] 7.222224 58.81740 55.86766 61.7671 55.88775 61.7470
> [58,] 7.333335 59.70372 56.75786 62.6496 56.77792 62.6295
> [59,] 7.444446 60.58901 57.64842 63.5296 57.66845 63.5096
> [60,] 7.555558 61.47327 58.53923 64.4073 58.55921 64.3873
> [61,] 7.666669 62.35650 59.43015 65.2829 59.45008 65.2629
> [62,] 7.777780 63.23870 60.32102 66.1564 60.34089 66.1365
> [63,] 7.888892 64.11988 61.21168 67.0281 61.23148 67.0083
> [64,] 8.000003 65.00002 62.10193 67.8981 62.12167 67.8784
> [65,] 8.111114 65.87914 62.99159 68.7667 63.01126 68.7470
> [66,] 8.222226 66.75723 63.88043 69.6340 63.90002 69.6144
> [67,] 8.333337 67.63429 64.76820 70.5004 64.78772 70.4809
> [68,] 8.444448 68.51032 65.65466 71.3660 65.67411 71.3465
> [69,] 8.555560 69.38532 66.53952 72.2311 66.55890 72.2117
> [70,] 8.666671 70.25929 67.42249 73.0961 67.44181 73.0768
> [71,] 8.777782 71.13224 68.30326 73.9612 68.32252 73.9420
> [72,] 8.888894 72.00415 69.18148 74.8268 69.20070 74.8076
> [73,] 9.000005 72.87504 70.05681 75.6933 70.07600 75.6741
> [74,] 9.111116 73.74490 70.92888 76.5609 70.94806 76.5417
> [75,] 9.222228 74.61373 71.79732 77.4301 71.81650 77.4109
> [76,] 9.333339 75.48153 72.66174 78.3013 72.68095 78.2821
> [77,] 9.444450 76.34830 73.52176 79.1748 73.54101 79.1556
> [78,] 9.555562 77.21404 74.37697 80.0511 74.39629 80.0318
> [79,] 9.666673 78.07875 75.22700 80.9305 75.24642 80.9111
> [80,] 9.777784 78.94244 76.07146 81.8134 76.09101 81.7939
> [81,] 9.888896 79.80509 76.90999 82.7002 76.92971 82.6805
> [82,] 10.000007 80.66672 77.74225 83.5912 77.76217 83.5713
> [83,] 10.111118 81.52732 78.56792 84.4867 78.58808 84.4666
> [84,] 10.222230 82.38689 79.38672 85.3871 79.40715 85.3666
> [85,] 10.333341 83.24543 80.19839 86.2925 80.21914 86.2717
> [86,] 10.444452 84.10294 81.00271 87.2032 81.02382 87.1821
> [87,] 10.555564 84.95942 81.79950 88.1194 81.82102 88.0978
> [88,] 10.666675 85.81488 82.58862 89.0411 82.61059 89.0192
> [89,] 10.777786 86.66930 83.36997 89.9686 83.39244 89.9462
> [90,] 10.888898 87.52270 84.14347 90.9019 84.16649 90.8789
> [91,] 11.000009 88.37507 84.90910 91.8410 84.93270 91.8174
> [92,] 11.111120 89.22641 85.66684 92.7860 85.69108 92.7617
> [93,] 11.222232 90.07672 86.41672 93.7367 86.44165 93.7118
> [94,] 11.333343 90.92600 87.15879 94.6932 87.18445 94.6676
> [95,] 11.444454 91.77425 87.89311 95.6554 87.91954 95.6290
> [96,] 11.555566 92.62148 88.61975 96.6232 88.64701 96.5959
> [97,] 11.666677 93.46767 89.33882 97.5965 89.36694 97.5684
> [98,] 11.777788 94.31284 90.05041 98.5753 90.07944 98.5462
> [99,] 11.888900 95.15697 90.75463 99.5593 90.78461 99.5293
> [100,] 12.000011 96.00008 91.45160 100.5486 91.48257 100.5176
Running ‘spline-ex.R’ [2s/3s]
Comparing ‘spline-ex.Rout’ to ‘spline-ex.Rout.save’ ... OK
Running ‘temp.R’ [3s/4s]
Comparing ‘temp.Rout’ to ‘temp.Rout.save’ ... OK
Running ‘wind.R’ [6s/7s]
Running the tests in ‘tests/ex1.R’ failed.
Complete output:
> #### OOps! Running this in 'CMD check' or in *R* __for the first time__
> #### ===== gives a wrong result (at the end) than when run a 2nd time
> ####-- problem disappears with introduction of if (psw) call ... in Fortran
>
> suppressMessages(library(cobs))
> options(digits = 6)
> if(!dev.interactive(orNone=TRUE)) pdf("ex1.pdf")
>
> source(system.file("util.R", package = "cobs"))
>
> ## Simple example from example(cobs)
> set.seed(908)
> x <- seq(-1,1, len = 50)
> f.true <- pnorm(2*x)
> y <- f.true + rnorm(50)/10
> ## specify constraints (boundary conditions)
> con <- rbind(c( 1,min(x),0),
+ c(-1,max(x),1),
+ c( 0, 0, 0.5))
> ## obtain the median *regression* B-spline using automatically selected knots
> coR <- cobs(x,y,constraint = "increase", pointwise = con)
qbsks2():
Performing general knot selection ...
Deleting unnecessary knots ...
> summaryCobs(coR)
List of 24
$ call : language cobs(x = x, y = y, constraint = "increase", pointwise = con)
$ tau : num 0.5
$ degree : num 2
$ constraint : chr "increase"
$ ic : chr "AIC"
$ pointwise : num [1:3, 1:3] 1 -1 0 -1 1 0 0 1 0.5
$ select.knots : logi TRUE
$ select.lambda: logi FALSE
$ x : num [1:50] -1 -0.959 -0.918 -0.878 -0.837 ...
$ y : num [1:50] 0.2254 0.0916 0.0803 -0.0272 -0.0454 ...
$ resid : num [1:50] 0.1976 0.063 0.0491 -0.0626 -0.0868 ...
$ fitted : num [1:50] 0.0278 0.0287 0.0312 0.0354 0.0414 ...
$ coef : num [1:4] 0.0278 0.0278 0.8154 1
$ knots : num [1:3] -1 -0.224 1
$ k0 : num 4
$ k : num 4
$ x.ps :Formal class 'matrix.csr' [package "SparseM"] with 4 slots
$ SSy : num 6.19
$ lambda : num 0
$ icyc : int 7
$ ifl : int 1
$ pp.lambda : NULL
$ pp.sic : NULL
$ i.mask : NULL
cb.lo ci.lo fit ci.up cb.up
1 -6.77514e-02 -0.029701622 0.0278152 0.0853320 0.123382
2 -6.41787e-02 -0.027468888 0.0280224 0.0835138 0.120224
3 -6.04433e-02 -0.024973163 0.0286442 0.0822615 0.117732
4 -5.65412e-02 -0.022212175 0.0296803 0.0815728 0.115902
5 -5.24674e-02 -0.019182756 0.0311310 0.0814447 0.114729
6 -4.82149e-02 -0.015880775 0.0329961 0.0818729 0.114207
7 -4.37751e-02 -0.012301110 0.0352757 0.0828524 0.114326
8 -3.91381e-02 -0.008437641 0.0379697 0.0843771 0.115077
9 -3.42918e-02 -0.004283290 0.0410782 0.0864397 0.116448
10 -2.92233e-02 0.000169901 0.0446012 0.0890325 0.118426
11 -2.39179e-02 0.004930665 0.0485387 0.0921467 0.120995
12 -1.83600e-02 0.010008360 0.0528906 0.0957728 0.124141
13 -1.25335e-02 0.015412811 0.0576570 0.0999012 0.127847
14 -6.42140e-03 0.021154129 0.0628378 0.1045216 0.132097
15 -6.81378e-06 0.027242531 0.0684332 0.1096238 0.136873
16 6.72715e-03 0.033688168 0.0744430 0.1151978 0.142159
17 1.37970e-02 0.040500961 0.0808672 0.1212335 0.147938
18 2.12185e-02 0.047690461 0.0877060 0.1277215 0.154193
19 2.90068e-02 0.055265726 0.0949592 0.1346527 0.160912
20 3.71760e-02 0.063235225 0.1026269 0.1420185 0.168078
21 4.57390e-02 0.071606758 0.1107090 0.1498113 0.175679
22 5.47075e-02 0.080387396 0.1192056 0.1580238 0.183704
23 6.40921e-02 0.089583438 0.1281167 0.1666500 0.192141
24 7.39018e-02 0.099200377 0.1374422 0.1756841 0.200983
25 8.41444e-02 0.109242876 0.1471823 0.1851216 0.210220
26 9.48262e-02 0.119714746 0.1573367 0.1949588 0.219847
27 1.05952e-01 0.130618921 0.1679057 0.2051925 0.229859
28 1.17526e-01 0.141957438 0.1788891 0.2158208 0.240253
29 1.29548e-01 0.153731401 0.1902870 0.2268426 0.251026
30 1.42021e-01 0.165940947 0.2020994 0.2382578 0.262178
31 1.54941e-01 0.178585191 0.2143262 0.2500672 0.273711
32 1.68306e-01 0.191662165 0.2269675 0.2622729 0.285629
33 1.82111e-01 0.205168744 0.2400233 0.2748778 0.297936
34 1.96348e-01 0.219100556 0.2534935 0.2878865 0.310639
35 2.11008e-01 0.233451886 0.2673782 0.3013046 0.323748
36 2.26079e-01 0.248215565 0.2816774 0.3151392 0.337276
37 2.41547e-01 0.263382876 0.2963910 0.3293992 0.351235
38 2.57393e-01 0.278943451 0.3115191 0.3440948 0.365645
39 2.73599e-01 0.294885220 0.3270617 0.3592382 0.380524
40 2.90023e-01 0.311080514 0.3429107 0.3747410 0.395798
41 3.06194e-01 0.327075735 0.3586411 0.3902065 0.411088
42 3.22074e-01 0.342831649 0.3742095 0.4055873 0.426345
43 3.37676e-01 0.358355597 0.3896158 0.4208761 0.441556
44 3.53012e-01 0.373655096 0.4048602 0.4360653 0.456709
45 3.68094e-01 0.388737688 0.4199426 0.4511475 0.471791
46 3.82936e-01 0.403610792 0.4348630 0.4661151 0.486790
47 3.97549e-01 0.418281590 0.4496214 0.4809611 0.501694
48 4.11944e-01 0.432756923 0.4642177 0.4956786 0.516491
49 4.26133e-01 0.447043216 0.4786521 0.5102611 0.531172
50 4.40124e-01 0.461146429 0.4929245 0.5247027 0.545725
51 4.53927e-01 0.475072016 0.5070350 0.5389979 0.560143
52 4.67551e-01 0.488824911 0.5209834 0.5531418 0.574416
53 4.81002e-01 0.502409521 0.5347698 0.5671300 0.588538
54 4.94287e-01 0.515829730 0.5483942 0.5809587 0.602501
55 5.07412e-01 0.529088909 0.5618566 0.5946243 0.616302
56 5.20381e-01 0.542189933 0.5751571 0.6081242 0.629933
57 5.33198e-01 0.555135196 0.5882955 0.6214558 0.643393
58 5.45867e-01 0.567926630 0.6012719 0.6346172 0.656677
59 5.58390e-01 0.580565721 0.6140864 0.6476070 0.669782
60 5.70769e-01 0.593053527 0.6267388 0.6604241 0.682708
61 5.83005e-01 0.605390690 0.6392293 0.6730679 0.695454
62 5.95098e-01 0.617577451 0.6515577 0.6855380 0.708017
63 6.07048e-01 0.629613656 0.6637242 0.6978347 0.720400
64 6.18854e-01 0.641498766 0.6757287 0.7099586 0.732603
65 6.30515e-01 0.653231865 0.6875711 0.7219104 0.744627
66 6.42028e-01 0.664811658 0.6992516 0.7336916 0.756475
67 6.53391e-01 0.676236478 0.7107701 0.7453037 0.768149
68 6.64600e-01 0.687504287 0.7221266 0.7567489 0.779653
69 6.75652e-01 0.698612675 0.7333211 0.7680295 0.790991
70 6.86541e-01 0.709558867 0.7443536 0.7791483 0.802166
71 6.97262e-01 0.720339721 0.7552241 0.7901084 0.813186
72 7.07810e-01 0.730951740 0.7659326 0.8009134 0.824055
73 7.18179e-01 0.741391078 0.7764791 0.8115671 0.834779
74 7.28361e-01 0.751653555 0.7868636 0.8220736 0.845367
75 7.38348e-01 0.761734678 0.7970861 0.8324375 0.855824
76 7.48134e-01 0.771629669 0.8071466 0.8426636 0.866160
77 7.57709e-01 0.781333498 0.8170452 0.8527568 0.876382
78 7.67065e-01 0.790840929 0.8267817 0.8627224 0.886499
79 7.76192e-01 0.800146569 0.8363562 0.8725659 0.896520
80 7.85083e-01 0.809244928 0.8457688 0.8822926 0.906455
81 7.93727e-01 0.818130488 0.8550193 0.8919081 0.916312
82 8.02116e-01 0.826797774 0.8641079 0.9014179 0.926100
83 8.10240e-01 0.835241429 0.8730344 0.9108274 0.935829
84 8.18091e-01 0.843456291 0.8817990 0.9201417 0.945507
85 8.25661e-01 0.851437463 0.8904015 0.9293656 0.955142
86 8.32942e-01 0.859180385 0.8988421 0.9385038 0.964742
87 8.39928e-01 0.866680887 0.9071207 0.9475605 0.974313
88 8.46612e-01 0.873935236 0.9152373 0.9565393 0.983862
89 8.52989e-01 0.880940170 0.9231918 0.9654435 0.993395
90 8.59054e-01 0.887692913 0.9309844 0.9742760 1.002915
91 8.64803e-01 0.894191180 0.9386150 0.9830389 1.012427
92 8.70233e-01 0.900433167 0.9460836 0.9917341 1.021934
93 8.75343e-01 0.906417527 0.9533902 1.0003629 1.031437
94 8.80130e-01 0.912143340 0.9605348 1.0089263 1.040939
95 8.84594e-01 0.917610075 0.9675174 1.0174248 1.050441
96 8.88735e-01 0.922817542 0.9743381 1.0258586 1.059942
97 8.92551e-01 0.927765853 0.9809967 1.0342275 1.069442
98 8.96045e-01 0.932455371 0.9874933 1.0425312 1.078941
99 8.99218e-01 0.936886669 0.9938279 1.0507692 1.088438
100 9.02069e-01 0.941060487 1.0000006 1.0589406 1.097932
knots :
[1] -1.00000 -0.22449 1.00000
coef :
[1] 0.0278152 0.0278152 0.8153868 1.0000006
> coR1 <- cobs(x,y,constraint = "increase", pointwise = con, degree = 1)
qbsks2():
Performing general knot selection ...
Deleting unnecessary knots ...
> summary(coR1)
COBS regression spline (degree = 1) from call:
cobs(x = x, y = y, constraint = "increase", degree = 1, pointwise = con)
{tau=0.5}-quantile; dimensionality of fit: 4 from {4}
x$knots[1:4]: -1.000002, -0.632653, 0.183673, 1.000002
with 3 pointwise constraints
coef[1:4]: 0.0504467, 0.0504467, 0.6305155, 1.0000009
R^2 = 93.83% ; empirical tau (over all): 21/50 = 0.42 (target tau= 0.5)
>
> ## compute the median *smoothing* B-spline using automatically chosen lambda
> coS <- cobs(x,y,constraint = "increase", pointwise = con,
+ lambda = -1, trace = 3)
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.
loo.design2(): -> Xeq 51 x 22 (nz = 151 =^= 0.13%)
Xieq 62 x 22 (nz = 224 =^= 0.16%)
........................
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.
> with(coS, cbind(pp.lambda, pp.sic, k0, ifl, icyc))
pp.lambda pp.sic k0 ifl icyc
[1,] 3.54019e-05 -2.64644 22 1 21
[2,] 6.92936e-05 -2.64644 22 1 21
[3,] 1.35631e-04 -2.64644 22 1 20
[4,] 2.65477e-04 -2.64644 22 1 22
[5,] 5.19629e-04 -2.64644 22 1 22
[6,] 1.01709e-03 -2.64644 22 1 23
[7,] 1.99080e-03 -2.68274 21 1 20
[8,] 3.89667e-03 -2.75212 19 1 18
[9,] 7.62711e-03 -2.73932 19 1 14
[10,] 1.49289e-02 -2.85261 16 1 13
[11,] 2.92209e-02 -2.97873 12 1 12
[12,] 5.71953e-02 -3.01058 11 1 12
[13,] 1.11951e-01 -3.04364 10 1 11
[14,] 2.19126e-01 -3.11242 8 1 12
[15,] 4.28904e-01 -3.17913 6 1 12
[16,] 8.39512e-01 -3.18824 5 1 11
[17,] 1.64321e+00 -3.01467 5 1 12
[18,] 3.21633e+00 -3.01380 4 1 11
[19,] 6.29545e+00 -3.01380 4 1 10
[20,] 1.23223e+01 -3.01380 4 1 11
[21,] 2.41190e+01 -3.01380 4 1 11
[22,] 4.72092e+01 -3.01380 4 1 10
[23,] 9.24046e+01 -3.01380 4 1 10
[24,] 1.80867e+02 -3.01380 4 1 10
[25,] 3.54019e+02 -3.01380 4 1 10
> with(coS, plot(pp.sic ~ pp.lambda, type = "b", log = "x", col=2,
+ main = deparse(call)))
> ##-> very nice minimum close to 1
>
> summaryCobs(coS)
List of 24
$ call : language cobs(x = x, y = y, constraint = "increase", lambda = -1, pointwise = con, trace = 3)
$ tau : num 0.5
$ degree : num 2
$ constraint : chr "increase"
$ ic : NULL
$ pointwise : num [1:3, 1:3] 1 -1 0 -1 1 0 0 1 0.5
$ select.knots : logi TRUE
$ select.lambda: logi TRUE
$ x : num [1:50] -1 -0.959 -0.918 -0.878 -0.837 ...
$ y : num [1:50] 0.2254 0.0916 0.0803 -0.0272 -0.0454 ...
$ resid : num [1:50] 0.2254 0.0829 0.062 -0.0562 -0.0862 ...
$ fitted : num [1:50] 0 0.00869 0.01837 0.02906 0.04075 ...
$ coef : num [1:22] 0 0.00819 0.03365 0.06662 0.10458 ...
$ knots : num [1:20] -1 -0.918 -0.796 -0.714 -0.592 ...
$ k0 : int [1:25] 22 22 22 22 22 22 21 19 19 16 ...
$ k : int 5
$ x.ps :Formal class 'matrix.csr' [package "SparseM"] with 4 slots
$ SSy : num 6.19
$ lambda : Named num 0.84
..- attr(*, "names")= chr "lambda"
$ icyc : int [1:25] 21 21 20 22 22 23 20 18 14 13 ...
$ ifl : int [1:25] 1 1 1 1 1 1 1 1 1 1 ...
$ pp.lambda : num [1:25] 0 0 0 0 0.001 0.001 0.002 0.004 0.008 0.015 ...
$ pp.sic : num [1:25] -2.65 -2.65 -2.65 -2.65 -2.65 ...
$ i.mask : logi [1:25] TRUE TRUE TRUE TRUE TRUE TRUE ...
cb.lo ci.lo fit ci.up cb.up
1 -0.07071332 -0.03907635 -3.77249e-07 0.0390756 0.0707126
2 -0.06555125 -0.03435600 4.17438e-03 0.0427048 0.0739000
3 -0.06016465 -0.02940203 8.59400e-03 0.0465900 0.0773526
4 -0.05455349 -0.02421442 1.32585e-02 0.0507314 0.0810704
5 -0.04871809 -0.01879334 1.81678e-02 0.0551289 0.0850537
6 -0.04265897 -0.01313909 2.33220e-02 0.0597831 0.0893029
7 -0.03637554 -0.00725134 2.87210e-02 0.0646934 0.0938176
8 -0.02986704 -0.00112966 3.43649e-02 0.0698595 0.0985969
9 -0.02313305 0.00522618 4.02537e-02 0.0752812 0.1036404
10 -0.01617351 0.01181620 4.63873e-02 0.0809584 0.1089481
11 -0.00898880 0.01864020 5.27658e-02 0.0868914 0.1145204
12 -0.00157983 0.02569768 5.93891e-02 0.0930806 0.1203581
13 0.00605308 0.03298846 6.62573e-02 0.0995262 0.1264615
14 0.01391000 0.04051257 7.33704e-02 0.1062282 0.1328307
15 0.02199057 0.04826981 8.07283e-02 0.1131867 0.1394660
16 0.03029461 0.05626010 8.83310e-02 0.1204020 0.1463675
17 0.03882336 0.06448412 9.61787e-02 0.1278732 0.1535339
18 0.04757769 0.07294234 1.04271e-01 0.1355999 0.1609646
19 0.05655804 0.08163500 1.12608e-01 0.1435819 0.1686589
20 0.06576441 0.09056212 1.21191e-01 0.1518192 0.1766169
21 0.07519637 0.09972344 1.30018e-01 0.1603120 0.1848391
22 0.08485262 0.10911826 1.39090e-01 0.1690610 0.1933266
23 0.09473211 0.11874598 1.48406e-01 0.1780668 0.2020807
24 0.10483493 0.12860668 1.57968e-01 0.1873294 0.2111011
25 0.11516076 0.13870015 1.67775e-01 0.1968489 0.2203882
26 0.12570956 0.14902638 1.77826e-01 0.2066253 0.2299421
27 0.13648327 0.15958645 1.88122e-01 0.2166576 0.2397608
28 0.14748286 0.17038090 1.98663e-01 0.2269453 0.2498433
29 0.15870881 0.18140998 2.09449e-01 0.2374880 0.2601892
30 0.17016110 0.19267368 2.20480e-01 0.2482859 0.2707984
31 0.18183922 0.20417172 2.31755e-01 0.2593391 0.2816716
32 0.19374227 0.21590361 2.43276e-01 0.2706482 0.2928095
33 0.20587062 0.22786955 2.55041e-01 0.2822129 0.3042118
34 0.21822524 0.24007008 2.67051e-01 0.2940328 0.3158776
35 0.23080666 0.25250549 2.79306e-01 0.3061075 0.3278063
36 0.24361488 0.26517577 2.91806e-01 0.3184370 0.3399979
37 0.25664938 0.27808064 3.04551e-01 0.3310217 0.3524530
38 0.26990862 0.29121926 3.17541e-01 0.3438624 0.3651730
39 0.28339034 0.30459037 3.30775e-01 0.3569602 0.3781603
40 0.29709467 0.31819405 3.44255e-01 0.3703152 0.3914146
41 0.31102144 0.33203019 3.57979e-01 0.3839275 0.4049363
42 0.32517059 0.34609876 3.71948e-01 0.3977971 0.4187252
43 0.33954481 0.36040126 3.86162e-01 0.4119224 0.4327789
44 0.35414537 0.37493839 4.00621e-01 0.4263028 0.4470958
45 0.36897279 0.38971043 4.15324e-01 0.4409381 0.4616757
46 0.38402708 0.40471738 4.30273e-01 0.4558281 0.4765184
47 0.39930767 0.41995895 4.45466e-01 0.4709732 0.4916245
48 0.41479557 0.43541678 4.60887e-01 0.4863568 0.5069780
49 0.43039487 0.45099622 4.76442e-01 0.5018872 0.5224885
50 0.44609197 0.46668362 4.92117e-01 0.5175506 0.5381422
51 0.46188684 0.48247895 5.07913e-01 0.5333471 0.5539392
52 0.47773555 0.49833835 5.23786e-01 0.5492329 0.5698357
53 0.49336687 0.51398935 5.39461e-01 0.5649325 0.5855550
54 0.50873469 0.52938518 5.54891e-01 0.5803975 0.6010480
55 0.52383955 0.54452615 5.70077e-01 0.5956277 0.6163143
56 0.53868141 0.55941225 5.85018e-01 0.6106231 0.6313539
57 0.55325974 0.57404316 5.99714e-01 0.6253839 0.6461673
58 0.56757320 0.58841816 6.14165e-01 0.6399109 0.6607558
59 0.58161907 0.60253574 6.28371e-01 0.6542056 0.6751223
60 0.59539741 0.61639593 6.42332e-01 0.6682680 0.6892665
61 0.60890835 0.62999881 6.56048e-01 0.6820980 0.7031884
62 0.62215175 0.64334429 6.69520e-01 0.6956957 0.7168882
63 0.63512996 0.65643368 6.82747e-01 0.7090597 0.7303634
64 0.64784450 0.66926783 6.95729e-01 0.7221893 0.7436126
65 0.66029589 0.68184700 7.08466e-01 0.7350841 0.7566352
66 0.67248408 0.69417118 7.20958e-01 0.7477442 0.7694313
67 0.68440855 0.70624008 7.33205e-01 0.7601699 0.7820014
68 0.69606829 0.71805313 7.45207e-01 0.7723617 0.7943465
69 0.70746295 0.72961016 7.56965e-01 0.7843198 0.8064670
70 0.71859343 0.74091165 7.68478e-01 0.7960438 0.8183620
71 0.72946023 0.75195789 7.79746e-01 0.8075332 0.8300309
72 0.74006337 0.76274887 7.90769e-01 0.8187883 0.8414738
73 0.75040233 0.77328433 8.01547e-01 0.8298091 0.8526911
74 0.76047612 0.78356369 8.12080e-01 0.8405963 0.8636839
75 0.77028266 0.79358583 8.22368e-01 0.8511510 0.8744542
76 0.77982200 0.80335076 8.32412e-01 0.8614732 0.8850020
77 0.78909446 0.81285866 8.42211e-01 0.8715627 0.8953269
78 0.79809990 0.82210946 8.51765e-01 0.8814196 0.9054292
79 0.80683951 0.83110382 8.61074e-01 0.8910433 0.9153076
80 0.81531459 0.83984244 8.70138e-01 0.9004329 0.9249608
81 0.82352559 0.84832559 8.78957e-01 0.9095884 0.9343884
82 0.83147249 0.85655324 8.87531e-01 0.9185095 0.9435903
83 0.83915483 0.86452515 8.95861e-01 0.9271968 0.9525671
84 0.84657171 0.87224082 9.03946e-01 0.9356505 0.9613196
85 0.85372180 0.87969951 9.11786e-01 0.9438715 0.9698492
86 0.86060525 0.88690131 9.19381e-01 0.9518597 0.9781558
87 0.86722242 0.89384640 9.26731e-01 0.9596149 0.9862389
88 0.87357322 0.90053476 9.33836e-01 0.9671371 0.9940986
89 0.87965804 0.90696658 9.40696e-01 0.9744261 1.0017347
90 0.88547781 0.91314239 9.47312e-01 0.9814814 1.0091460
91 0.89103290 0.91906239 9.53683e-01 0.9883028 1.0163323
92 0.89632328 0.92472655 9.59808e-01 0.9948904 1.0232937
93 0.90134850 0.93013464 9.65689e-01 1.0012443 1.0300304
94 0.90610776 0.93528622 9.71326e-01 1.0073650 1.0365434
95 0.91060065 0.94018104 9.76717e-01 1.0132527 1.0428331
96 0.91482784 0.94481950 9.81863e-01 1.0189071 1.0488987
97 0.91878971 0.94920179 9.86765e-01 1.0243279 1.0547400
98 0.92248624 0.95332789 9.91422e-01 1.0295152 1.0603569
99 0.92591703 0.95719761 9.95833e-01 1.0344692 1.0657498
100 0.92908136 0.96081053 1.00000e+00 1.0391902 1.0709194
knots :
[1] -1.0000020 -0.9183673 -0.7959184 -0.7142857 -0.5918367 -0.5102041
[7] -0.3877551 -0.2653061 -0.1836735 -0.0612245 0.0204082 0.1428571
[13] 0.2244898 0.3469388 0.4693878 0.5510204 0.6734694 0.7551020
[19] 0.8775510 1.0000020
coef :
[1] -4.01161e-07 8.18714e-03 3.36534e-02 6.66159e-02 1.04576e-01
[6] 1.50032e-01 2.00486e-01 2.70027e-01 3.35473e-01 4.05918e-01
[11] 4.83858e-01 5.64259e-01 6.37163e-01 7.05069e-01 7.77561e-01
[16] 8.30474e-01 8.78390e-01 9.18810e-01 9.54232e-01 9.87743e-01
[21] 1.00000e+00 5.99960e-01
>
> plot(x, y, main = "cobs(x,y, constraint=\"increase\", pointwise = *)")
> matlines(x, cbind(fitted(coR), fitted(coR1), fitted(coS)),
+ col = 2:4, lty=1)
>
> ##-- real data example (still n = 50)
> data(cars)
> attach(cars)
> co1 <- cobs(speed, dist, "increase")
qbsks2():
Performing general knot selection ...
Deleting unnecessary knots ...
> co1.1 <- cobs(speed, dist, "increase", knots.add = TRUE)
qbsks2():
Performing general knot selection ...
Deleting unnecessary knots ...
Searching for missing knots ...
> co1.2 <- cobs(speed, dist, "increase", knots.add = TRUE, repeat.delete.add = TRUE)
qbsks2():
Performing general knot selection ...
Deleting unnecessary knots ...
Searching for missing knots ...
> ## These three all give the same -- only remaining knots (outermost data):
> ic <- which("call" == names(co1))
> stopifnot(all.equal(co1[-ic], co1.1[-ic]),
+ all.equal(co1[-ic], co1.2[-ic]))
> 1 - sum(co1 $ resid ^2) / sum((dist - mean(dist))^2) # R^2 = 64.2%
[1] 0.642288
>
> co2 <- cobs(speed, dist, "increase", lambda = -1)# 6 warnings
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
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.10.11-amd64"
version
"#1 SMP PREEMPT_DYNAMIC Debian 6.10.11-1 (2024-09-22)"
nodename
"gimli2"
machine
"x86_64"
login
"hornik"
user
"hornik"
effective_user
"hornik"
> 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%)
Warning message:
In cobs(x, y, trace = 3) :
drqssbc2(): Not all flags are normal (== 1), ifl : 21
> 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 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 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%)
Warning message:
In cobs(x, y, constraint = c("increase"), trace = 3) :
drqssbc2(): Not all flags are normal (== 1), ifl : 21
> 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 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 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%)
Warning message:
In cobs(x, y, constraint = c("concave"), trace = 3) :
drqssbc2(): Not all flags are normal (== 1), ifl : 21
>
> 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 4 (nz = 150 =^= 0.75%)
Xieq 5 x 4 (nz = 15 =^= 0.75%)
loo.design2(): -> Xeq 50 x 4 (nz = 150 =^= 0.75%)
Xieq 5 x 4 (nz = 15 =^= 0.75%)
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%)
Warning message:
In cobs(x, y, constraint = c("inc", "concave"), trace = 3) :
drqssbc2(): Not all flags are normal (== 1), ifl : 20
> ## here, it's the same as just "i":
> all.equal(fitted(c2i), fitted(c2IC))
[1] TRUE
>
> 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%)
Warning message:
In cobs(x, y, degree = 1, trace = 3) :
drqssbc2(): Not all flags are normal (== 1), ifl : 22
> 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%)
Warning message:
In cobs(x, y, degree = 1, constraint = c("increase"), trace = 3) :
drqssbc2(): Not all flags are normal (== 1), ifl : 22
> 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 4 (nz = 100 =^= 0.5%)
Xieq 2 x 4 (nz = 6 =^= 0.75%)
l1.design2(): -> Xeq 50 x 4 (nz = 100 =^= 0.5%)
Xieq 2 x 4 (nz = 6 =^= 0.75%)
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%)
Warning message:
In cobs(x, y, degree = 1, constraint = c("concave"), trace = 3) :
drqssbc2(): Not all flags are normal (== 1), ifl : 21
>
> 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 4 (nz = 100 =^= 0.5%)
Xieq 5 x 4 (nz = 12 =^= 0.6%)
l1.design2(): -> Xeq 50 x 4 (nz = 100 =^= 0.5%)
Xieq 5 x 4 (nz = 12 =^= 0.6%)
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%)
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 cobs(x, y, degree = 1, constraint = c("inc", "concave"), trace = 3) :
drqssbc2(): Not all flags are normal (== 1), ifl : 21
>
> 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 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%)
Warning message:
In cobs(x, y, pointwise = con, trace = 3) :
drqssbc2(): Not all flags are normal (== 1), ifl : 21
>
> ## 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 4 (nz = 150 =^= 0.75%)
Xieq 4 x 4 (nz = 12 =^= 0.75%)
loo.design2(): -> Xeq 50 x 4 (nz = 150 =^= 0.75%)
Xieq 4 x 4 (nz = 12 =^= 0.75%)
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%)
Warning message:
In cobs(x, y, constraint = "concave", pointwise = con, trace = 3) :
drqssbc2(): Not all flags are normal (== 1), ifl : 21
>
> ## 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 4 (nz = 150 =^= 0.75%)
Xieq 7 x 4 (nz = 21 =^= 0.75%)
loo.design2(): -> Xeq 50 x 4 (nz = 150 =^= 0.75%)
Xieq 7 x 4 (nz = 21 =^= 0.75%)
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%)
Warning message:
In cobs(x, y, constraint = c("inc", "concave"), pointwise = con, :
drqssbc2(): Not all flags are normal (== 1), ifl : 20
> ## 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 : 22
> 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%)
Warning message:
In cobs(x, y, degree = 1, constraint = "increase", pointwise = con, :
drqssbc2(): Not all flags are normal (== 1), ifl : 22
> 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 4 (nz = 100 =^= 0.5%)
Xieq 4 x 4 (nz = 10 =^= 0.62%)
l1.design2(): -> Xeq 50 x 4 (nz = 100 =^= 0.5%)
Xieq 4 x 4 (nz = 10 =^= 0.62%)
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%)
Warning message:
In cobs(x, y, degree = 1, constraint = "concave", pointwise = con, :
drqssbc2(): Not all flags are normal (== 1), ifl : 21
>
> 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 4 (nz = 100 =^= 0.5%)
Xieq 7 x 4 (nz = 16 =^= 0.57%)
l1.design2(): -> Xeq 50 x 4 (nz = 100 =^= 0.5%)
Xieq 7 x 4 (nz = 16 =^= 0.57%)
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%)
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 cobs(x, y, degree = 1, constraint = c("inc", "concave"), pointwise = con, :
drqssbc2(): Not all flags are normal (== 1), ifl : 21
>
> ## 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
$ c2 : num [1:3] -1 -0.449 2
$ c2IC : num [1:3] -1 -0.449 2
$ c2c : num [1:3] -1 -0.449 2
$ c2i : num [1:3] -1 -0.449 2
$ cp2 : num [1:3] -1 -0.449 2
$ cp2IC: num [1:3] -1 -0.449 2
$ cp2c : num [1:3] -1 -0.449 2
$ cp2i : num [1:3] -1 -0.449 2
$ c1 : num [1:5] -1 -0.449 0.163 0.776 2
$ c1IC : num [1:5] -1 -0.449 0.163 0.776 2
$ c1c : 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
$ cp1IC: num [1:5] -1 -0.449 0.163 0.776 2
$ cp1c : 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.95341697, c1IC = 0.95341697, c1c = 0.95341697, c1i = 0.95341697,
c2 = 0.94864721, c2IC = 0.94864721, c2c = 0.94864721, c2i = 0.94864721,
cp1 = 0.95341697, cp1IC = 0.95341697, cp1c = 0.95341697, cp1i = 0.95341697,
cp2 = 0.94864721, cp2IC = 0.94864721, cp2c = 0.94864721, cp2i = 0.94864721
)
> all.equal(Rsqrs[iR2], gotRsqrs[iR2], tolerance=0)# 2.6277e-9 (Lnx F 38); 2.6898e-9 (M1 mac)
[1] "Mean relative difference: 0.0179511"
> all.equal(Rsqrs[nR2], gotRsqrs[nR2], tolerance=0)# differ; drastically only for 'noLD'
[1] "Mean relative difference: 0.0330278"
> ## 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.0179511
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.001 0.002
>
> 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.415 0.027 0.526
>
> 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.032 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.714 0.003 2.145
>
> 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
Flavor: r-devel-linux-x86_64-debian-clang
Version: 1.3-8
Check: tests
Result: ERROR
Running ‘0_pt-ex.R’ [2s/2s]
Running ‘ex1.R’ [5s/6s]
Running ‘ex2-long.R’ [5s/6s]
Running ‘ex3.R’ [2s/2s]
Comparing ‘ex3.Rout’ to ‘ex3.Rout.save’ ...15,16c15,16
< Warning messages:
< 1: In cobs(weight, height, knots = weight, nknots = length(weight)) :
---
> Warning message:
> In cobs(weight, height, knots = weight, nknots = length(weight)) :
19,20d18
< 2: In cobs(weight, height, knots = weight, nknots = length(weight)) :
< drqssbc2(): Not all flags are normal (== 1), ifl : 18
Running ‘multi-constr.R’ [3s/4s]
Running ‘roof.R’ [3s/3s]
Comparing ‘roof.Rout’ to ‘roof.Rout.save’ ... OK
Running ‘small-ex.R’ [2s/2s]
Comparing ‘small-ex.Rout’ to ‘small-ex.Rout.save’ ...24,25d23
< Warning message:
< In cobs(x, y) : drqssbc2(): Not all flags are normal (== 1), ifl : 18
Running ‘spline-ex.R’ [2s/2s]
Comparing ‘spline-ex.Rout’ to ‘spline-ex.Rout.save’ ... OK
Running ‘temp.R’ [2s/2s]
Comparing ‘temp.Rout’ to ‘temp.Rout.save’ ...445,446d444
< Warning in cobs(year, temp, knots.add = TRUE, degree = 1, constraint = "none", :
< The algorithm has not converged after 100 iterations
456,457d453
< Warning in cobs(year, temp, nknots = 9, knots.add = TRUE, degree = 1, constraint = "none", :
< The algorithm has not converged after 100 iterations
463,464c459,460
< coef[1:5]: -0.39324840, -0.28115087, 0.05916295, -0.07465159, 0.31227753
< R^2 = 73.22% ; empirical tau (over all): 63/113 = 0.5575221 (target tau= 0.5)
---
> coef[1:5]: -0.40655906, -0.31473700, 0.05651823, -0.05681818, 0.28681956
> R^2 = 72.56% ; empirical tau (over all): 54/113 = 0.4778761 (target tau= 0.5)
470,473d465
<
< **** ERROR in algorithm: ifl = 18
<
<
476,478d467
< Warning message:
< In cobs(year, temp, nknots = length(a50$knots), knots = a50$knot, :
< drqssbc2(): Not all flags are normal (== 1), ifl : 18
481,482d469
< Warning in cobs(year, temp, nknots = length(a50$knots), knots = a50$knot, :
< The algorithm has not converged after 100 iterations
488,490c475
< [1] 1 2 9 10 17 18 20 21 22 23 26 27 35 36 42 47 48 49 52
< [20] 53 58 59 61 62 63 64 65 68 73 74 78 79 80 81 82 83 84 88
< [39] 90 91 94 98 100 101 102 104 108 109 111 112
---
> [1] 10 18 21 22 47 61 68 74 78 79 102 111
492,495c477
< [1] 3 4 5 6 7 8 11 12 13 14 15 16 19 24 25 28 29 30 31
< [20] 32 33 34 37 38 39 40 41 43 44 45 46 50 51 54 55 56 57 60
< [39] 66 67 69 70 71 72 75 76 77 85 86 87 89 92 93 95 96 97 99
< [58] 103 105 106 107 110 113
---
> [1] 5 8 25 38 39 50 54 77 85 97 113
Running ‘wind.R’ [6s/7s]
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.9.7-amd64"
version
"#1 SMP PREEMPT_DYNAMIC Debian 6.9.7-1 (2024-06-27)"
nodename
"gimli1"
machine
"x86_64"
login
"hornik"
user
"hornik"
effective_user
"hornik"
> 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%)
Warning message:
In cobs(x, y, trace = 3) :
drqssbc2(): Not all flags are normal (== 1), ifl : 18
> 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%)
Warning in cobs(x, y, degree = 1, trace = 3) :
The algorithm has not converged after 100 iterations
> 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] "Mean relative difference: 0.0215671"
>
> ## 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
NULL
> ## 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%)
>
> ## 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%)
> 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.95341697, c1IC = 0.92974549, c1c = 0.92974549, c1i = 0.95079126,
c2 = 0.94864721, 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
)
> all.equal(Rsqrs[iR2], gotRsqrs[iR2], tolerance=0)# 2.6277e-9 (Lnx F 38); 2.6898e-9 (M1 mac)
[1] "Mean relative difference: 0.000374227"
> all.equal(Rsqrs[nR2], gotRsqrs[nR2], tolerance=0)# differ; drastically only for 'noLD'
[1] "Mean relative difference: 2.18514e-09"
> ## 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.000374227
Execution halted
Flavor: r-devel-linux-x86_64-debian-gcc
Version: 1.3-8
Check: tests
Result: ERROR
Running ‘0_pt-ex.R’
Running ‘ex1.R’ [13s/13s]
Running ‘ex2-long.R’ [12s/12s]
Running ‘ex3.R’
Comparing ‘ex3.Rout’ to ‘ex3.Rout.save’ ...15,16c15,16
< Warning messages:
< 1: In cobs(weight, height, knots = weight, nknots = length(weight)) :
---
> Warning message:
> In cobs(weight, height, knots = weight, nknots = length(weight)) :
19,20d18
< 2: In cobs(weight, height, knots = weight, nknots = length(weight)) :
< drqssbc2(): Not all flags are normal (== 1), ifl : 18
Running ‘multi-constr.R’
Running ‘roof.R’
Comparing ‘roof.Rout’ to ‘roof.Rout.save’ ... OK
Running ‘small-ex.R’
Comparing ‘small-ex.Rout’ to ‘small-ex.Rout.save’ ...24,25d23
< Warning message:
< In cobs(x, y) : drqssbc2(): Not all flags are normal (== 1), ifl : 18
Running ‘spline-ex.R’
Comparing ‘spline-ex.Rout’ to ‘spline-ex.Rout.save’ ... OK
Running ‘temp.R’
Comparing ‘temp.Rout’ to ‘temp.Rout.save’ ...445,446d444
< Warning in cobs(year, temp, knots.add = TRUE, degree = 1, constraint = "none", :
< The algorithm has not converged after 100 iterations
456,457d453
< Warning in cobs(year, temp, nknots = 9, knots.add = TRUE, degree = 1, constraint = "none", :
< The algorithm has not converged after 100 iterations
463,464c459,460
< coef[1:5]: -0.39324840, -0.28115087, 0.05916295, -0.07465159, 0.31227753
< R^2 = 73.22% ; empirical tau (over all): 63/113 = 0.5575221 (target tau= 0.5)
---
> coef[1:5]: -0.40655906, -0.31473700, 0.05651823, -0.05681818, 0.28681956
> R^2 = 72.56% ; empirical tau (over all): 54/113 = 0.4778761 (target tau= 0.5)
470,473d465
<
< **** ERROR in algorithm: ifl = 18
<
<
476,478d467
< Warning message:
< In cobs(year, temp, nknots = length(a50$knots), knots = a50$knot, :
< drqssbc2(): Not all flags are normal (== 1), ifl : 18
481,482d469
< Warning in cobs(year, temp, nknots = length(a50$knots), knots = a50$knot, :
< The algorithm has not converged after 100 iterations
488,490c475
< [1] 1 2 9 10 17 18 20 21 22 23 26 27 35 36 42 47 48 49 52
< [20] 53 58 59 61 62 63 64 65 68 73 74 78 79 80 81 82 83 84 88
< [39] 90 91 94 98 100 101 102 104 108 109 111 112
---
> [1] 10 18 21 22 47 61 68 74 78 79 102 111
492,495c477
< [1] 3 4 5 6 7 8 11 12 13 14 15 16 19 24 25 28 29 30 31
< [20] 32 33 34 37 38 39 40 41 43 44 45 46 50 51 54 55 56 57 60
< [39] 66 67 69 70 71 72 75 76 77 85 86 87 89 92 93 95 96 97 99
< [58] 103 105 106 107 110 113
---
> [1] 5 8 25 38 39 50 54 77 85 97 113
Running ‘wind.R’ [14s/15s]
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%)
Warning message:
In cobs(x, y, trace = 3) :
drqssbc2(): Not all flags are normal (== 1), ifl : 18
> 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%)
Warning in cobs(x, y, degree = 1, trace = 3) :
The algorithm has not converged after 100 iterations
> 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] "Mean relative difference: 0.0215671"
>
> ## 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
NULL
> ## 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%)
>
> ## 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%)
> 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.95341697, c1IC = 0.92974549, c1c = 0.92974549, c1i = 0.95079126,
c2 = 0.94864721, 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
)
> all.equal(Rsqrs[iR2], gotRsqrs[iR2], tolerance=0)# 2.6277e-9 (Lnx F 38); 2.6898e-9 (M1 mac)
[1] "Mean relative difference: 0.000374227"
> all.equal(Rsqrs[nR2], gotRsqrs[nR2], tolerance=0)# differ; drastically only for 'noLD'
[1] "Mean relative difference: 2.18514e-09"
> ## 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.000374227
Execution halted
Flavor: r-devel-linux-x86_64-fedora-clang
Version: 1.3-8
Check: tests
Result: ERROR
Running ‘0_pt-ex.R’
Running ‘ex1.R’ [12s/40s]
Running ‘ex2-long.R’ [12s/30s]
Running ‘ex3.R’ [4s/14s]
Comparing ‘ex3.Rout’ to ‘ex3.Rout.save’ ...15,16c15,16
< Warning messages:
< 1: In cobs(weight, height, knots = weight, nknots = length(weight)) :
---
> Warning message:
> In cobs(weight, height, knots = weight, nknots = length(weight)) :
19,20d18
< 2: In cobs(weight, height, knots = weight, nknots = length(weight)) :
< drqssbc2(): Not all flags are normal (== 1), ifl : 18
Running ‘multi-constr.R’ [7s/16s]
Running ‘roof.R’ [6s/14s]
Comparing ‘roof.Rout’ to ‘roof.Rout.save’ ... OK
Running ‘small-ex.R’
Comparing ‘small-ex.Rout’ to ‘small-ex.Rout.save’ ...24,25d23
< Warning message:
< In cobs(x, y) : drqssbc2(): Not all flags are normal (== 1), ifl : 18
Running ‘spline-ex.R’
Comparing ‘spline-ex.Rout’ to ‘spline-ex.Rout.save’ ... OK
Running ‘temp.R’ [5s/20s]
Comparing ‘temp.Rout’ to ‘temp.Rout.save’ ...445,446d444
< Warning in cobs(year, temp, knots.add = TRUE, degree = 1, constraint = "none", :
< The algorithm has not converged after 100 iterations
456,457d453
< Warning in cobs(year, temp, nknots = 9, knots.add = TRUE, degree = 1, constraint = "none", :
< The algorithm has not converged after 100 iterations
463,464c459,460
< coef[1:5]: -0.39324840, -0.28115087, 0.05916295, -0.07465159, 0.31227753
< R^2 = 73.22% ; empirical tau (over all): 63/113 = 0.5575221 (target tau= 0.5)
---
> coef[1:5]: -0.40655906, -0.31473700, 0.05651823, -0.05681818, 0.28681956
> R^2 = 72.56% ; empirical tau (over all): 54/113 = 0.4778761 (target tau= 0.5)
470,473d465
<
< **** ERROR in algorithm: ifl = 18
<
<
476,478d467
< Warning message:
< In cobs(year, temp, nknots = length(a50$knots), knots = a50$knot, :
< drqssbc2(): Not all flags are normal (== 1), ifl : 18
481,482d469
< Warning in cobs(year, temp, nknots = length(a50$knots), knots = a50$knot, :
< The algorithm has not converged after 100 iterations
488,490c475
< [1] 1 2 9 10 17 18 20 21 22 23 26 27 35 36 42 47 48 49 52
< [20] 53 58 59 61 62 63 64 65 68 73 74 78 79 80 81 82 83 84 88
< [39] 90 91 94 98 100 101 102 104 108 109 111 112
---
> [1] 10 18 21 22 47 61 68 74 78 79 102 111
492,495c477
< [1] 3 4 5 6 7 8 11 12 13 14 15 16 19 24 25 28 29 30 31
< [20] 32 33 34 37 38 39 40 41 43 44 45 46 50 51 54 55 56 57 60
< [39] 66 67 69 70 71 72 75 76 77 85 86 87 89 92 93 95 96 97 99
< [58] 103 105 106 107 110 113
---
> [1] 5 8 25 38 39 50 54 77 85 97 113
Running ‘wind.R’ [14s/51s]
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%)
Warning message:
In cobs(x, y, trace = 3) :
drqssbc2(): Not all flags are normal (== 1), ifl : 18
> 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%)
Warning in cobs(x, y, degree = 1, trace = 3) :
The algorithm has not converged after 100 iterations
> 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] "Mean relative difference: 0.0215671"
>
> ## 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
NULL
> ## 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%)
>
> ## 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%)
> 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.95341697, c1IC = 0.92974549, c1c = 0.92974549, c1i = 0.95079126,
c2 = 0.94864721, 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
)
> all.equal(Rsqrs[iR2], gotRsqrs[iR2], tolerance=0)# 2.6277e-9 (Lnx F 38); 2.6898e-9 (M1 mac)
[1] "Mean relative difference: 0.000374227"
> all.equal(Rsqrs[nR2], gotRsqrs[nR2], tolerance=0)# differ; drastically only for 'noLD'
[1] "Mean relative difference: 2.18514e-09"
> ## 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.000374227
Execution halted
Flavor: r-devel-linux-x86_64-fedora-gcc
Version: 1.3-8
Check: tests
Result: ERROR
Running '0_pt-ex.R' [2s]
Running 'ex1.R' [6s]
Running 'ex2-long.R' [6s]
Running 'ex3.R' [2s]
Comparing 'ex3.Rout' to 'ex3.Rout.save' ...15,16c15,16
< Warning messages:
< 1: In cobs(weight, height, knots = weight, nknots = length(weight)) :
---
> Warning message:
> In cobs(weight, height, knots = weight, nknots = length(weight)) :
19,20d18
< 2: In cobs(weight, height, knots = weight, nknots = length(weight)) :
< drqssbc2(): Not all flags are normal (== 1), ifl : 18
Running 'multi-constr.R' [3s]
Running 'roof.R' [3s]
Comparing 'roof.Rout' to 'roof.Rout.save' ... OK
Running 'small-ex.R' [2s]
Comparing 'small-ex.Rout' to 'small-ex.Rout.save' ...24,25d23
< Warning message:
< In cobs(x, y) : drqssbc2(): Not all flags are normal (== 1), ifl : 18
Running 'spline-ex.R' [2s]
Comparing 'spline-ex.Rout' to 'spline-ex.Rout.save' ... OK
Running 'temp.R' [2s]
Comparing 'temp.Rout' to 'temp.Rout.save' ...445,446d444
< Warning in cobs(year, temp, knots.add = TRUE, degree = 1, constraint = "none", :
< The algorithm has not converged after 100 iterations
456,457d453
< Warning in cobs(year, temp, nknots = 9, knots.add = TRUE, degree = 1, constraint = "none", :
< The algorithm has not converged after 100 iterations
463,464c459,460
< coef[1:5]: -0.39324840, -0.28115087, 0.05916295, -0.07465159, 0.31227753
< R^2 = 73.22% ; empirical tau (over all): 63/113 = 0.5575221 (target tau= 0.5)
---
> coef[1:5]: -0.40655906, -0.31473700, 0.05651823, -0.05681818, 0.28681956
> R^2 = 72.56% ; empirical tau (over all): 54/113 = 0.4778761 (target tau= 0.5)
470,473d465
<
< **** ERROR in algorithm: ifl = 18
<
<
476,478d467
< Warning message:
< In cobs(year, temp, nknots = length(a50$knots), knots = a50$knot, :
< drqssbc2(): Not all flags are normal (== 1), ifl : 18
481,482d469
< Warning in cobs(year, temp, nknots = length(a50$knots), knots = a50$knot, :
< The algorithm has not converged after 100 iterations
488,490c475
< [1] 1 2 9 10 17 18 20 21 22 23 26 27 35 36 42 47 48 49 52
< [20] 53 58 59 61 62 63 64 65 68 73 74 78 79 80 81 82 83 84 88
< [39] 90 91 94 98 100 101 102 104 108 109 111 112
---
> [1] 10 18 21 22 47 61 68 74 78 79 102 111
492,495c477
< [1] 3 4 5 6 7 8 11 12 13 14 15 16 19 24 25 28 29 30 31
< [20] 32 33 34 37 38 39 40 41 43 44 45 46 50 51 54 55 56 57 60
< [39] 66 67 69 70 71 72 75 76 77 85 86 87 89 92 93 95 96 97 99
< [58] 103 105 106 107 110 113
---
> [1] 5 8 25 38 39 50 54 77 85 97 113
Running 'wind.R' [8s]
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 release version nodename machine
"Windows" "Server x64" "build 20348" "CRANWIN3" "x86-64"
login user effective_user udomain
"CRAN" "CRAN" "CRAN" "CRANWIN3"
> 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%)
Warning message:
In cobs(x, y, trace = 3) :
drqssbc2(): Not all flags are normal (== 1), ifl : 18
> 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%)
Warning in cobs(x, y, degree = 1, trace = 3) :
The algorithm has not converged after 100 iterations
> 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] "Mean relative difference: 0.0215671"
>
> ## 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
NULL
> ## 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%)
>
> ## 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%)
> 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.95341697, c1IC = 0.92974549, c1c = 0.92974549, c1i = 0.95079126,
c2 = 0.94864721, 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
)
> all.equal(Rsqrs[iR2], gotRsqrs[iR2], tolerance=0)# 2.6277e-9 (Lnx F 38); 2.6898e-9 (M1 mac)
[1] "Mean relative difference: 0.000374227"
> all.equal(Rsqrs[nR2], gotRsqrs[nR2], tolerance=0)# differ; drastically only for 'noLD'
[1] "Mean relative difference: 2.18514e-09"
> ## 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.000374227
Execution halted
Flavor: r-devel-windows-x86_64
These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.