R Under development (unstable) (2024-08-15 r87022 ucrt) -- "Unsuffered Consequences" Copyright (C) 2024 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > ## Copyright (C) 2012, 2015 Marius Hofert, Ivan Kojadinovic, Martin Maechler, and Jun Yan > ## > ## This program is free software; you can redistribute it and/or modify it under > ## the terms of the GNU General Public License as published by the Free Software > ## Foundation; either version 3 of the License, or (at your option) any later > ## version. > ## > ## This program is distributed in the hope that it will be useful, but WITHOUT > ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS > ## FOR A PARTICULAR PURPOSE. See the GNU General Public License for more > ## details. > ## > ## You should have received a copy of the GNU General Public License along with > ## this program; if not, see . > > ###__ Estimation methods for Archimedean Copulas __ > > require(copula) Loading required package: copula > > (doExtras <- copula:::doExtras()) [1] FALSE > > ## From source(system.file("test-tools.R", package = "Matrix")) : > showProc.time <- local({ + pct <- proc.time() + function() { ## CPU elapsed __since last called__ + ot <- pct ; pct <<- proc.time() + cat('Time elapsed: ', (pct - ot)[1:3],'\n') + } + }) > > ### We have other implicit estimation "tests" in demo(estimation.gof) > ## i.e., ../demo/estimation.gof.R > ## ~~~~~~~~~~~~~~~~~~~~~~~~ > > ## Frank & mle ---- log-density is *REALLY* not good enough: > ## ----------- > ## at least not for "GoF" where we may have "tail-dependent" > ## observations which are highly improbable under Frank and would > ## "need" a very large theta : > > p.log.f <- function(t0, t1, n.th = 256, + thet.vec = seq(t0, t1, length.out= n.th), + d, u0, uu = rbind(rep(u0, length.out = d)), + cop = "Frank", legend.x = "bottomright", + col = mapply(adjustcolor, col=c(1,3,2,4,1,2), alpha=c(2, 5:7, 4,2)/10), + lwd = c(5,4,3,1,5,5), lty = 1:6) + { + stopifnot(d == as.integer(d), length(d) == 1, d >= 2, + NCOL(uu) == d, is.numeric(thet.vec), length(thet.vec) > 1) + cop <- onacopulaL(getAcop(cop), list(NA, 1:d)) + + stopifnot(is.function(.f. <- cop@copula@dacopula)) + ## meths <- c("negI-s-Stirling", "negI-s-Eulerian") + meths <- local({m <- eval(formals(polylog)$method);m[grep("^negI-s", m)]}) + l.arr <- sapply(meths, function(meth) + { + sapply(c(FALSE,TRUE), function(Li.log) + sapply(thet.vec, .f., u = uu, log=TRUE, + method = meth, Li.log.arg=Li.log)) + }, simplify = "array") + dim(l.arr) <- local({D <- dim(l.arr); c(D[1], prod(D[-1]))}) + ## n x 4 + nms <- c(t(outer(substring(meths, 8), + paste("Li_n(", c("", "log(arg)"), ")", sep=""), + ##paste("Li.log.A=",c(FALSE,TRUE),sep=""), + paste, sep="_"))) + colnames(l.arr) <- nms + + tit <- bquote(.(cop@copula@name)*" copula - log density " ~~ + log~f(theta, bold(u)[d==.(d)])) + matplot(thet.vec, l.arr, type = "l", lwd=lwd, col=col, lty=lty, + main = tit, xlab=expression(theta), ylab = expression(log~f(theta, .))) + legend(legend.x, nms, col=col, lwd=lwd, lty=lty, + title = "polylog(.., method = \"*\")", inset = .01) + mkU <- function(u, n1 = 4, n2 = 2) { + n <- length(u <- c(u)) + cu <- if(n > (n1+n2)) ## use "..." + c(format(head(u, n1)), "...", format(tail(u, n2))) else format(u) + paste("(", paste(cu, collapse = ", "), ")", sep="") + } + mtext(bquote(bold(u)[d==.(d)] == .(mkU(uu))), line=1/8) + invisible(structure(cbind(theta = thet.vec, l.arr), + u = uu, dacopula = .f.)) + } > > unattr <- function(obj) { mostattributes(obj) <- NULL ; obj } > show. <- function(pmat) unattr(pmat)[, c(1,3,5,7)] > > if(!dev.interactive(orNone=TRUE)) pdf("estim-ex.pdf") > > th <- seq(1,60, by=1/4) > show.(m1 <- p.log.f(thet.vec = th, d = 5, u0 = .987)) theta Stirling_Li_n(log(arg)) Eulerian_Li_n(log(arg)) [1,] 1.00 4.071523 4.071523 [2,] 1.25 4.586197 4.586197 [3,] 1.50 5.038015 5.038015 [4,] 1.75 5.441197 5.441197 [5,] 2.00 5.805090 5.805090 [6,] 2.25 6.136270 6.136270 [7,] 2.50 6.439601 6.439601 [8,] 2.75 6.718831 6.718831 [9,] 3.00 6.976946 6.976946 [10,] 3.25 7.216378 7.216378 [11,] 3.50 7.439159 7.439159 [12,] 3.75 7.647008 7.647008 [13,] 4.00 7.841402 7.841402 [14,] 4.25 8.023624 8.023624 [15,] 4.50 8.194793 8.194793 [16,] 4.75 8.355900 8.355900 [17,] 5.00 8.507818 8.507818 [18,] 5.25 8.651327 8.651327 [19,] 5.50 8.787122 8.787122 [20,] 5.75 8.915827 8.915827 [21,] 6.00 9.038000 9.038000 [22,] 6.25 9.154146 9.154146 [23,] 6.50 9.264720 9.264720 [24,] 6.75 9.370133 9.370133 [25,] 7.00 9.470757 9.470757 [26,] 7.25 9.566930 9.566930 [27,] 7.50 9.658960 9.658960 [28,] 7.75 9.747128 9.747128 [29,] 8.00 9.831687 9.831687 [30,] 8.25 9.912872 9.912872 [31,] 8.50 9.990896 9.990896 [32,] 8.75 10.065956 10.065956 [33,] 9.00 10.138231 10.138231 [34,] 9.25 10.207887 10.207887 [35,] 9.50 10.275075 10.275075 [36,] 9.75 10.339938 10.339938 [37,] 10.00 10.402604 10.402604 [38,] 10.25 10.463195 10.463195 [39,] 10.50 10.521822 10.521822 [40,] 10.75 10.578587 10.578587 [41,] 11.00 10.633589 10.633589 [42,] 11.25 10.686915 10.686915 [43,] 11.50 10.738650 10.738650 [44,] 11.75 10.788871 10.788871 [45,] 12.00 10.837652 10.837652 [46,] 12.25 10.885060 10.885060 [47,] 12.50 10.931159 10.931159 [48,] 12.75 10.976009 10.976009 [49,] 13.00 11.019666 11.019666 [50,] 13.25 11.062183 11.062183 [51,] 13.50 11.103609 11.103609 [52,] 13.75 11.143991 11.143991 [53,] 14.00 11.183374 11.183374 [54,] 14.25 11.221798 11.221798 [55,] 14.50 11.259304 11.259304 [56,] 14.75 11.295927 11.295927 [57,] 15.00 11.331704 11.331704 [58,] 15.25 11.366667 11.366667 [59,] 15.50 11.400849 11.400849 [60,] 15.75 11.434278 11.434278 [61,] 16.00 11.466983 11.466983 [62,] 16.25 11.498991 11.498991 [63,] 16.50 11.530328 11.530328 [64,] 16.75 11.561018 11.561018 [65,] 17.00 11.591085 11.591085 [66,] 17.25 11.620549 11.620549 [67,] 17.50 11.649433 11.649433 [68,] 17.75 11.677756 11.677756 [69,] 18.00 11.705537 11.705537 [70,] 18.25 11.732795 11.732795 [71,] 18.50 11.759548 11.759548 [72,] 18.75 11.785811 11.785811 [73,] 19.00 11.811601 11.811601 [74,] 19.25 11.836932 11.836932 [75,] 19.50 11.861821 11.861821 [76,] 19.75 11.886280 11.886280 [77,] 20.00 11.910323 11.910323 [78,] 20.25 11.933963 11.933963 [79,] 20.50 11.957212 11.957212 [80,] 20.75 11.980082 11.980082 [81,] 21.00 12.002585 12.002585 [82,] 21.25 12.024730 12.024730 [83,] 21.50 12.046530 12.046530 [84,] 21.75 12.067993 12.067993 [85,] 22.00 12.089130 12.089130 [86,] 22.25 12.109950 12.109950 [87,] 22.50 12.130461 12.130461 [88,] 22.75 12.150673 12.150673 [89,] 23.00 12.170593 12.170593 [90,] 23.25 12.190231 12.190231 [91,] 23.50 12.209593 12.209593 [92,] 23.75 12.228686 12.228686 [93,] 24.00 12.247519 12.247519 [94,] 24.25 12.266098 12.266098 [95,] 24.50 12.284430 12.284430 [96,] 24.75 12.302521 12.302521 [97,] 25.00 12.320377 12.320377 [98,] 25.25 12.338004 12.338004 [99,] 25.50 12.355409 12.355409 [100,] 25.75 12.372597 12.372597 [101,] 26.00 12.389573 12.389573 [102,] 26.25 12.406342 12.406342 [103,] 26.50 12.422910 12.422910 [104,] 26.75 12.439282 12.439282 [105,] 27.00 12.455462 12.455462 [106,] 27.25 12.471455 12.471455 [107,] 27.50 12.487266 12.487266 [108,] 27.75 12.502898 12.502898 [109,] 28.00 12.518356 12.518356 [110,] 28.25 12.533644 12.533644 [111,] 28.50 12.548766 12.548766 [112,] 28.75 12.563725 12.563725 [113,] 29.00 12.578527 12.578527 [114,] 29.25 12.593173 12.593173 [115,] 29.50 12.607667 12.607667 [116,] 29.75 12.622014 12.622014 [117,] 30.00 12.636216 12.636216 [118,] 30.25 12.650276 12.650276 [119,] 30.50 12.664197 12.664197 [120,] 30.75 12.677984 12.677984 [121,] 31.00 12.691637 12.691637 [122,] 31.25 12.705161 12.705161 [123,] 31.50 12.718557 12.718557 [124,] 31.75 12.731830 12.731830 [125,] 32.00 12.744981 12.744981 [126,] 32.25 12.758012 12.758012 [127,] 32.50 12.770927 12.770927 [128,] 32.75 12.783727 12.783727 [129,] 33.00 12.796416 12.796416 [130,] 33.25 12.808995 12.808995 [131,] 33.50 12.821466 12.821466 [132,] 33.75 12.833832 12.833832 [133,] 34.00 12.846095 12.846095 [134,] 34.25 12.858256 12.858256 [135,] 34.50 12.870319 12.870319 [136,] 34.75 12.882284 12.882284 [137,] 35.00 12.894153 12.894153 [138,] 35.25 12.905929 12.905929 [139,] 35.50 12.917614 12.917614 [140,] 35.75 12.929208 12.929208 [141,] 36.00 12.940714 12.940714 [142,] 36.25 12.952134 12.952134 [143,] 36.50 12.963468 12.963468 [144,] 36.75 12.974719 12.974719 [145,] 37.00 12.985888 12.985888 [146,] 37.25 12.996977 12.996977 [147,] 37.50 13.007986 13.007986 [148,] 37.75 13.018919 13.018919 [149,] 38.00 13.029775 13.029775 [150,] 38.25 13.040556 13.040556 [151,] 38.50 13.051264 13.051264 [152,] 38.75 13.061900 13.061900 [153,] 39.00 13.072465 13.072465 [154,] 39.25 13.082961 13.082961 [155,] 39.50 13.093388 13.093388 [156,] 39.75 13.103748 13.103748 [157,] 40.00 13.114041 13.114041 [158,] 40.25 13.124270 13.124270 [159,] 40.50 13.134436 13.134436 [160,] 40.75 13.144538 13.144538 [161,] 41.00 13.154579 13.154579 [162,] 41.25 13.164559 13.164559 [163,] 41.50 13.174480 13.174480 [164,] 41.75 13.184341 13.184341 [165,] 42.00 13.194146 13.194146 [166,] 42.25 13.203893 13.203893 [167,] 42.50 13.213585 13.213585 [168,] 42.75 13.223222 13.223222 [169,] 43.00 13.232804 13.232804 [170,] 43.25 13.242334 13.242334 [171,] 43.50 13.251811 13.251811 [172,] 43.75 13.261237 13.261237 [173,] 44.00 13.270612 13.270612 [174,] 44.25 13.279938 13.279938 [175,] 44.50 13.289214 13.289214 [176,] 44.75 13.298442 13.298442 [177,] 45.00 13.307622 13.307622 [178,] 45.25 13.316756 13.316756 [179,] 45.50 13.325843 13.325843 [180,] 45.75 13.334885 13.334885 [181,] 46.00 13.343882 13.343882 [182,] 46.25 13.352835 13.352835 [183,] 46.50 13.361745 13.361745 [184,] 46.75 13.370612 13.370612 [185,] 47.00 13.379436 13.379436 [186,] 47.25 13.388220 13.388220 [187,] 47.50 13.396962 13.396962 [188,] 47.75 13.405664 13.405664 [189,] 48.00 13.414326 13.414326 [190,] 48.25 13.422949 13.422949 [191,] 48.50 13.431533 13.431533 [192,] 48.75 13.440080 13.440080 [193,] 49.00 13.448589 13.448589 [194,] 49.25 13.457061 13.457061 [195,] 49.50 13.465496 13.465496 [196,] 49.75 13.473896 13.473896 [197,] 50.00 13.482260 13.482260 [198,] 50.25 13.490589 13.490589 [199,] 50.50 13.498883 13.498883 [200,] 50.75 13.507144 13.507144 [201,] 51.00 13.515371 13.515371 [202,] 51.25 13.523565 13.523565 [203,] 51.50 13.531726 13.531726 [204,] 51.75 13.539855 13.539855 [205,] 52.00 13.547952 13.547952 [206,] 52.25 13.556018 13.556018 [207,] 52.50 13.564053 13.564053 [208,] 52.75 13.572058 13.572058 [209,] 53.00 13.580032 13.580032 [210,] 53.25 13.587976 13.587976 [211,] 53.50 13.595892 13.595892 [212,] 53.75 13.603778 13.603778 [213,] 54.00 13.611635 13.611635 [214,] 54.25 13.619465 13.619465 [215,] 54.50 13.627266 13.627266 [216,] 54.75 13.635040 13.635040 [217,] 55.00 13.642787 13.642787 [218,] 55.25 13.650507 13.650507 [219,] 55.50 13.658200 13.658200 [220,] 55.75 13.665867 13.665867 [221,] 56.00 13.673508 13.673508 [222,] 56.25 13.681124 13.681124 [223,] 56.50 13.688715 13.688715 [224,] 56.75 13.696280 13.696280 [225,] 57.00 13.703821 13.703821 [226,] 57.25 13.711338 13.711338 [227,] 57.50 13.718830 13.718830 [228,] 57.75 13.726299 13.726299 [229,] 58.00 13.733744 13.733744 [230,] 58.25 13.741166 13.741166 [231,] 58.50 13.748566 13.748566 [232,] 58.75 13.755942 13.755942 [233,] 59.00 13.763296 13.763296 [234,] 59.25 13.770628 13.770628 [235,] 59.50 13.777938 13.777938 [236,] 59.75 13.785227 13.785227 [237,] 60.00 13.792494 13.792494 asymp-w_Li_n(log(arg)) [1,] 4.071520 [2,] 4.586196 [3,] 5.038015 [4,] 5.441197 [5,] 5.805090 [6,] 6.136270 [7,] 6.439601 [8,] 6.718831 [9,] 6.976946 [10,] 7.216378 [11,] 7.439159 [12,] 7.647008 [13,] 7.841402 [14,] 8.023624 [15,] 8.194793 [16,] 8.355900 [17,] 8.507818 [18,] 8.651327 [19,] 8.787122 [20,] 8.915827 [21,] 9.038000 [22,] 9.154146 [23,] 9.264720 [24,] 9.370133 [25,] 9.470757 [26,] 9.566930 [27,] 9.658960 [28,] 9.747128 [29,] 9.831687 [30,] 9.912872 [31,] 9.990896 [32,] 10.065956 [33,] 10.138231 [34,] 10.207887 [35,] 10.275075 [36,] 10.339938 [37,] 10.402604 [38,] 10.463195 [39,] 10.521822 [40,] 10.578587 [41,] 10.633589 [42,] 10.686915 [43,] 10.738650 [44,] 10.788871 [45,] 10.837652 [46,] 10.885060 [47,] 10.931159 [48,] 10.976009 [49,] 11.019666 [50,] 11.062183 [51,] 11.103609 [52,] 11.143991 [53,] 11.183374 [54,] 11.221798 [55,] 11.259304 [56,] 11.295927 [57,] 11.331704 [58,] 11.366667 [59,] 11.400849 [60,] 11.434278 [61,] 11.466983 [62,] 11.498991 [63,] 11.530328 [64,] 11.561018 [65,] 11.591085 [66,] 11.620549 [67,] 11.649433 [68,] 11.677756 [69,] 11.705537 [70,] 11.732795 [71,] 11.759548 [72,] 11.785811 [73,] 11.811601 [74,] 11.836932 [75,] 11.861821 [76,] 11.886280 [77,] 11.910323 [78,] 11.933963 [79,] 11.957212 [80,] 11.980082 [81,] 12.002585 [82,] 12.024730 [83,] 12.046530 [84,] 12.067993 [85,] 12.089130 [86,] 12.109950 [87,] 12.130461 [88,] 12.150673 [89,] 12.170593 [90,] 12.190231 [91,] 12.209593 [92,] 12.228686 [93,] 12.247519 [94,] 12.266098 [95,] 12.284430 [96,] 12.302521 [97,] 12.320377 [98,] 12.338004 [99,] 12.355409 [100,] 12.372597 [101,] 12.389573 [102,] 12.406342 [103,] 12.422910 [104,] 12.439282 [105,] 12.455462 [106,] 12.471455 [107,] 12.487266 [108,] 12.502898 [109,] 12.518356 [110,] 12.533644 [111,] 12.548766 [112,] 12.563725 [113,] 12.578527 [114,] 12.593173 [115,] 12.607667 [116,] 12.622014 [117,] 12.636216 [118,] 12.650276 [119,] 12.664197 [120,] 12.677984 [121,] 12.691637 [122,] 12.705161 [123,] 12.718557 [124,] 12.731830 [125,] 12.744981 [126,] 12.758012 [127,] 12.770927 [128,] 12.783727 [129,] 12.796416 [130,] 12.808995 [131,] 12.821466 [132,] 12.833832 [133,] 12.846095 [134,] 12.858256 [135,] 12.870319 [136,] 12.882284 [137,] 12.894153 [138,] 12.905929 [139,] 12.917614 [140,] 12.929208 [141,] 12.940714 [142,] 12.952134 [143,] 12.963468 [144,] 12.974719 [145,] 12.985888 [146,] 12.996977 [147,] 13.007986 [148,] 13.018919 [149,] 13.029775 [150,] 13.040556 [151,] 13.051264 [152,] 13.061900 [153,] 13.072465 [154,] 13.082961 [155,] 13.093388 [156,] 13.103748 [157,] 13.114041 [158,] 13.124270 [159,] 13.134436 [160,] 13.144538 [161,] 13.154579 [162,] 13.164559 [163,] 13.174480 [164,] 13.184341 [165,] 13.194146 [166,] 13.203893 [167,] 13.213585 [168,] 13.223222 [169,] 13.232804 [170,] 13.242334 [171,] 13.251811 [172,] 13.261237 [173,] 13.270612 [174,] 13.279938 [175,] 13.289214 [176,] 13.298442 [177,] 13.307622 [178,] 13.316756 [179,] 13.325843 [180,] 13.334885 [181,] 13.343882 [182,] 13.352835 [183,] 13.361745 [184,] 13.370612 [185,] 13.379436 [186,] 13.388220 [187,] 13.396962 [188,] 13.405664 [189,] 13.414326 [190,] 13.422949 [191,] 13.431533 [192,] 13.440080 [193,] 13.448589 [194,] 13.457061 [195,] 13.465496 [196,] 13.473896 [197,] 13.482260 [198,] 13.490589 [199,] 13.498883 [200,] 13.507144 [201,] 13.515371 [202,] 13.523565 [203,] 13.531726 [204,] 13.539855 [205,] 13.547952 [206,] 13.556018 [207,] 13.564053 [208,] 13.572058 [209,] 13.580032 [210,] 13.587976 [211,] 13.595892 [212,] 13.603778 [213,] 13.611635 [214,] 13.619465 [215,] 13.627266 [216,] 13.635040 [217,] 13.642787 [218,] 13.650507 [219,] 13.658200 [220,] 13.665867 [221,] 13.673508 [222,] 13.681124 [223,] 13.688715 [224,] 13.696280 [225,] 13.703821 [226,] 13.711338 [227,] 13.718830 [228,] 13.726299 [229,] 13.733744 [230,] 13.741166 [231,] 13.748566 [232,] 13.755942 [233,] 13.763296 [234,] 13.770628 [235,] 13.777938 [236,] 13.785227 [237,] 13.792494 > ## wow! asymptotic is quite good even here! > > ## using MC (instead of polylog(.)) also ``works'' {the function is noisy ..}: > l.th.MC <- sapply(th, attr(m1, "dacopula"), u = attr(m1, "u"), + n.MC = 2000, log=TRUE) ## 2000: be speedy in test > lines(th, l.th.MC, lwd=3, col=adjustcolor("sandybrown", .5)) > legend("right", "n.MC = 2000", lwd=3, col=adjustcolor("sandybrown", .5), + inset=.01) > rm(th) > > showProc.time() Time elapsed: 0.5 0 0.5 > > ## Extend the range: > show.(m2 <- p.log.f(1, 200, n.th=401, d = 5, u0 = .987)) theta Stirling_Li_n(log(arg)) Eulerian_Li_n(log(arg)) [1,] 1.0000 4.071523 4.071523 [2,] 1.4975 5.033756 5.033756 [3,] 1.9950 5.798152 5.798152 [4,] 2.4925 6.430868 6.430868 [5,] 2.9900 6.966994 6.966994 [6,] 3.4875 7.428388 7.428388 [7,] 3.9850 7.830094 7.830094 [8,] 4.4825 8.183150 8.183150 [9,] 4.9800 8.495984 8.495984 [10,] 5.4775 8.775200 8.775200 [11,] 5.9750 9.026062 9.026062 [12,] 6.4725 9.252817 9.252817 [13,] 6.9700 9.458923 9.458923 [14,] 7.4675 9.647221 9.647221 [15,] 7.9650 9.820057 9.820057 [16,] 8.4625 9.979386 9.979386 [17,] 8.9600 10.126847 10.126847 [18,] 9.4575 10.263821 10.263821 [19,] 9.9550 10.391481 10.391481 [20,] 10.4525 10.510829 10.510829 [21,] 10.9500 10.622725 10.622725 [22,] 11.4475 10.727914 10.727914 [23,] 11.9450 10.827040 10.827040 [24,] 12.4425 10.920669 10.920669 [25,] 12.9400 11.009294 11.009294 [26,] 13.4375 11.093352 11.093352 [27,] 13.9350 11.173228 11.173228 [28,] 14.4325 11.249266 11.249266 [29,] 14.9300 11.321770 11.321770 [30,] 15.4275 11.391015 11.391015 [31,] 15.9250 11.457246 11.457246 [32,] 16.4225 11.520684 11.520684 [33,] 16.9200 11.581530 11.581530 [34,] 17.4175 11.639964 11.639964 [35,] 17.9150 11.696151 11.696151 [36,] 18.4125 11.750241 11.750241 [37,] 18.9100 11.802370 11.802370 [38,] 19.4075 11.852663 11.852663 [39,] 19.9050 11.901235 11.901235 [40,] 20.4025 11.948191 11.948191 [41,] 20.9000 11.993627 11.993627 [42,] 21.3975 12.037633 12.037633 [43,] 21.8950 12.080292 12.080292 [44,] 22.3925 12.121678 12.121678 [45,] 22.8900 12.161864 12.161864 [46,] 23.3875 12.200913 12.200913 [47,] 23.8850 12.238888 12.238888 [48,] 24.3825 12.275844 12.275844 [49,] 24.8800 12.311835 12.311835 [50,] 25.3775 12.346908 12.346908 [51,] 25.8750 12.381111 12.381111 [52,] 26.3725 12.414486 12.414486 [53,] 26.8700 12.447072 12.447072 [54,] 27.3675 12.478909 12.478909 [55,] 27.8650 12.510030 12.510030 [56,] 28.3625 12.540469 12.540469 [57,] 28.8600 12.570257 12.570257 [58,] 29.3575 12.599424 12.599424 [59,] 29.8550 12.627996 12.627996 [60,] 30.3525 12.656000 12.656000 [61,] 30.8500 12.683461 12.683461 [62,] 31.3475 12.710400 12.710400 [63,] 31.8450 12.736841 12.736841 [64,] 32.3425 12.762804 12.762804 [65,] 32.8400 12.788308 12.788308 [66,] 33.3375 12.813372 12.813372 [67,] 33.8350 12.838013 12.838013 [68,] 34.3325 12.862248 12.862248 [69,] 34.8300 12.886092 12.886092 [70,] 35.3275 12.909561 12.909561 [71,] 35.8250 12.932669 12.932669 [72,] 36.3225 12.955429 12.955429 [73,] 36.8200 12.977854 12.977854 [74,] 37.3175 12.999957 12.999957 [75,] 37.8150 13.021748 13.021748 [76,] 38.3125 13.043240 13.043240 [77,] 38.8100 13.064442 13.064442 [78,] 39.3075 13.085365 13.085365 [79,] 39.8050 13.106018 13.106018 [80,] 40.3025 13.126410 13.126410 [81,] 40.8000 13.146551 13.146551 [82,] 41.2975 13.166448 13.166448 [83,] 41.7950 13.186110 13.186110 [84,] 42.2925 13.205545 13.205545 [85,] 42.7900 13.224758 13.224758 [86,] 43.2875 13.243759 13.243759 [87,] 43.7850 13.262553 13.262553 [88,] 44.2825 13.281146 13.281146 [89,] 44.7800 13.299546 13.299546 [90,] 45.2775 13.317758 13.317758 [91,] 45.7750 13.335787 13.335787 [92,] 46.2725 13.353639 13.353639 [93,] 46.7700 13.371319 13.371319 [94,] 47.2675 13.388833 13.388833 [95,] 47.7650 13.406185 13.406185 [96,] 48.2625 13.423379 13.423379 [97,] 48.7600 13.440421 13.440421 [98,] 49.2575 13.457314 13.457314 [99,] 49.7550 13.474063 13.474063 [100,] 50.2525 13.490672 13.490672 [101,] 50.7500 13.507144 13.507144 [102,] 51.2475 13.523483 13.523483 [103,] 51.7450 13.539693 13.539693 [104,] 52.2425 13.555777 13.555777 [105,] 52.7400 13.571738 13.571738 [106,] 53.2375 13.587580 13.587580 [107,] 53.7350 13.603305 13.603305 [108,] 54.2325 13.618918 13.618918 [109,] 54.7300 13.634419 13.634419 [110,] 55.2275 13.649813 13.649813 [111,] 55.7250 13.665102 13.665102 [112,] 56.2225 13.680288 13.680288 [113,] 56.7200 13.695374 13.695374 [114,] 57.2175 13.710362 13.710362 [115,] 57.7150 13.725255 13.725255 [116,] 58.2125 13.740054 13.740054 [117,] 58.7100 13.754763 13.754763 [118,] 59.2075 13.769383 13.769383 [119,] 59.7050 13.783916 13.783916 [120,] 60.2025 13.798364 13.798364 [121,] 60.7000 13.812729 13.812729 [122,] 61.1975 13.827013 13.827013 [123,] 61.6950 13.841218 13.841218 [124,] 62.1925 13.855345 13.855345 [125,] 62.6900 13.869396 13.869396 [126,] 63.1875 13.883373 13.883373 [127,] 63.6850 13.897276 13.897276 [128,] 64.1825 13.911109 13.911109 [129,] 64.6800 13.924871 13.924871 [130,] 65.1775 13.938566 13.938566 [131,] 65.6750 13.952193 13.952193 [132,] 66.1725 13.965754 13.965754 [133,] 66.6700 13.979252 13.979252 [134,] 67.1675 13.992686 13.992686 [135,] 67.6650 14.006058 14.006058 [136,] 68.1625 14.019370 14.019370 [137,] 68.6600 14.032622 14.032622 [138,] 69.1575 14.045815 14.045815 [139,] 69.6550 14.058952 14.058952 [140,] 70.1525 14.072032 14.072032 [141,] 70.6500 14.085057 14.085057 [142,] 71.1475 14.098027 14.098027 [143,] 71.6450 14.110945 14.110945 [144,] 72.1425 14.123810 14.123810 [145,] 72.6400 14.136623 14.136623 [146,] 73.1375 14.149387 14.149387 [147,] 73.6350 14.162100 14.162100 [148,] 74.1325 14.174765 14.174765 [149,] 74.6300 14.187382 14.187382 [150,] 75.1275 14.199952 14.199952 [151,] 75.6250 14.212475 14.212475 [152,] 76.1225 14.224953 14.224953 [153,] 76.6200 14.237386 14.237386 [154,] 77.1175 14.249775 14.249775 [155,] 77.6150 14.262120 14.262120 [156,] 78.1125 14.274423 14.274423 [157,] 78.6100 14.286683 14.286683 [158,] 79.1075 14.298902 14.298902 [159,] 79.6050 14.311080 14.311080 [160,] 80.1025 14.323218 14.323218 [161,] 80.6000 14.335317 14.335317 [162,] 81.0975 14.347376 14.347376 [163,] 81.5950 14.359397 14.359397 [164,] 82.0925 14.371380 14.371380 [165,] 82.5900 14.383325 14.383325 [166,] 83.0875 14.395234 14.395234 [167,] 83.5850 14.407107 14.407107 [168,] 84.0825 14.418943 14.418943 [169,] 84.5800 14.430745 14.430745 [170,] 85.0775 14.442511 14.442511 [171,] 85.5750 14.454243 14.454243 [172,] 86.0725 14.465941 14.465941 [173,] 86.5700 14.477606 14.477606 [174,] 87.0675 14.489237 14.489237 [175,] 87.5650 14.500836 14.500836 [176,] 88.0625 14.512403 14.512403 [177,] 88.5600 14.523938 14.523938 [178,] 89.0575 14.535442 14.535442 [179,] 89.5550 14.546914 14.546914 [180,] 90.0525 14.558356 14.558356 [181,] 90.5500 14.569767 14.569767 [182,] 91.0475 14.581149 14.581149 [183,] 91.5450 14.592501 14.592501 [184,] 92.0425 14.603824 14.603824 [185,] 92.5400 14.615118 14.615118 [186,] 93.0375 14.626383 14.626383 [187,] 93.5350 14.637620 14.637620 [188,] 94.0325 14.648829 14.648829 [189,] 94.5300 14.660011 14.660011 [190,] 95.0275 14.671165 14.671165 [191,] 95.5250 14.682292 14.682292 [192,] 96.0225 14.693393 14.693393 [193,] 96.5200 14.704466 14.704466 [194,] 97.0175 14.715514 14.715514 [195,] 97.5150 14.726536 14.726536 [196,] 98.0125 14.737532 14.737532 [197,] 98.5100 14.748503 14.748503 [198,] 99.0075 14.759448 14.759448 [199,] 99.5050 14.770369 14.770369 [200,] 100.0025 14.781265 14.781265 [201,] 100.5000 14.792137 14.792137 [202,] 100.9975 14.802984 14.802984 [203,] 101.4950 14.813807 14.813807 [204,] 101.9925 14.824607 14.824607 [205,] 102.4900 14.835383 14.835383 [206,] 102.9875 14.846136 14.846136 [207,] 103.4850 14.856865 14.856865 [208,] 103.9825 14.867572 14.867572 [209,] 104.4800 14.878256 14.878256 [210,] 104.9775 14.888918 14.888918 [211,] 105.4750 14.899557 14.899557 [212,] 105.9725 14.910174 14.910174 [213,] 106.4700 14.920769 14.920769 [214,] 106.9675 14.931343 14.931343 [215,] 107.4650 14.941894 14.941894 [216,] 107.9625 14.952425 14.952425 [217,] 108.4600 14.962934 14.962934 [218,] 108.9575 14.973422 14.973422 [219,] 109.4550 14.983889 14.983889 [220,] 109.9525 14.994336 14.994336 [221,] 110.4500 15.004761 15.004761 [222,] 110.9475 15.015167 15.015167 [223,] 111.4450 15.025552 15.025552 [224,] 111.9425 15.035917 15.035917 [225,] 112.4400 15.046262 15.046262 [226,] 112.9375 15.056587 15.056587 [227,] 113.4350 15.066892 15.066892 [228,] 113.9325 15.077178 15.077178 [229,] 114.4300 15.087445 15.087445 [230,] 114.9275 15.097692 15.097692 [231,] 115.4250 15.107920 15.107920 [232,] 115.9225 15.118129 15.118129 [233,] 116.4200 15.128319 15.128319 [234,] 116.9175 15.138490 15.138490 [235,] 117.4150 15.148643 15.148643 [236,] 117.9125 15.158777 15.158777 [237,] 118.4100 15.168892 15.168892 [238,] 118.9075 15.178990 15.178990 [239,] 119.4050 15.189069 15.189069 [240,] 119.9025 15.199130 15.199130 [241,] 120.4000 15.209173 15.209173 [242,] 120.8975 15.219198 15.219198 [243,] 121.3950 15.229206 15.229206 [244,] 121.8925 15.239196 15.239196 [245,] 122.3900 15.249168 15.249168 [246,] 122.8875 15.259123 15.259123 [247,] 123.3850 15.269060 15.269060 [248,] 123.8825 15.278981 15.278981 [249,] 124.3800 15.288884 15.288884 [250,] 124.8775 15.298770 15.298770 [251,] 125.3750 15.308639 15.308639 [252,] 125.8725 15.318491 15.318491 [253,] 126.3700 15.328326 15.328326 [254,] 126.8675 15.338145 15.338145 [255,] 127.3650 15.347947 15.347947 [256,] 127.8625 15.357733 15.357733 [257,] 128.3600 15.367502 15.367502 [258,] 128.8575 15.377255 15.377255 [259,] 129.3550 15.386991 15.386991 [260,] 129.8525 15.396712 15.396712 [261,] 130.3500 15.406416 15.406416 [262,] 130.8475 15.416104 15.416104 [263,] 131.3450 15.425776 15.425776 [264,] 131.8425 15.435433 15.435433 [265,] 132.3400 15.445073 15.445073 [266,] 132.8375 15.454698 15.454698 [267,] 133.3350 15.464307 15.464307 [268,] 133.8325 15.473901 15.473901 [269,] 134.3300 15.483479 15.483479 [270,] 134.8275 15.493041 15.493041 [271,] 135.3250 15.502589 15.502589 [272,] 135.8225 15.512121 15.512121 [273,] 136.3200 15.521637 15.521637 [274,] 136.8175 15.531139 15.531139 [275,] 137.3150 15.540625 15.540625 [276,] 137.8125 15.550096 15.550096 [277,] 138.3100 15.559553 15.559553 [278,] 138.8075 15.568994 15.568994 [279,] 139.3050 15.578420 15.578420 [280,] 139.8025 15.587832 15.587832 [281,] 140.3000 15.597229 15.597229 [282,] 140.7975 15.606611 15.606611 [283,] 141.2950 15.615979 15.615979 [284,] 141.7925 15.625332 15.625332 [285,] 142.2900 15.634670 15.634670 [286,] 142.7875 15.643994 15.643994 [287,] 143.2850 15.653304 15.653304 [288,] 143.7825 15.662599 15.662599 [289,] 144.2800 15.671880 15.671880 [290,] 144.7775 15.681147 15.681147 [291,] 145.2750 15.690400 15.690400 [292,] 145.7725 15.699638 15.699638 [293,] 146.2700 15.708862 15.708862 [294,] 146.7675 15.718073 15.718073 [295,] 147.2650 15.727269 15.727269 [296,] 147.7625 15.736451 15.736451 [297,] 148.2600 15.745619 15.745619 [298,] 148.7575 15.754774 15.754774 [299,] 149.2550 15.763915 15.763915 [300,] 149.7525 15.773042 15.773042 [301,] 150.2500 15.782155 15.782155 [302,] 150.7475 15.791255 15.791255 [303,] 151.2450 15.800341 15.800341 [304,] 151.7425 15.809413 15.809413 [305,] 152.2400 15.818472 15.818472 [306,] 152.7375 15.827517 15.827517 [307,] 153.2350 15.836549 15.836549 [308,] 153.7325 15.845568 15.845568 [309,] 154.2300 15.854573 15.854573 [310,] 154.7275 15.863565 15.863565 [311,] 155.2250 15.872543 15.872543 [312,] 155.7225 15.881509 15.881509 [313,] 156.2200 15.890461 15.890461 [314,] 156.7175 15.899400 15.899400 [315,] 157.2150 15.908326 15.908326 [316,] 157.7125 15.917238 15.917238 [317,] 158.2100 15.926138 15.926138 [318,] 158.7075 15.935025 15.935025 [319,] 159.2050 15.943899 15.943899 [320,] 159.7025 15.952759 15.952759 [321,] 160.2000 15.961607 15.961607 [322,] 160.6975 15.970442 15.970442 [323,] 161.1950 15.979265 15.979265 [324,] 161.6925 15.988074 15.988074 [325,] 162.1900 15.996871 15.996871 [326,] 162.6875 16.005655 16.005655 [327,] 163.1850 16.014426 16.014426 [328,] 163.6825 16.023185 16.023185 [329,] 164.1800 16.031931 16.031931 [330,] 164.6775 16.040665 16.040665 [331,] 165.1750 16.049386 16.049386 [332,] 165.6725 16.058094 16.058094 [333,] 166.1700 16.066790 16.066790 [334,] 166.6675 16.075474 16.075474 [335,] 167.1650 16.084145 16.084145 [336,] 167.6625 16.092803 16.092803 [337,] 168.1600 16.101450 16.101450 [338,] 168.6575 16.110084 16.110084 [339,] 169.1550 16.118706 16.118706 [340,] 169.6525 16.127315 16.127315 [341,] 170.1500 16.135913 16.135913 [342,] 170.6475 16.144498 16.144498 [343,] 171.1450 16.153071 16.153071 [344,] 171.6425 16.161632 16.161632 [345,] 172.1400 16.170180 16.170180 [346,] 172.6375 16.178717 16.178717 [347,] 173.1350 16.187242 16.187242 [348,] 173.6325 16.195754 16.195754 [349,] 174.1300 16.204255 16.204255 [350,] 174.6275 16.212743 16.212743 [351,] 175.1250 16.221220 16.221220 [352,] 175.6225 16.229685 16.229685 [353,] 176.1200 16.238138 16.238138 [354,] 176.6175 16.246579 16.246579 [355,] 177.1150 16.255008 16.255008 [356,] 177.6125 16.263426 16.263426 [357,] 178.1100 16.271832 16.271832 [358,] 178.6075 16.280226 16.280226 [359,] 179.1050 16.288608 16.288608 [360,] 179.6025 16.296979 16.296979 [361,] 180.1000 16.305338 16.305338 [362,] 180.5975 Inf 16.313685 [363,] 181.0950 Inf 16.322021 [364,] 181.5925 Inf 16.330345 [365,] 182.0900 Inf 16.338657 [366,] 182.5875 Inf 16.346959 [367,] 183.0850 Inf 16.355248 [368,] 183.5825 Inf 16.363526 [369,] 184.0800 Inf 16.371793 [370,] 184.5775 Inf 16.380048 [371,] 185.0750 Inf 16.388292 [372,] 185.5725 Inf 16.396524 [373,] 186.0700 Inf 16.404745 [374,] 186.5675 Inf 16.412955 [375,] 187.0650 Inf 16.421154 [376,] 187.5625 Inf 16.429341 [377,] 188.0600 Inf 16.437517 [378,] 188.5575 Inf 16.445681 [379,] 189.0550 Inf 16.453835 [380,] 189.5525 Inf 16.461977 [381,] 190.0500 Inf 16.470108 [382,] 190.5475 Inf 16.478228 [383,] 191.0450 Inf 16.486337 [384,] 191.5425 Inf 16.494435 [385,] 192.0400 Inf 16.502521 [386,] 192.5375 Inf 16.510597 [387,] 193.0350 Inf 16.518662 [388,] 193.5325 Inf 16.526715 [389,] 194.0300 Inf 16.534758 [390,] 194.5275 Inf 16.542790 [391,] 195.0250 Inf 16.550810 [392,] 195.5225 Inf 16.558820 [393,] 196.0200 Inf 16.566819 [394,] 196.5175 Inf 16.574807 [395,] 197.0150 Inf 16.582784 [396,] 197.5125 Inf 16.590751 [397,] 198.0100 Inf 16.598706 [398,] 198.5075 Inf 16.606651 [399,] 199.0050 Inf 16.614585 [400,] 199.5025 Inf 16.622508 [401,] 200.0000 Inf 16.630421 asymp-w_Li_n(log(arg)) [1,] 4.071520 [2,] 5.033756 [3,] 5.798152 [4,] 6.430868 [5,] 6.966994 [6,] 7.428388 [7,] 7.830094 [8,] 8.183150 [9,] 8.495984 [10,] 8.775200 [11,] 9.026062 [12,] 9.252817 [13,] 9.458923 [14,] 9.647221 [15,] 9.820057 [16,] 9.979386 [17,] 10.126847 [18,] 10.263821 [19,] 10.391481 [20,] 10.510829 [21,] 10.622725 [22,] 10.727914 [23,] 10.827040 [24,] 10.920669 [25,] 11.009294 [26,] 11.093352 [27,] 11.173228 [28,] 11.249266 [29,] 11.321770 [30,] 11.391015 [31,] 11.457246 [32,] 11.520684 [33,] 11.581530 [34,] 11.639964 [35,] 11.696151 [36,] 11.750241 [37,] 11.802370 [38,] 11.852663 [39,] 11.901235 [40,] 11.948191 [41,] 11.993627 [42,] 12.037633 [43,] 12.080292 [44,] 12.121678 [45,] 12.161864 [46,] 12.200913 [47,] 12.238888 [48,] 12.275844 [49,] 12.311835 [50,] 12.346908 [51,] 12.381111 [52,] 12.414486 [53,] 12.447072 [54,] 12.478909 [55,] 12.510030 [56,] 12.540469 [57,] 12.570257 [58,] 12.599424 [59,] 12.627996 [60,] 12.656000 [61,] 12.683461 [62,] 12.710400 [63,] 12.736841 [64,] 12.762804 [65,] 12.788308 [66,] 12.813372 [67,] 12.838013 [68,] 12.862248 [69,] 12.886092 [70,] 12.909561 [71,] 12.932669 [72,] 12.955429 [73,] 12.977854 [74,] 12.999957 [75,] 13.021748 [76,] 13.043240 [77,] 13.064442 [78,] 13.085365 [79,] 13.106018 [80,] 13.126410 [81,] 13.146551 [82,] 13.166448 [83,] 13.186110 [84,] 13.205545 [85,] 13.224758 [86,] 13.243759 [87,] 13.262553 [88,] 13.281146 [89,] 13.299546 [90,] 13.317758 [91,] 13.335787 [92,] 13.353639 [93,] 13.371319 [94,] 13.388833 [95,] 13.406185 [96,] 13.423379 [97,] 13.440421 [98,] 13.457314 [99,] 13.474063 [100,] 13.490672 [101,] 13.507144 [102,] 13.523483 [103,] 13.539693 [104,] 13.555777 [105,] 13.571738 [106,] 13.587580 [107,] 13.603305 [108,] 13.618918 [109,] 13.634419 [110,] 13.649813 [111,] 13.665102 [112,] 13.680288 [113,] 13.695374 [114,] 13.710362 [115,] 13.725255 [116,] 13.740054 [117,] 13.754763 [118,] 13.769383 [119,] 13.783916 [120,] 13.798364 [121,] 13.812729 [122,] 13.827013 [123,] 13.841218 [124,] 13.855345 [125,] 13.869396 [126,] 13.883373 [127,] 13.897276 [128,] 13.911109 [129,] 13.924871 [130,] 13.938566 [131,] 13.952193 [132,] 13.965754 [133,] 13.979252 [134,] 13.992686 [135,] 14.006058 [136,] 14.019370 [137,] 14.032622 [138,] 14.045815 [139,] 14.058952 [140,] 14.072032 [141,] 14.085057 [142,] 14.098027 [143,] 14.110945 [144,] 14.123810 [145,] 14.136623 [146,] 14.149387 [147,] 14.162100 [148,] 14.174765 [149,] 14.187382 [150,] 14.199952 [151,] 14.212475 [152,] 14.224953 [153,] 14.237386 [154,] 14.249775 [155,] 14.262120 [156,] 14.274423 [157,] 14.286683 [158,] 14.298902 [159,] 14.311080 [160,] 14.323218 [161,] 14.335317 [162,] 14.347376 [163,] 14.359397 [164,] 14.371380 [165,] 14.383325 [166,] 14.395234 [167,] 14.407107 [168,] 14.418943 [169,] 14.430745 [170,] 14.442511 [171,] 14.454243 [172,] 14.465941 [173,] 14.477606 [174,] 14.489237 [175,] 14.500836 [176,] 14.512403 [177,] 14.523938 [178,] 14.535442 [179,] 14.546914 [180,] 14.558356 [181,] 14.569767 [182,] 14.581149 [183,] 14.592501 [184,] 14.603824 [185,] 14.615118 [186,] 14.626383 [187,] 14.637620 [188,] 14.648829 [189,] 14.660011 [190,] 14.671165 [191,] 14.682292 [192,] 14.693393 [193,] 14.704466 [194,] 14.715514 [195,] 14.726536 [196,] 14.737532 [197,] 14.748503 [198,] 14.759448 [199,] 14.770369 [200,] 14.781265 [201,] 14.792137 [202,] 14.802984 [203,] 14.813807 [204,] 14.824607 [205,] 14.835383 [206,] 14.846136 [207,] 14.856865 [208,] 14.867572 [209,] 14.878256 [210,] 14.888918 [211,] 14.899557 [212,] 14.910174 [213,] 14.920769 [214,] 14.931343 [215,] 14.941894 [216,] 14.952425 [217,] 14.962934 [218,] 14.973422 [219,] 14.983889 [220,] 14.994336 [221,] 15.004761 [222,] 15.015167 [223,] 15.025552 [224,] 15.035917 [225,] 15.046262 [226,] 15.056587 [227,] 15.066892 [228,] 15.077178 [229,] 15.087445 [230,] 15.097692 [231,] 15.107920 [232,] 15.118129 [233,] 15.128319 [234,] 15.138490 [235,] 15.148643 [236,] 15.158777 [237,] 15.168892 [238,] 15.178990 [239,] 15.189069 [240,] 15.199130 [241,] 15.209173 [242,] 15.219198 [243,] 15.229206 [244,] 15.239196 [245,] 15.249168 [246,] 15.259123 [247,] 15.269060 [248,] 15.278981 [249,] 15.288884 [250,] 15.298770 [251,] 15.308639 [252,] 15.318491 [253,] 15.328326 [254,] 15.338145 [255,] 15.347947 [256,] 15.357733 [257,] 15.367502 [258,] 15.377255 [259,] 15.386991 [260,] 15.396712 [261,] 15.406416 [262,] 15.416104 [263,] 15.425776 [264,] 15.435433 [265,] 15.445073 [266,] 15.454698 [267,] 15.464307 [268,] 15.473901 [269,] 15.483479 [270,] 15.493041 [271,] 15.502589 [272,] 15.512121 [273,] 15.521637 [274,] 15.531139 [275,] 15.540625 [276,] 15.550096 [277,] 15.559553 [278,] 15.568994 [279,] 15.578420 [280,] 15.587832 [281,] 15.597229 [282,] 15.606611 [283,] 15.615979 [284,] 15.625332 [285,] 15.634670 [286,] 15.643994 [287,] 15.653304 [288,] 15.662599 [289,] 15.671880 [290,] 15.681147 [291,] 15.690400 [292,] 15.699638 [293,] 15.708862 [294,] 15.718073 [295,] 15.727269 [296,] 15.736451 [297,] 15.745619 [298,] 15.754774 [299,] 15.763915 [300,] 15.773042 [301,] 15.782155 [302,] 15.791255 [303,] 15.800341 [304,] 15.809413 [305,] 15.818472 [306,] 15.827517 [307,] 15.836549 [308,] 15.845568 [309,] 15.854573 [310,] 15.863565 [311,] 15.872543 [312,] 15.881509 [313,] 15.890461 [314,] 15.899400 [315,] 15.908326 [316,] 15.917238 [317,] 15.926138 [318,] 15.935025 [319,] 15.943899 [320,] 15.952759 [321,] 15.961607 [322,] 15.970442 [323,] 15.979265 [324,] 15.988074 [325,] 15.996871 [326,] 16.005655 [327,] 16.014426 [328,] 16.023185 [329,] 16.031931 [330,] 16.040665 [331,] 16.049386 [332,] 16.058094 [333,] 16.066790 [334,] 16.075474 [335,] 16.084145 [336,] 16.092803 [337,] 16.101450 [338,] 16.110084 [339,] 16.118706 [340,] 16.127315 [341,] 16.135913 [342,] 16.144498 [343,] 16.153071 [344,] 16.161632 [345,] 16.170180 [346,] 16.178717 [347,] 16.187242 [348,] 16.195754 [349,] 16.204255 [350,] 16.212743 [351,] 16.221220 [352,] 16.229685 [353,] 16.238138 [354,] 16.246579 [355,] 16.255008 [356,] 16.263426 [357,] 16.271832 [358,] 16.280226 [359,] 16.288608 [360,] 16.296979 [361,] 16.305338 [362,] 16.313685 [363,] 16.322021 [364,] 16.330345 [365,] 16.338657 [366,] 16.346959 [367,] 16.355248 [368,] 16.363526 [369,] 16.371793 [370,] 16.380048 [371,] 16.388292 [372,] 16.396524 [373,] 16.404745 [374,] 16.412955 [375,] 16.421154 [376,] 16.429341 [377,] 16.437517 [378,] 16.445681 [379,] 16.453835 [380,] 16.461977 [381,] 16.470108 [382,] 16.478228 [383,] 16.486337 [384,] 16.494435 [385,] 16.502521 [386,] 16.510597 [387,] 16.518662 [388,] 16.526715 [389,] 16.534758 [390,] 16.542790 [391,] 16.550810 [392,] 16.558820 [393,] 16.566819 [394,] 16.574807 [395,] 16.582784 [396,] 16.590751 [397,] 16.598706 [398,] 16.606651 [399,] 16.614585 [400,] 16.622508 [401,] 16.630421 > > if(doExtras) { + ## Extend the range even more -- change u0 + show.(m3 <- p.log.f(10, 500, d = 5, u0 = .96)) + ## and more + show.(m3.2 <- p.log.f(10, 800, d = 5, u0 = .96))#-> breakdown at ~ 775.. + } > showProc.time() Time elapsed: 0.25 0 0.25 > > > ## higher d: > if(doExtras) { + show.(m4 <- p.log.f(10, 500, d = 12, u0 = 0.95)) + m5 <- p.log.f(10, 500, d = 12, u0 = 0.92) + m6 <- p.log.f(1, 200, d = 12, u0 = c(0.8, 0.9), legend.x="topright") + m7 <- p.log.f(1, 400, d = 12, u0 = c(0.88, 0.9)) + } > > ## whereas this now also overflows for Eulerian *AND* asymp. + log.arg: > show.(mm <- p.log.f(10, 1000, d = 12, u0 = 0.9)) theta Stirling_Li_n(log(arg)) Eulerian_Li_n(log(arg)) [1,] 10.00000 17.94816 17.94816 [2,] 13.88235 19.73673 19.73673 [3,] 17.76471 21.35567 21.35567 [4,] 21.64706 22.84110 22.84110 [5,] 25.52941 24.21004 24.21004 [6,] 29.41176 25.47414 25.47414 [7,] 33.29412 26.64318 26.64318 [8,] 37.17647 27.72609 27.72609 [9,] 41.05882 28.73110 28.73110 [10,] 44.94118 29.66588 29.66588 [11,] 48.82353 30.53743 30.53743 [12,] 52.70588 31.35211 31.35211 [13,] 56.58824 32.11565 32.11565 [14,] 60.47059 32.83320 32.83320 [15,] 64.35294 33.50929 33.50929 [16,] 68.23529 34.14798 34.14798 [17,] 72.11765 34.75283 34.75283 [18,] 76.00000 Inf 35.32700 [19,] 79.88235 Inf 35.87327 [20,] 83.76471 Inf 36.39409 [21,] 87.64706 Inf 36.89164 [22,] 91.52941 Inf 37.36786 [23,] 95.41176 Inf 37.82444 [24,] 99.29412 Inf 38.26291 [25,] 103.17647 Inf 38.68464 [26,] 107.05882 Inf 39.09084 [27,] 110.94118 Inf 39.48260 [28,] 114.82353 Inf 39.86090 [29,] 118.70588 Inf 40.22664 [30,] 122.58824 Inf 40.58062 [31,] 126.47059 Inf 40.92357 [32,] 130.35294 Inf 41.25615 [33,] 134.23529 Inf 41.57898 [34,] 138.11765 Inf 41.89260 [35,] 142.00000 Inf 42.19753 [36,] 145.88235 Inf 42.49424 [37,] 149.76471 Inf 42.78315 [38,] 153.64706 Inf 43.06467 [39,] 157.52941 Inf 43.33916 [40,] 161.41176 Inf 43.60697 [41,] 165.29412 Inf 43.86842 [42,] 169.17647 Inf 44.12379 [43,] 173.05882 Inf 44.37338 [44,] 176.94118 Inf 44.61742 [45,] 180.82353 Inf 44.85617 [46,] 184.70588 Inf 45.08984 [47,] 188.58824 Inf 45.31865 [48,] 192.47059 Inf 45.54280 [49,] 196.35294 Inf 45.76248 [50,] 200.23529 Inf 45.97785 [51,] 204.11765 Inf 46.18909 [52,] 208.00000 Inf 46.39635 [53,] 211.88235 Inf 46.59977 [54,] 215.76471 Inf 46.79950 [55,] 219.64706 Inf 46.99567 [56,] 223.52941 Inf 47.18840 [57,] 227.41176 Inf 47.37781 [58,] 231.29412 Inf 47.56402 [59,] 235.17647 Inf 47.74713 [60,] 239.05882 Inf 47.92723 [61,] 242.94118 Inf 48.10444 [62,] 246.82353 Inf 48.27884 [63,] 250.70588 Inf 48.45051 [64,] 254.58824 Inf 48.61955 [65,] 258.47059 Inf 48.78603 [66,] 262.35294 Inf 48.95003 [67,] 266.23529 Inf 49.11161 [68,] 270.11765 Inf 49.27086 [69,] 274.00000 Inf 49.42784 [70,] 277.88235 Inf 49.58260 [71,] 281.76471 Inf 49.73522 [72,] 285.64706 Inf 49.88575 [73,] 289.52941 Inf 50.03425 [74,] 293.41176 Inf 50.18077 [75,] 297.29412 Inf 50.32537 [76,] 301.17647 Inf 50.46809 [77,] 305.05882 Inf 50.60898 [78,] 308.94118 Inf 50.74809 [79,] 312.82353 Inf 50.88546 [80,] 316.70588 Inf 51.02114 [81,] 320.58824 Inf 51.15516 [82,] 324.47059 Inf 51.28757 [83,] 328.35294 Inf 51.41841 [84,] 332.23529 Inf 51.54771 [85,] 336.11765 Inf 51.67550 [86,] 340.00000 Inf 51.80183 [87,] 343.88235 Inf 51.92672 [88,] 347.76471 Inf 52.05022 [89,] 351.64706 Inf 52.17234 [90,] 355.52941 Inf 52.29312 [91,] 359.41176 Inf 52.41258 [92,] 363.29412 Inf 52.53077 [93,] 367.17647 Inf 52.64770 [94,] 371.05882 Inf 52.76339 [95,] 374.94118 Inf 52.87789 [96,] 378.82353 Inf 52.99120 [97,] 382.70588 Inf 53.10336 [98,] 386.58824 Inf 53.21439 [99,] 390.47059 Inf 53.32431 [100,] 394.35294 Inf 53.43314 [101,] 398.23529 Inf 53.54090 [102,] 402.11765 Inf 53.64762 [103,] 406.00000 Inf 53.75331 [104,] 409.88235 Inf 53.85800 [105,] 413.76471 Inf 53.96170 [106,] 417.64706 Inf 54.06443 [107,] 421.52941 Inf 54.16621 [108,] 425.41176 Inf 54.26706 [109,] 429.29412 Inf 54.36699 [110,] 433.17647 Inf 54.46603 [111,] 437.05882 Inf 54.56417 [112,] 440.94118 Inf 54.66145 [113,] 444.82353 Inf 54.75788 [114,] 448.70588 Inf 54.85347 [115,] 452.58824 Inf 54.94824 [116,] 456.47059 Inf 55.04220 [117,] 460.35294 Inf 55.13536 [118,] 464.23529 Inf 55.22773 [119,] 468.11765 Inf 55.31934 [120,] 472.00000 Inf 55.41020 [121,] 475.88235 Inf 55.50031 [122,] 479.76471 Inf 55.58968 [123,] 483.64706 Inf 55.67834 [124,] 487.52941 Inf 55.76628 [125,] 491.41176 Inf 55.85353 [126,] 495.29412 Inf 55.94010 [127,] 499.17647 Inf 56.02598 [128,] 503.05882 Inf 56.11121 [129,] 506.94118 Inf 56.19577 [130,] 510.82353 Inf 56.27969 [131,] 514.70588 Inf 56.36298 [132,] 518.58824 Inf 56.44564 [133,] 522.47059 Inf 56.52768 [134,] 526.35294 Inf 56.60912 [135,] 530.23529 Inf 56.68996 [136,] 534.11765 Inf 56.77021 [137,] 538.00000 Inf 56.84987 [138,] 541.88235 Inf 56.92897 [139,] 545.76471 Inf 57.00750 [140,] 549.64706 Inf 57.08547 [141,] 553.52941 Inf 57.16289 [142,] 557.41176 Inf 57.23977 [143,] 561.29412 Inf 57.31612 [144,] 565.17647 Inf 57.39195 [145,] 569.05882 Inf 57.46725 [146,] 572.94118 Inf 57.54204 [147,] 576.82353 Inf 57.61633 [148,] 580.70588 Inf 57.69012 [149,] 584.58824 Inf 57.76341 [150,] 588.47059 Inf 57.83622 [151,] 592.35294 Inf 57.90856 [152,] 596.23529 Inf 57.98042 [153,] 600.11765 Inf 58.05181 [154,] 604.00000 Inf 58.12274 [155,] 607.88235 Inf 58.19322 [156,] 611.76471 Inf 58.26325 [157,] 615.64706 Inf 58.33284 [158,] 619.52941 Inf 58.40199 [159,] 623.41176 Inf 58.47071 [160,] 627.29412 Inf 58.53900 [161,] 631.17647 Inf 58.60687 [162,] 635.05882 Inf 58.67432 [163,] 638.94118 Inf 58.74136 [164,] 642.82353 Inf 58.80800 [165,] 646.70588 Inf 58.87424 [166,] 650.58824 Inf 58.94007 [167,] 654.47059 Inf 59.00552 [168,] 658.35294 Inf 59.07058 [169,] 662.23529 Inf 59.13526 [170,] 666.11765 Inf 59.19956 [171,] 670.00000 Inf 59.26348 [172,] 673.88235 Inf 59.32704 [173,] 677.76471 Inf 59.39023 [174,] 681.64706 Inf 59.45306 [175,] 685.52941 Inf 59.51553 [176,] 689.41176 Inf 59.57765 [177,] 693.29412 Inf 59.63943 [178,] 697.17647 Inf 59.70085 [179,] 701.05882 Inf 59.76194 [180,] 704.94118 Inf 59.82269 [181,] 708.82353 Inf 59.88310 [182,] 712.70588 Inf 59.94319 [183,] 716.58824 Inf 60.00294 [184,] 720.47059 Inf 60.06238 [185,] 724.35294 Inf 60.12149 [186,] 728.23529 Inf 60.18029 [187,] 732.11765 Inf 60.23878 [188,] 736.00000 Inf 60.29696 [189,] 739.88235 Inf 60.35483 [190,] 743.76471 Inf 60.41240 [191,] 747.64706 Inf 60.46967 [192,] 751.52941 Inf 60.52664 [193,] 755.41176 Inf 60.58332 [194,] 759.29412 Inf 60.63971 [195,] 763.17647 Inf 60.69581 [196,] 767.05882 Inf 60.75163 [197,] 770.94118 Inf 60.80716 [198,] 774.82353 Inf 60.86242 [199,] 778.70588 Inf 60.91740 [200,] 782.58824 Inf 60.97210 [201,] 786.47059 Inf 61.02654 [202,] 790.35294 Inf 61.08070 [203,] 794.23529 Inf 61.13461 [204,] 798.11765 Inf 61.18824 [205,] 802.00000 Inf 61.24162 [206,] 805.88235 Inf 61.29474 [207,] 809.76471 Inf 61.34761 [208,] 813.64706 Inf 61.40024 [209,] 817.52941 Inf 61.45329 [210,] 821.41176 Inf 61.49269 [211,] 825.29412 Inf 62.34770 [212,] 829.17647 NaN Inf [213,] 833.05882 NaN Inf [214,] 836.94118 NaN Inf [215,] 840.82353 NaN Inf [216,] 844.70588 NaN Inf [217,] 848.58824 NaN Inf [218,] 852.47059 NaN Inf [219,] 856.35294 NaN Inf [220,] 860.23529 NaN Inf [221,] 864.11765 NaN Inf [222,] 868.00000 NaN Inf [223,] 871.88235 NaN Inf [224,] 875.76471 NaN Inf [225,] 879.64706 NaN Inf [226,] 883.52941 NaN Inf [227,] 887.41176 NaN Inf [228,] 891.29412 NaN Inf [229,] 895.17647 NaN Inf [230,] 899.05882 NaN Inf [231,] 902.94118 NaN Inf [232,] 906.82353 NaN Inf [233,] 910.70588 NaN Inf [234,] 914.58824 NaN Inf [235,] 918.47059 NaN Inf [236,] 922.35294 NaN Inf [237,] 926.23529 NaN Inf [238,] 930.11765 NaN Inf [239,] 934.00000 NaN Inf [240,] 937.88235 NaN Inf [241,] 941.76471 NaN Inf [242,] 945.64706 NaN Inf [243,] 949.52941 NaN Inf [244,] 953.41176 NaN Inf [245,] 957.29412 NaN Inf [246,] 961.17647 NaN Inf [247,] 965.05882 NaN Inf [248,] 968.94118 NaN Inf [249,] 972.82353 NaN Inf [250,] 976.70588 NaN Inf [251,] 980.58824 NaN Inf [252,] 984.47059 NaN Inf [253,] 988.35294 NaN Inf [254,] 992.23529 NaN Inf [255,] 996.11765 NaN Inf [256,] 1000.00000 NaN Inf asymp-w_Li_n(log(arg)) [1,] 17.94816 [2,] 19.73673 [3,] 21.35567 [4,] 22.84110 [5,] 24.21004 [6,] 25.47414 [7,] 26.64318 [8,] 27.72609 [9,] 28.73110 [10,] 29.66588 [11,] 30.53743 [12,] 31.35211 [13,] 32.11565 [14,] 32.83320 [15,] 33.50929 [16,] 34.14798 [17,] 34.75283 [18,] 35.32700 [19,] 35.87327 [20,] 36.39409 [21,] 36.89164 [22,] 37.36786 [23,] 37.82444 [24,] 38.26291 [25,] 38.68464 [26,] 39.09084 [27,] 39.48260 [28,] 39.86090 [29,] 40.22664 [30,] 40.58062 [31,] 40.92357 [32,] 41.25615 [33,] 41.57898 [34,] 41.89260 [35,] 42.19753 [36,] 42.49424 [37,] 42.78315 [38,] 43.06467 [39,] 43.33916 [40,] 43.60697 [41,] 43.86842 [42,] 44.12379 [43,] 44.37338 [44,] 44.61742 [45,] 44.85617 [46,] 45.08984 [47,] 45.31865 [48,] 45.54280 [49,] 45.76248 [50,] 45.97785 [51,] 46.18909 [52,] 46.39635 [53,] 46.59977 [54,] 46.79950 [55,] 46.99567 [56,] 47.18840 [57,] 47.37781 [58,] 47.56402 [59,] 47.74713 [60,] 47.92723 [61,] 48.10444 [62,] 48.27884 [63,] 48.45051 [64,] 48.61955 [65,] 48.78603 [66,] 48.95003 [67,] 49.11161 [68,] 49.27086 [69,] 49.42784 [70,] 49.58260 [71,] 49.73522 [72,] 49.88575 [73,] 50.03425 [74,] 50.18077 [75,] 50.32537 [76,] 50.46809 [77,] 50.60898 [78,] 50.74809 [79,] 50.88546 [80,] 51.02114 [81,] 51.15516 [82,] 51.28757 [83,] 51.41841 [84,] 51.54771 [85,] 51.67550 [86,] 51.80183 [87,] 51.92672 [88,] 52.05022 [89,] 52.17234 [90,] 52.29312 [91,] 52.41258 [92,] 52.53077 [93,] 52.64770 [94,] 52.76339 [95,] 52.87789 [96,] 52.99120 [97,] 53.10336 [98,] 53.21439 [99,] 53.32431 [100,] 53.43314 [101,] 53.54090 [102,] 53.64762 [103,] 53.75331 [104,] 53.85800 [105,] 53.96170 [106,] 54.06443 [107,] 54.16621 [108,] 54.26706 [109,] 54.36699 [110,] 54.46603 [111,] 54.56417 [112,] 54.66145 [113,] 54.75788 [114,] 54.85347 [115,] 54.94824 [116,] 55.04220 [117,] 55.13536 [118,] 55.22773 [119,] 55.31934 [120,] 55.41020 [121,] 55.50031 [122,] 55.58968 [123,] 55.67834 [124,] 55.76628 [125,] 55.85353 [126,] 55.94010 [127,] 56.02598 [128,] 56.11121 [129,] 56.19577 [130,] 56.27969 [131,] 56.36298 [132,] 56.44564 [133,] 56.52768 [134,] 56.60912 [135,] 56.68996 [136,] 56.77021 [137,] 56.84987 [138,] 56.92897 [139,] 57.00750 [140,] 57.08547 [141,] 57.16289 [142,] 57.23977 [143,] 57.31612 [144,] 57.39195 [145,] 57.46725 [146,] 57.54204 [147,] 57.61633 [148,] 57.69012 [149,] 57.76341 [150,] 57.83622 [151,] 57.90856 [152,] 57.98042 [153,] 58.05181 [154,] 58.12274 [155,] 58.19322 [156,] 58.26325 [157,] 58.33284 [158,] 58.40199 [159,] 58.47071 [160,] 58.53900 [161,] 58.60687 [162,] 58.67432 [163,] 58.74136 [164,] 58.80800 [165,] 58.87424 [166,] 58.94007 [167,] 59.00552 [168,] 59.07058 [169,] 59.13526 [170,] 59.19956 [171,] 59.26348 [172,] 59.32704 [173,] 59.39023 [174,] 59.45306 [175,] 59.51553 [176,] 59.57765 [177,] 59.63943 [178,] 59.70085 [179,] 59.76194 [180,] 59.82269 [181,] 59.88310 [182,] 59.94319 [183,] 60.00294 [184,] 60.06238 [185,] 60.12149 [186,] 60.18029 [187,] 60.23878 [188,] 60.29696 [189,] 60.35483 [190,] 60.41240 [191,] 60.46967 [192,] 60.52664 [193,] 60.58332 [194,] 60.63971 [195,] 60.69581 [196,] 60.75163 [197,] 60.80716 [198,] 60.86242 [199,] 60.91740 [200,] 60.97210 [201,] 61.02654 [202,] 61.08070 [203,] 61.13461 [204,] 61.18824 [205,] 61.24162 [206,] 61.29474 [207,] 61.34761 [208,] 61.40024 [209,] 61.45329 [210,] 61.49269 [211,] 62.34770 [212,] Inf [213,] Inf [214,] Inf [215,] Inf [216,] Inf [217,] Inf [218,] Inf [219,] Inf [220,] Inf [221,] Inf [222,] Inf [223,] Inf [224,] Inf [225,] Inf [226,] Inf [227,] Inf [228,] Inf [229,] Inf [230,] Inf [231,] Inf [232,] Inf [233,] Inf [234,] Inf [235,] Inf [236,] Inf [237,] Inf [238,] Inf [239,] Inf [240,] Inf [241,] Inf [242,] Inf [243,] Inf [244,] Inf [245,] Inf [246,] Inf [247,] Inf [248,] Inf [249,] Inf [250,] Inf [251,] Inf [252,] Inf [253,] Inf [254,] Inf [255,] Inf [256,] Inf There were 45 warnings (use warnings() to see them) > > ##--> investigation shows that w is *also* underflowing to 0 (!!!) : > > myTracer <- quote({ + cat(sprintf("Tracing %s: variables at entry are\n", + ## the -4 is purely by trial and error -- better? + deparse(sys.call(-4)[1]))) + print(ls.str()) + }) > trace(polylog, tracer = myTracer, print=FALSE) Tracing function "polylog" in package "copula" [1] "polylog" > ## Tracing function "polylog" in package "copula" > > f <- attr(mm,"dacopula") > > f(attr(mm,"u"), 825, log=TRUE, method="negI-s-asymp") Tracing polylog(): variables at entry are asymp.w.order : num 0 is.log.z : logi TRUE is.logmlog : logi FALSE logarithm : logi TRUE method : chr "negI-s-asymp" n.sum : s : num -11 z : num -4.15e-322 [1] 61.48259 > ## .... > ## z : num -4.15e-322 > ## [1] 61.48259 > f(attr(mm,"u"), 800, log=TRUE, method="negI-s-asymp") Tracing polylog(): variables at entry are asymp.w.order : num 0 is.log.z : logi TRUE is.logmlog : logi FALSE logarithm : logi TRUE method : chr "negI-s-asymp" n.sum : s : num -11 z : num -2.44e-312 [1] 61.21416 > ## z : num -2.44e-312 > > if(FALSE) + debug(f) > ## and then realize that in f(), i.e., Frank's dacopula() final line, > ## (d-1)*log(theta) + Li. - theta*u.sum - lu > ## the terms Li. - (theta * u.sum) > ## 'more or less cancel' -- so if we could compute a > ## *different* rescaled polylog() .. maybe we can get better > > untrace(polylog) Untracing function "polylog" in package "copula" > > if(doExtras) + ## and similarly here: + ll <- p.log.f(1, 12000, d = 12, u0 = 0.08) > > showProc.time() Time elapsed: 0.39 0 0.4 > > ###------- Diagonal MLE --- dDiag() , etc ---------------------------------- > > ## Some basic dDiag() checks {that failed earlier}: > stopifnot(identical(0, dDiag(0, onacopulaL("AMH", list(0.5, 1:4)))), + all.equal(dDiag(.9, onacopulaL("Frank", list(44, 1:3))), + 1.008252, tolerance= 1e-5), + all.equal(dDiag(.9, onacopulaL("Joe", list(44, 1:3))), + 1.025283, tolerance= 1e-5), + TRUE) > > demo("dDiag-plots", package = "copula") demo(dDiag-plots) ---- ~~~~~~~~~~~ > ## Copyright (C) 2012 Marius Hofert, Ivan Kojadinovic, Martin Maechler, and Jun Yan > ## > ## This program is free software; you can redistribute it and/or modify it under > ## the terms of the GNU General Public License as published by the Free Software > ## Foundation; either version 3 of the License, or (at your option) any later > ## version. > ## > ## This program is distributed in the hope that it will be useful, but WITHOUT > ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS > ## FOR A PARTICULAR PURPOSE. See the GNU General Public License for more > ## details. > ## > ## You should have received a copy of the GNU General Public License along with > ## this program; if not, see . > > > p.dDiag <- + function(family, d, n = 1000, log = FALSE, + tau = c(0.01, + ## First tau = 0.01 for all families + if(family=="AMH") c((1:3)/10, 0.33) else + c((1:9)/10, .95, 0.99, .999)), + lty = 1:3, + cols = adjustcolor(colorRampPalette(c("red", "orange", "blue"), + space = "Lab")(length(tau)), 0.8)) + { + stopifnot(length(d) == 1, d == as.integer(d), d >= 2, n >= 10) + cop <- getAcop(family) + th. <- cop@iTau(tau) + u <- seq(0,1, length = n) + mainTit <- + paste("Diagonal densities of ", family, "\n d = ",d, + if(log)", log = TRUE", sep="") + c.tau <- format(tau) + c.th <- format(th., scientific=FALSE, digits = 4) + + ## yMat <- sapply(th., function(theta) + ## dDiag(u, cop=onacopulaL(family, list(theta, 1:d)), + ## log=log)) + ## The "generic" Archimedean dDiag(): + yM.a <- sapply(th., function(theta) + copula:::dDiagA(u, d=d, cop= setTheta(cop, theta), log=log)) + ## The dDiag-*s*lot : + yM.s <- sapply(th., function(theta) + cop@dDiag(u, theta=theta, d=d, log=log)) + + p.oneMat <- function(yMat, label) { + thisTit <- paste(label, mainTit, sep=": ") + non.fin <- !is.finite(yMat) + ## Rather: in log-scale, -Inf may be ok : + if(log) non.fin <- non.fin & (is.na(yMat) | yMat != -Inf) + nF.mat <- + if(any(has.nF <- apply(non.fin, 2, any))) { + i.non.fin <- apply(non.fin, 2, which.max) + i.non.fin[!has.nF] <- NA + cat(thisTit,":\n non-finite values -- for (theta, u >= *):\n") + print(cbind(theta = th.[has.nF], u = u[i.non.fin[has.nF]])) + } + matplot(u, yMat, type = "l", col = cols, lty = lty, + ylab="dDiag(u, *)", main= thisTit) + if(any(has.nF)) { ## mark the end points + i1 <- i.non.fin - 1 + points(u[i1], yMat[cbind(i1, seq_along(tau))], + col = cols, cex = 1.5, lwd=2, pch = 4) + } + abline(h=0, lty=3) + lleg <- lapply(seq_along(tau), function(j) { + cc <- if(has.nF[j]) { + i <- i.non.fin[j] + sprintf(" f(%.3g) = %g", u[i], yMat[i,j]) + } + substitute(list(tau == TAU, theta == TH) ~ COMM, + list(TAU= c.tau[j], TH= c.th[j], COMM= cc)) + + }) + legend(if(log)"bottomright" else "topleft", + do.call(expression, lleg), + lty = lty, col=cols, bty="n") + }## end{ p.oneMat() } + + p.oneMat(yM.a, "dDiagA()") + p.oneMat(yM.s, "cop @ dDiag()") + + invisible(list(d = d, tau=tau, theta=th., u = u, dDiag.a = yM.a, dDiag.s = yM.s)) + } > if(!exists("doExtras") || !is.logical(doExtras)) + doExtras <- interactive() || nzchar(Sys.getenv("R_copula_check_extra")) > doExtras [1] FALSE > par(ask = dev.interactive(orNone = TRUE)) > r.dDiag.3 <- lapply(.ac.longNames, + function(family) p.dDiag(family, d = 3)) dDiagA(): Diagonal densities of Clayton d = 3 : non-finite values -- for (theta, u >= *): theta u [1,] 198 0.001001001 [2,] 1998 0.001001001 dDiagA(): Diagonal densities of Frank d = 3 : non-finite values -- for (theta, u >= *): theta u [1,] 3998.354 0.1781782 dDiagA(): Diagonal densities of Gumbel d = 3 : non-finite values -- for (theta, u >= *): theta u [1,] 1.010101 1.000000000 [2,] 1.111111 1.000000000 [3,] 1.250000 1.000000000 [4,] 1.428571 1.000000000 [5,] 1.666667 1.000000000 [6,] 2.000000 1.000000000 [7,] 2.500000 1.000000000 [8,] 3.333333 1.000000000 [9,] 5.000000 1.000000000 [10,] 10.000000 1.000000000 [11,] 20.000000 1.000000000 [12,] 100.000000 1.000000000 [13,] 1000.000000 0.001001001 dDiagA(): Diagonal densities of Joe d = 3 : non-finite values -- for (theta, u >= *): theta u [1,] 1.017448 1.0000000 [2,] 1.194410 1.0000000 [3,] 1.443824 1.0000000 [4,] 1.772108 1.0000000 [5,] 2.219066 1.0000000 [6,] 2.856238 1.0000000 [7,] 3.826639 1.0000000 [8,] 5.463756 1.0000000 [9,] 8.767709 1.0000000 [10,] 18.738688 1.0000000 [11,] 38.724328 1.0000000 [12,] 198.712959 0.9739740 [13,] 1998.710414 0.3023023 > r.dDiag.4.L <- lapply(.ac.longNames, + function(family) p.dDiag(family, d = 4, log=TRUE)) dDiagA(): Diagonal densities of Frank d = 4, log = TRUE : non-finite values -- for (theta, u >= *): theta u [1,] 3998.354 0.1871872 dDiagA(): Diagonal densities of Gumbel d = 4, log = TRUE : non-finite values -- for (theta, u >= *): theta u [1,] 1.010101 1.0000000 [2,] 1.111111 1.0000000 [3,] 1.250000 1.0000000 [4,] 1.428571 1.0000000 [5,] 1.666667 1.0000000 [6,] 2.000000 1.0000000 [7,] 2.500000 1.0000000 [8,] 3.333333 1.0000000 [9,] 5.000000 1.0000000 [10,] 10.000000 1.0000000 [11,] 20.000000 1.0000000 [12,] 100.000000 1.0000000 [13,] 1000.000000 0.6226226 dDiagA(): Diagonal densities of Joe d = 4, log = TRUE : non-finite values -- for (theta, u >= *): theta u [1,] 1.017448 1.0000000 [2,] 1.194410 1.0000000 [3,] 1.443824 1.0000000 [4,] 1.772108 1.0000000 [5,] 2.219066 1.0000000 [6,] 2.856238 1.0000000 [7,] 3.826639 1.0000000 [8,] 5.463756 1.0000000 [9,] 8.767709 1.0000000 [10,] 18.738688 1.0000000 [11,] 38.724328 1.0000000 [12,] 198.712959 0.9769770 [13,] 1998.710414 0.3113113 > r.dDiag.15 <- lapply(.ac.longNames, + function(family) p.dDiag(family, d = 15)) dDiagA(): Diagonal densities of Clayton d = 15 : non-finite values -- for (theta, u >= *): theta u [1,] 198 0.001001001 [2,] 1998 0.001001001 dDiagA(): Diagonal densities of Frank d = 15 : non-finite values -- for (theta, u >= *): theta u [1,] 3998.354 0.1791792 dDiagA(): Diagonal densities of Gumbel d = 15 : non-finite values -- for (theta, u >= *): theta u [1,] 1.010101 1.000000000 [2,] 1.111111 1.000000000 [3,] 1.250000 1.000000000 [4,] 1.428571 1.000000000 [5,] 1.666667 1.000000000 [6,] 2.000000 1.000000000 [7,] 2.500000 1.000000000 [8,] 3.333333 1.000000000 [9,] 5.000000 1.000000000 [10,] 10.000000 1.000000000 [11,] 20.000000 1.000000000 [12,] 100.000000 1.000000000 [13,] 1000.000000 0.001001001 dDiagA(): Diagonal densities of Joe d = 15 : non-finite values -- for (theta, u >= *): theta u [1,] 1.017448 1.0000000 [2,] 1.194410 1.0000000 [3,] 1.443824 1.0000000 [4,] 1.772108 1.0000000 [5,] 2.219066 1.0000000 [6,] 2.856238 1.0000000 [7,] 3.826639 1.0000000 [8,] 5.463756 1.0000000 [9,] 8.767709 1.0000000 [10,] 18.738688 1.0000000 [11,] 38.724328 1.0000000 [12,] 198.712959 0.9739740 [13,] 1998.710414 0.3023023 > r.dDiag.75 <- lapply(.ac.longNames, + function(family) p.dDiag(family, d = 75)) dDiagA(): Diagonal densities of Clayton d = 75 : non-finite values -- for (theta, u >= *): theta u [1,] 198 0.001001001 [2,] 1998 0.001001001 dDiagA(): Diagonal densities of Frank d = 75 : non-finite values -- for (theta, u >= *): theta u [1,] 3998.354 0.1791792 dDiagA(): Diagonal densities of Gumbel d = 75 : non-finite values -- for (theta, u >= *): theta u [1,] 1.010101 1.000000000 [2,] 1.111111 1.000000000 [3,] 1.250000 1.000000000 [4,] 1.428571 1.000000000 [5,] 1.666667 1.000000000 [6,] 2.000000 1.000000000 [7,] 2.500000 1.000000000 [8,] 3.333333 1.000000000 [9,] 5.000000 1.000000000 [10,] 10.000000 1.000000000 [11,] 20.000000 1.000000000 [12,] 100.000000 1.000000000 [13,] 1000.000000 0.001001001 dDiagA(): Diagonal densities of Joe d = 75 : non-finite values -- for (theta, u >= *): theta u [1,] 1.017448 1.0000000 [2,] 1.194410 1.0000000 [3,] 1.443824 1.0000000 [4,] 1.772108 1.0000000 [5,] 2.219066 1.0000000 [6,] 2.856238 1.0000000 [7,] 3.826639 1.0000000 [8,] 5.463756 1.0000000 [9,] 8.767709 1.0000000 [10,] 18.738688 1.0000000 [11,] 38.724328 1.0000000 [12,] 198.712959 0.9739740 [13,] 1998.710414 0.3023023 > r.dDiag.200.L <- lapply(.ac.longNames, + function(family) p.dDiag(family, d = 200, log=TRUE)) dDiagA(): Diagonal densities of Frank d = 200, log = TRUE : non-finite values -- for (theta, u >= *): theta u [1,] 3998.354 0.1871872 dDiagA(): Diagonal densities of Gumbel d = 200, log = TRUE : non-finite values -- for (theta, u >= *): theta u [1,] 1.010101 1.0000000 [2,] 1.111111 1.0000000 [3,] 1.250000 1.0000000 [4,] 1.428571 1.0000000 [5,] 1.666667 1.0000000 [6,] 2.000000 1.0000000 [7,] 2.500000 1.0000000 [8,] 3.333333 1.0000000 [9,] 5.000000 1.0000000 [10,] 10.000000 1.0000000 [11,] 20.000000 1.0000000 [12,] 100.000000 1.0000000 [13,] 1000.000000 0.6226226 dDiagA(): Diagonal densities of Joe d = 200, log = TRUE : non-finite values -- for (theta, u >= *): theta u [1,] 1.017448 1.0000000 [2,] 1.194410 1.0000000 [3,] 1.443824 1.0000000 [4,] 1.772108 1.0000000 [5,] 2.219066 1.0000000 [6,] 2.856238 1.0000000 [7,] 3.826639 1.0000000 [8,] 5.463756 1.0000000 [9,] 8.767709 1.0000000 [10,] 18.738688 1.0000000 [11,] 38.724328 1.0000000 [12,] 198.712959 0.9769770 [13,] 1998.710414 0.3113113 > if(doExtras) + demo("dDiag-plots-part-2") > ## ------------------ ../demo/dDiag-plots-part-2.R > ## FIXME? get 'echo', 'ask' , 'verbose' from *our own* caller -- how? > ## -----------> ../demo/dDiag-plots.R > > showProc.time() Time elapsed: 1.97 0.08 2.04 > > ## A problem to solve: > curve(copAMH@dDiag(exp(-x), theta=.9, d=4, log=TRUE), 0, 500, n=2001, + col = 2, lwd= 2, + main = "FIXME (*not* urgent): dDiag(u, , log=TRUE) for small u") > > showProc.time() Time elapsed: 0.02 0 0.02 > > ##----------- Test all estimation methods ---------------- > > ## all methods > (estMeth <- eval(formals(enacopula)$method)) [1] "mle" "smle" "dmle" "mde.chisq.CvM" [5] "mde.chisq.KS" "mde.gamma.CvM" "mde.gamma.KS" "tau.tau.mean" [9] "tau.theta.mean" "beta" > > ## Use selected (the best 6) estimation methods: > ## (estMeth <- c("mle", "smle", "dmle", "tau.tau.mean", "mde.chisq.CvM", "mde.chisq.KS")) > > n <- if(doExtras) 128 else 8 > d <- 6 ; theta <- 2 > (cop <- onacopula("Clayton", C(theta, 1:d))) Nested Archimedean copula ("outer_nacopula" of dim. 6), with slot 'comp' = (1, 2, 3, 4, 5, 6) and root 'copula' = Archimedean copula ("acopula"), family "Clayton", theta= (2) and *no* child copulas > set.seed(1) > U <- rnacopula(n, cop) > ## Show that 'racopula()' is indeed identical {and faster for the case of non-nested} > stopifnot(identical(U, + {set.seed(1); copula:::racopula(n, copClayton, theta=theta, d=d)})) > > rr <- sapply(estMeth, function(e) { + enacopula(U, cop, e, n.MC = if(e == "smle") 1000 else 0) + }) > > round(cbind(rr, bias = rr - theta), 3) rr bias mle 2.997 0.997 smle 2.922 0.922 dmle 1.284 -0.716 mde.chisq.CvM 2.846 0.846 mde.chisq.KS 2.725 0.725 mde.gamma.CvM 2.710 0.710 mde.gamma.KS 0.000 -2.000 tau.tau.mean 1.750 -0.250 tau.theta.mean 2.086 0.086 beta 0.843 -1.157 > > showProc.time() Time elapsed: 0.06 0.02 0.08 > > ### 2-level nested Archimedean copulas : > if(doExtras) { + ## was simply (!) : demo("dnac-demo") + v.dnac <- vignette("dNAC", package="copula") + R.dnac <- + if(nzchar(v.dnac$R)) + with(v.dnac, file.path(Dir, "doc", R)) + else { ## default: knitr::rmarkdown() unfortunately does *not* provide R + stopifnot(file.exists(ff <- with(v.dnac, file.path(Dir, "doc", File)))) + oR <- file.path(tempdir(), sub("Rmd$", "R", basename(ff))) + knitr::purl(ff, output = oR) + oR + } + cat("R file: ", R.dnac, "\n -- now source()ing it :\n") + source(R.dnac, echo = TRUE, max.deparse.length = Inf, + keep.source = TRUE, encoding = getOption("encoding")) + } > > > showProc.time() Time elapsed: 0 0 0 > > proc.time() user system elapsed 4.26 0.25 4.50