R Under development (unstable) (2023-06-09 r84528 ucrt) -- "Unsuffered Consequences" Copyright (C) 2023 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > library(robustbase) > > > ### Back-compatibility check of consistency & finite sample correction for Qn() : > > set.seed(153) > x <- sort(c(rnorm(80), rt(20, df = 1))) > > ix <- c(27, 57, 15, 1, 26, 13, 23, 70, 9, 54, 6, 12, 8, 80, 11, 69, + 41, 53, 10, 37, 50, 21, 48, 51, 71, 17, 30, 16, 45, 28, 55, 5, + 59, 77, 61, 40, 63, 42, 72, 66) > QnA. <- c(0, 0.72307295, 1.2926498, 1.596857, 1.0979815, 0.84209457, + 1.0719335, 0.88620416, 1.0905118, 0.99056842, 1.2229216, 1.0626517, + 1.1738174, 1.1433873, 1.2071829, 1.1562513, 1.2182886, 1.1587793, + 1.1585524, 1.1555462, 1.1376428, 1.0532134, 1.0447343, 1.0200998, + 1.0495224, 1.0120569, 1.0094172, 0.9749928, 0.9530458, 0.92767184, + 0.90922667, 0.98987601, 0.98223857, 1.0053697, 0.98792848, 0.951908, + 0.92226488, 0.92312857, 0.92313406, 0.92733413) > QnAx <- sapply(seq_along(ix), function(n) Qn(x[ix[1:n]])) > stopifnot( all.equal(QnA., QnAx) ) > > > ### -------------- Plots ----------------------------------------------- > > if(!dev.interactive(orNone=TRUE)) pdf("Qn-Sn-plots.pdf") > > plot(QnA., type="o", main="Qn()") > abline(h = 1, lty=2) > > n <- 1:50 > (qnn <- sapply(n, function(n)Qn(1:n, const=1))) [1] 0 1 1 1 1 2 1 2 2 2 2 2 2 3 3 3 3 3 3 4 3 4 4 4 4 4 4 5 4 5 5 5 5 5 5 6 6 6 [39] 6 6 6 6 6 7 7 7 7 7 7 8 > plot(n, qnn, type = 'b', col = 2, + ylab = "Qn", main = "Qn(1:n) [unscaled]") > > (snn <- sapply(n, function(n)Sn(1:n, const=1))) [1] 0 1 1 1 1 2 2 2 2 3 3 3 3 4 4 4 4 5 5 5 5 6 6 6 6 [26] 7 7 7 7 8 8 8 8 9 9 9 9 10 10 10 10 11 11 11 11 12 12 12 12 13 > plot(n, snn, type = 'b', col = 2, + ylab = "Sn", main = "Sn(1:n) [unscaled]") > > matplot(n, cbind(qnn, snn), type = 'b', + ylab = "Qn & Sn", main = "Qn(1:n) & Sn(1:n) [unscaled]") > legend("topleft", c("Qn", "Sn"), col=1:2, lty=1:2, bty="n", pch=paste(1:2)) > > (sdn <- c(1, sapply(n[-1], function(n)sd(1:n)/n))) [1] 1.0000000 0.3535534 0.3333333 0.3227486 0.3162278 0.3118048 0.3086067 [8] 0.3061862 0.3042903 0.3027650 0.3015113 0.3004626 0.2995723 0.2988072 [15] 0.2981424 0.2975595 0.2970443 0.2965855 0.2961744 0.2958040 0.2954684 [22] 0.2951630 0.2948839 0.2946278 0.2943920 0.2941742 0.2939724 0.2937848 [29] 0.2936101 0.2934469 0.2932942 0.2931510 0.2930164 0.2928896 0.2927700 [36] 0.2926570 0.2925501 0.2924488 0.2923527 0.2922613 0.2921744 0.2920915 [43] 0.2920125 0.2919371 0.2918650 0.2917960 0.2917300 0.2916667 0.2916059 [50] 0.2915476 > ## sd(1) => NA > > for(Sample in c(function(n) ppoints(n), + function(n) qnorm(ppoints(n)))) + { + ##mult.fig(2) : + op <- par(mfrow=c(2,1), mgp = c(1.5,.6,0), mar = .1 + c(4,4,2,1)) + for(N in c(50, 200)) { + n <- 1:N + sdn <- c(1, sapply(n[-1], function(m)sd(Sample(m)))) + r <- cbind(Qn = sapply(n, function(m)Qn(Sample(m))), + Sn = sapply(n, function(m)Sn(Sample(m)))) / sdn + matplot(n, r, type = 'b', col = 2:3, lty = 1, ylab = "Qn & Sn", + main = "Qn(Sample(n)) & Sn(..) [consistently scaled]") + legend(.85*N, 0.4, c("Qn()", "Sn()"), col = 2:3, lty = 1, + pch = paste(1:2)) + abline(h=1, col = "gray", lty = 2) + } + par(op) + + ## Hmm, the above does not look 100% consistent to *my* eyes... + ## Investigate: + + matplot(n, r, ylim = c(0.9, 1.1), type = 'b', col = 2:3, lty = 1) + abline(h=1, col = "gray", lty = 2) + + matplot(n, r^2, ylim = c(0.7, 1.3), type = 'b', col = 2:3, lty = 1) + abline(h=1, col = "gray", lty = 2) + } > rownames(r) <- paste(n) > r Qn Sn 1 0.0000000 0.0000000 2 1.2533141 1.2531372 3 2.2050485 2.2075026 4 0.9586584 0.9576951 5 1.4121898 1.2148172 6 1.2537914 1.0929275 7 1.2502496 1.1724933 8 1.1183680 1.2179731 9 1.1839692 1.2193954 10 1.0622839 1.1303545 11 1.1838478 1.2373579 12 1.1084824 1.0441905 13 1.1783568 1.0461282 14 1.0280173 1.1381402 15 1.1300531 1.1455153 16 1.0975640 0.9808160 17 1.1069023 0.9852123 18 1.0668214 1.0624166 19 1.1170674 1.0693949 20 1.0632780 1.0099979 21 1.0793114 1.0373703 22 1.1081293 1.0148426 23 1.0980137 1.0210812 24 1.0760013 1.0482209 25 1.0712945 1.0508545 26 1.0655762 1.0784132 27 1.0480275 1.1022088 28 1.0700198 1.0122606 29 1.1094825 1.0152751 30 1.0522957 1.0597715 31 1.0743482 1.0637003 32 1.0853126 0.9856182 33 1.0608374 0.9887047 34 1.0651112 1.0292893 35 1.0691718 1.0330266 36 1.0500552 1.0155475 37 1.0400758 1.0312545 38 1.0682034 1.0053387 39 1.0724800 1.0121787 40 1.0455886 1.0250685 41 1.0609770 1.0273275 42 1.0643505 1.0578421 43 1.0450362 1.0613117 44 1.0639112 1.0045793 45 1.0524585 1.0068645 46 1.0466486 1.0361589 47 1.0540086 1.0388060 48 1.0555688 0.9875950 49 1.0459186 0.9898554 50 1.0582257 1.0173856 51 1.0530772 1.0199364 52 1.0487478 1.0181890 53 1.0394412 1.0292047 54 1.0534751 1.0054990 55 1.0489389 1.0161168 56 1.0431361 1.0154413 57 1.0579340 1.0172547 58 1.0455372 1.0401705 59 1.0421755 1.0422140 60 1.0488290 1.0011312 61 1.0348723 1.0029317 62 1.0358527 1.0247613 63 1.0469162 1.0267544 64 1.0328647 0.9886494 65 1.0373891 0.9926067 66 1.0405167 1.0112521 67 1.0468842 1.0131881 68 1.0315401 1.0197079 69 1.0406152 1.0239818 70 1.0374683 1.0099878 71 1.0339173 1.0182336 72 1.0356644 1.0101629 73 1.0339045 1.0116582 74 1.0262912 1.0297406 75 1.0338120 1.0313730 76 1.0304585 0.9991718 77 1.0294667 1.0006501 78 1.0353933 1.0180455 79 1.0315025 1.0196432 80 1.0290092 0.9917780 81 1.0277752 0.9987728 82 1.0261869 1.0075101 83 1.0259820 1.0090701 84 1.0300727 1.0166174 85 1.0317531 1.0178957 86 1.0224292 1.0128137 87 1.0263644 1.0195537 88 1.0314665 1.0068284 89 1.0245860 1.0080956 90 1.0256323 1.0230268 91 1.0326019 1.0243849 92 1.0243610 0.9979076 93 1.0265952 0.9991594 94 1.0257519 1.0136175 95 1.0276671 1.0149505 96 1.0263910 0.9970979 97 1.0268222 1.0029782 98 1.0251603 1.0049884 99 1.0209779 1.0062947 100 1.0259809 1.0126494 101 1.0241265 1.0137579 102 1.0216934 1.0147559 103 1.0276353 1.0204551 104 1.0249992 1.0045303 105 1.0239035 1.0056281 106 1.0261457 1.0183429 107 1.0221293 1.0195053 108 1.0228166 0.9970240 109 1.0235131 0.9981087 110 1.0201276 1.0104780 111 1.0223323 1.0116215 112 1.0247535 1.0009575 113 1.0238656 1.0060296 114 1.0225941 1.0031735 115 1.0210953 1.0042970 116 1.0245629 1.0097845 117 1.0213179 1.0107616 118 1.0224002 1.0161725 119 1.0262676 1.0211094 120 1.0206794 1.0028503 121 1.0220803 1.0038178 122 1.0207793 1.0148890 123 1.0212114 1.0159048 124 1.0207892 0.9963715 125 1.0223094 0.9973281 126 1.0179975 1.0081359 127 1.0199437 1.0091369 128 1.0190798 1.0038855 129 1.0220290 1.0083447 130 1.0171853 1.0018045 131 1.0209487 1.0031645 132 1.0193469 1.0076187 133 1.0184645 1.0084916 134 1.0205855 1.0172514 135 1.0180659 1.0194240 136 1.0158370 1.0015684 137 1.0174794 1.0024329 138 1.0162069 1.0122365 139 1.0160160 1.0131386 140 1.0177024 0.9958697 141 1.0201002 0.9967251 142 1.0176871 1.0063214 143 1.0160048 1.0072116 144 1.0180068 1.0061828 145 1.0154433 1.0101613 146 1.0159359 1.0016303 147 1.0189714 1.0055554 148 1.0151345 1.0059238 149 1.0165125 1.0067123 150 1.0172866 1.0156821 151 1.0156318 1.0165025 152 1.0164800 1.0005579 153 1.0174555 1.0013392 154 1.0161708 1.0101356 155 1.0138415 1.0109467 156 1.0171954 0.9954718 157 1.0182371 0.9964004 158 1.0151745 1.0048742 159 1.0159541 1.0056757 160 1.0157937 1.0080333 161 1.0135156 1.0103514 162 1.0173009 1.0039306 163 1.0163459 1.0074784 164 1.0139466 1.0045613 165 1.0146181 1.0052800 166 1.0153397 1.0133993 167 1.0153876 1.0141440 168 1.0162602 0.9997410 169 1.0176006 1.0004534 170 1.0163195 1.0084303 171 1.0136678 1.0091671 172 1.0161586 0.9957005 173 1.0140386 0.9989947 174 1.0150578 1.0036931 175 1.0178367 1.0044219 176 1.0145431 1.0080196 177 1.0152377 1.0086849 178 1.0173132 1.0058220 179 1.0147503 1.0090587 180 1.0141714 1.0034420 181 1.0152680 1.0041022 182 1.0144619 1.0115183 183 1.0131775 1.0122000 184 1.0160414 0.9990667 185 1.0165892 0.9997214 186 1.0121585 1.0070184 187 1.0146730 1.0076934 188 1.0135477 0.9981359 189 1.0131875 1.0011603 190 1.0143600 1.0027108 191 1.0146898 1.0033790 192 1.0115525 1.0066812 193 1.0137139 1.0072961 194 1.0120277 1.0074046 195 1.0135076 1.0103804 196 1.0127864 1.0025061 197 1.0130898 1.0031165 198 1.0118025 1.0099416 199 1.0126815 1.0105700 200 1.0124098 0.9985008 > > proc.time() user system elapsed 0.31 0.07 0.37