https://github.com/cran/nacopula
Tip revision: 161411bb86f97e5a8bd89091cd61d03a33c2761a authored by Martin Maechler on 06 February 2012, 00:00:00 UTC
version 0.8-0
version 0.8-0
Tip revision: 161411b
retstable-ex.R
## Copyright (C) 2010-2012 Marius Hofert and Martin Maechler
##
## 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 <http://www.gnu.org/licenses/>.
#### Generating Exponentially Tilted Stable Random Vars (r ET stable)
#### ================================================================
#### Experiments with retstable*() versions
###
### More (computer intensive) experiments are in
### ../demo/retstable.R
### ~~~~~~~~~~~~~~~~~~~
require(nacopula)
source(system.file("Rsource", "utils.R", package="nacopula"))
##--> 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
## 983 times "MH", 17 x "LD"
set.seed(17); showSys.time(rET.H <- retstable (alpha, V0, method= "MH"))
set.seed(17); showSys.time(rET.D <- retstable (alpha, V0, method= "LD"))
set.seed(17); showSys.time(rET.R <- retstableR(alpha, V0))
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,])
## "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