R Under development (unstable) (2023-12-04 r85659 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. > ## 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 . > > > ### Generating Exponentially Tilted Stable Random Vars (r ET stable) > ### ================================================================ > ### Experiments with retstable*() versions > ### > ### More (computer intensive) experiments are in > ### ../demo/retstable.R > ### ~~~~~~~~~~~~~~~~~~~ > > require(copula) Loading required package: copula > source(system.file("Rsource", "utils.R", package="copula")) Loading required package: tools > ##--> tryCatch.W.E(), canGet() > > ## it works for 0-length V0 as well: > .N <- numeric(0) ; stopifnot(identical(.N, retstable(1/4, .N))) > > ## This is from "next version of Matrix" test-tools-1.R: > showSys.time <- function(expr) { + ## prepend 'Time' for R CMD Rdiff + st <- system.time(expr) + writeLines(paste("Time", capture.output(print(st)))) + invisible(st) + } > > ### using both retstableR() and retstable() > set.seed(1) > alpha <- .2 > V0 <- rgamma(2^12, 1) > set.seed(17); showSys.time(rET <- retstable (alpha, V0)) ## method = default: here takes Time user system elapsed Time 0.02 0.00 0.01 > ## 983 times "MH", 17 x "LD" > set.seed(17); showSys.time(rET.H <- retstable (alpha, V0, method= "MH")) Time user system elapsed Time 0 0 0 > set.seed(17); showSys.time(rET.D <- retstable (alpha, V0, method= "LD")) Time user system elapsed Time 0.02 0.00 0.02 > set.seed(17); showSys.time(rET.R <- retstableR(alpha, V0)) Time user system elapsed Time 0.27 0.00 0.26 > T <- function(r) r^(1/8) # log() is too much > bp <- boxplot(T(rET), T(rET.H), T(rET.D), T(rET.R), + notch=TRUE, col = "thistle") > (meds <- bp$stats[3,]) [1] 0.6268907 0.6131906 0.6164955 0.6215219 > > ## "H0": The 4 groups are not different -- here for the medians: > stopifnot(bp$conf[1,] < meds & meds < bp$conf[2,], + bp$stats > 0, + abs(bp$stats[2,] - 0.4035) < 0.007, ## first Quartiles + abs(bp$stats[4,] - 0.8085) < 0.006) ## 3rd Quartiles > > proc.time() user system elapsed 1.78 0.20 1.92