diff --git a/NEWS.md b/NEWS.md index a52ecdf..10ba254 100644 --- a/NEWS.md +++ b/NEWS.md @@ -3,6 +3,12 @@ * Add `math_lcm()` and `math_gcd()` to compute the lowest common multiple and the greatest common divisor. * Add `interval_hdr()` and `interval_credible()` to compute the credible intervals. +## Bugfixes & changes +* `jackknife()` gained a new argument to apply a function on the leave-one-out values (`f`). + +## Internals +* Add `with_seed()` to evaluate an expression with a temporarily seed. + # arkhe 1.1.0 ## New classes and methods * Add `needs()` to check for the availability of a package. diff --git a/R/AllGenerics.R b/R/AllGenerics.R index 942992b..c89a853 100644 --- a/R/AllGenerics.R +++ b/R/AllGenerics.R @@ -494,17 +494,17 @@ setGeneric( #' `do`) as argument. #' @param ... Extra arguments to be passed to `do`. #' @return -#' If `f` is not `NULL`, `bootstrap()` returns the result of `f` applied to -#' the `n` values of `do`. -#' -#' If `f` is `NULL`, `bootstrap()` returns a named `numeric` `vector` with the -#' following elements: +#' If `f` is `NULL` (the default), `bootstrap()` returns a named `numeric` +#' vector with the following elements: #' \describe{ #' \item{`original`}{The observed value of `do` applied to `object`.} #' \item{`mean`}{The bootstrap estimate of mean of `do`.} #' \item{`bias`}{The bootstrap estimate of bias of `do`.} #' \item{`error`}{he bootstrap estimate of standard error of `do`.} #' } +#' +#' If `f` is a `function`, `bootstrap()` returns the result of `f` applied to +#' the `n` values of `do`. #' @example inst/examples/ex-resample.R #' @author N. Frerebeau #' @docType methods @@ -523,14 +523,20 @@ setGeneric( #' @param do A [`function`] that takes `object` as an argument and returns a #' single numeric value. #' @param ... Extra arguments to be passed to `do`. +#' @param f A [`function`] that takes a single numeric vector (the leave-one-out +#' values of `do`) as argument. #' @return -#' Returns a named `numeric` `vector` with the following elements: +#' If `f` is `NULL` (the default), `jackknife()` returns a named `numeric` +#' vector with the following elements: #' \describe{ #' \item{`original`}{The observed value of `do` applied to `object`.} #' \item{`mean`}{The jackknife estimate of mean of `do`.} #' \item{`bias`}{The jackknife estimate of bias of `do`.} #' \item{`error`}{he jackknife estimate of standard error of `do`.} #' } +#' +#' If `f` is a `function`, `jackknife()` returns the result of `f` applied to +#' the leave-one-out values of `do`. #' @example inst/examples/ex-resample.R #' @author N. Frerebeau #' @docType methods diff --git a/R/statistics.R b/R/statistics.R index 4072916..d9b2e02 100644 --- a/R/statistics.R +++ b/R/statistics.R @@ -164,8 +164,9 @@ setMethod( spl <- sample(object, size = length(object) * n, replace = TRUE) replicates <- t(matrix(spl, nrow = n)) values <- apply(X = replicates, MARGIN = 2, FUN = do, ...) - values <- if (is.function(f)) f(values) else summary_bootstrap(values, hat) - values + + if (is.function(f)) return(f(values)) + summary_bootstrap(values, hat) } ) @@ -187,7 +188,7 @@ summary_bootstrap <- function(x, hat) { setMethod( f = "jackknife", signature = c(object = "numeric"), - definition = function(object, do, ...) { + definition = function(object, do, ..., f = NULL) { n <- length(object) hat <- do(object, ...) @@ -199,8 +200,9 @@ setMethod( FUN.VALUE = double(1), x = object, do = do, ... ) - values <- summary_jackknife(values, hat) - values + + if (is.function(f)) return(f(values)) + summary_jackknife(values, hat) } ) diff --git a/R/utilities.R b/R/utilities.R index 53c9b66..5b3a403 100644 --- a/R/utilities.R +++ b/R/utilities.R @@ -1,5 +1,43 @@ # HELPERS +## https://stackoverflow.com/questions/56191862/where-do-i-specify-random-seed-for-tests-in-r-package +#' Evaluate an Expression with a Temporarily Seed +#' +#' @param expr An [`expression`] to be evaluated. +#' @param seed A single value to be passed to [set.seed()]. +#' @param envir The [environment][environment()] in which `expr` should be +#' evaluated. +#' @param rounding A [`logical`] scalar: should the default discrete uniform +#' generation method in \R versions prior to 3.6.0 be used? Usefull for unit +#' testing. +#' @param ... Further arguments to be passed to [set.seed()]. +#' @return +#' The results of `expr` evaluated. +#' @seealso [set.seed()] +#' @keywords internal +with_seed <- function(expr, seed, ..., envir = parent.frame(), rounding = TRUE) { + expr <- substitute(expr) + ## Save and restore the random number generator (RNG) state + env <- globalenv() + old_seed <- env$.Random.seed + on.exit({ + if (is.null(old_seed)) { + rm(list = ".Random.seed", envir = env, inherits = FALSE) + } else { + assign(".Random.seed", value = old_seed, envir = env, inherits = FALSE) + } + }) + ## Keep the results the same for R versions prior to 3.6 + if (isTRUE(rounding) && getRversion() >= "3.6") { + ## Set sample.kind = "Rounding" to reproduce the old sampling + ## Suppress warning "non-uniform 'Rounding' sampler used" + suppressWarnings(set.seed(seed, sample.kind = "Rounding")) + } else { + set.seed(seed) + } + eval(expr, envir = envir) +} + # Helpers ====================================================================== #' Helpers #' diff --git a/inst/examples/ex-resample.R b/inst/examples/ex-resample.R index fdf32b3..b2b71b2 100644 --- a/inst/examples/ex-resample.R +++ b/inst/examples/ex-resample.R @@ -7,5 +7,11 @@ bootstrap(x, do = mean, n = 100) quant <- function(x) { quantile(x, probs = c(0.25, 0.75)) } bootstrap(x, n = 100, do = mean, f = quant) +## Get the n bootstrap values +bootstrap(x, n = 100, do = mean, f = function(x) { x }) + ## Jackknife jackknife(x, do = mean) # Sample mean + +## Get the leave-one-out values instead of summary +jackknife(x, do = mean, f = function(x) { x }) diff --git a/man/bootstrap.Rd b/man/bootstrap.Rd index 33774dc..a938e35 100644 --- a/man/bootstrap.Rd +++ b/man/bootstrap.Rd @@ -26,17 +26,17 @@ replications.} \code{do}) as argument.} } \value{ -If \code{f} is not \code{NULL}, \code{bootstrap()} returns the result of \code{f} applied to -the \code{n} values of \code{do}. - -If \code{f} is \code{NULL}, \code{bootstrap()} returns a named \code{numeric} \code{vector} with the -following elements: +If \code{f} is \code{NULL} (the default), \code{bootstrap()} returns a named \code{numeric} +vector with the following elements: \describe{ \item{\code{original}}{The observed value of \code{do} applied to \code{object}.} \item{\code{mean}}{The bootstrap estimate of mean of \code{do}.} \item{\code{bias}}{The bootstrap estimate of bias of \code{do}.} \item{\code{error}}{he bootstrap estimate of standard error of \code{do}.} } + +If \code{f} is a \code{function}, \code{bootstrap()} returns the result of \code{f} applied to +the \code{n} values of \code{do}. } \description{ Samples randomly from the elements of \code{object} with replacement. @@ -51,8 +51,14 @@ bootstrap(x, do = mean, n = 100) quant <- function(x) { quantile(x, probs = c(0.25, 0.75)) } bootstrap(x, n = 100, do = mean, f = quant) +## Get the n bootstrap values +bootstrap(x, n = 100, do = mean, f = function(x) { x }) + ## Jackknife jackknife(x, do = mean) # Sample mean + +## Get the leave-one-out values instead of summary +jackknife(x, do = mean, f = function(x) { x }) } \seealso{ Other resampling methods: diff --git a/man/jackknife.Rd b/man/jackknife.Rd index 9ac4bdd..d6dd287 100644 --- a/man/jackknife.Rd +++ b/man/jackknife.Rd @@ -9,7 +9,7 @@ \usage{ jackknife(object, ...) -\S4method{jackknife}{numeric}(object, do, ...) +\S4method{jackknife}{numeric}(object, do, ..., f = NULL) } \arguments{ \item{object}{A \code{\link{numeric}} vector.} @@ -18,15 +18,22 @@ jackknife(object, ...) \item{do}{A \code{\link{function}} that takes \code{object} as an argument and returns a single numeric value.} + +\item{f}{A \code{\link{function}} that takes a single numeric vector (the leave-one-out +values of \code{do}) as argument.} } \value{ -Returns a named \code{numeric} \code{vector} with the following elements: +If \code{f} is \code{NULL} (the default), \code{jackknife()} returns a named \code{numeric} +vector with the following elements: \describe{ \item{\code{original}}{The observed value of \code{do} applied to \code{object}.} \item{\code{mean}}{The jackknife estimate of mean of \code{do}.} \item{\code{bias}}{The jackknife estimate of bias of \code{do}.} \item{\code{error}}{he jackknife estimate of standard error of \code{do}.} } + +If \code{f} is a \code{function}, \code{jackknife()} returns the result of \code{f} applied to +the leave-one-out values of \code{do}. } \description{ Jackknife Estimation @@ -41,8 +48,14 @@ bootstrap(x, do = mean, n = 100) quant <- function(x) { quantile(x, probs = c(0.25, 0.75)) } bootstrap(x, n = 100, do = mean, f = quant) +## Get the n bootstrap values +bootstrap(x, n = 100, do = mean, f = function(x) { x }) + ## Jackknife jackknife(x, do = mean) # Sample mean + +## Get the leave-one-out values instead of summary +jackknife(x, do = mean, f = function(x) { x }) } \seealso{ Other resampling methods: diff --git a/man/with_seed.Rd b/man/with_seed.Rd new file mode 100644 index 0000000..fb32c91 --- /dev/null +++ b/man/with_seed.Rd @@ -0,0 +1,32 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/utilities.R +\name{with_seed} +\alias{with_seed} +\title{Evaluate an Expression with a Temporarily Seed} +\usage{ +with_seed(expr, seed, ..., envir = parent.frame(), rounding = TRUE) +} +\arguments{ +\item{expr}{An \code{\link{expression}} to be evaluated.} + +\item{seed}{A single value to be passed to \code{\link[=set.seed]{set.seed()}}.} + +\item{...}{Further arguments to be passed to \code{\link[=set.seed]{set.seed()}}.} + +\item{envir}{The \link[=environment]{environment} in which \code{expr} should be +evaluated.} + +\item{rounding}{A \code{\link{logical}} scalar: should the default discrete uniform +generation method in \R versions prior to 3.6.0 be used? Usefull for unit +testing.} +} +\value{ +The results of \code{expr} evaluated. +} +\description{ +Evaluate an Expression with a Temporarily Seed +} +\seealso{ +\code{\link[=set.seed]{set.seed()}} +} +\keyword{internal} diff --git a/tests/testthat/_snaps/statistics.md b/tests/testthat/_snaps/statistics.md index e9d14ec..06efd4f 100644 --- a/tests/testthat/_snaps/statistics.md +++ b/tests/testthat/_snaps/statistics.md @@ -61,3 +61,57 @@ [3,] 0.06413649 0.1558635 [4,] 0.27661853 0.4133815 +# Bootstrap + + Code + bootstrap_summary + Output + original mean bias error + 0.07651681 0.05439040 -0.02212641 0.19517006 + +--- + + Code + bootstrap_values + Output + [1] 0.274803536 0.271761759 0.131918811 -0.108350186 -0.407722040 + [6] 0.474370950 0.096520144 0.140542686 0.048469812 0.045673600 + [11] 0.098911690 0.174660229 0.407774218 0.076149892 0.233883643 + [16] -0.158214771 0.296334314 -0.236746678 -0.027278798 0.140680055 + [21] 0.084317514 0.131248858 -0.126159227 0.351657449 0.130019851 + [26] -0.064196657 0.049937514 -0.166241173 -0.085116371 0.003131881 + [31] 0.159010387 0.126015372 0.175198395 0.127748767 -0.047995525 + [36] 0.245223426 -0.029654335 0.153690734 0.372854685 -0.061229774 + [41] 0.102818533 0.037641950 0.125330950 -0.082734254 -0.081537996 + [46] 0.003646969 0.096332439 0.116993888 -0.176113259 0.345299218 + [51] 0.595379927 -0.311485051 0.093283378 0.342125925 0.001140435 + [56] 0.246825067 -0.173706701 0.113309462 0.338231445 -0.425915621 + [61] 0.446095061 0.156535443 -0.266725129 0.175778042 0.201231074 + [66] -0.065400705 0.015916353 0.257033280 -0.111331877 0.073891696 + [71] -0.004623160 -0.326442795 0.048737537 0.005317883 -0.053255377 + [76] 0.081287636 -0.013228710 0.344409439 -0.125705295 -0.181941858 + [81] 0.118077490 -0.154583474 0.035114092 0.026580327 -0.119802041 + [86] -0.197818974 0.107870115 0.253120777 -0.199949361 0.050689044 + [91] 0.002231045 0.309152956 -0.390878169 0.027390498 0.026426616 + [96] 0.223192535 -0.022840641 -0.146577129 -0.014807361 0.038401655 + +# Jackknife + + Code + jackknife(x, do = mean) + Output + original mean bias error + 0.07651681 0.07651681 0.00000000 0.18647363 + +--- + + Code + jackknife(x, do = mean, f = function(x) { + x + }) + Output + [1] 0.04972670 0.04320369 0.08629682 0.10441228 0.04865520 0.17622590 + [7] 0.04738093 0.09508001 0.09549979 0.12892938 0.08666232 -0.01510399 + [13] 0.06103728 0.05316420 0.12004569 0.03754928 0.12719441 0.09799546 + [19] 0.02155913 0.06482171 + diff --git a/tests/testthat/helper.R b/tests/testthat/helper.R deleted file mode 100644 index 506528a..0000000 --- a/tests/testthat/helper.R +++ /dev/null @@ -1,21 +0,0 @@ -## Helpers for testthat - -## Set the seed for just a code block -## https://stackoverflow.com/questions/56191862/where-do-i-specify-random-seed-for-tests-in-r-package -with_seed <- function(seed, code) { - code <- substitute(code) - ## Save and restore the random number generator (RNG) state - if (exists(".Random.seed", mode = "numeric")) { - orig.seed <- .Random.seed - on.exit(.Random.seed <<- orig.seed) - } - ## Keep the results the same for R versions prior to 3.6 - if (getRversion() >= "3.6") { - ## Set sample.kind = "Rounding" to reproduce the old sampling - ## Suppress warning "non-uniform 'Rounding' sampler used" - suppressWarnings(set.seed(seed, sample.kind = "Rounding")) - } else { - set.seed(seed) - } - eval.parent(code) -} diff --git a/tests/testthat/test-statistics.R b/tests/testthat/test-statistics.R index d1cd2ec..00b45d1 100644 --- a/tests/testthat/test-statistics.R +++ b/tests/testthat/test-statistics.R @@ -21,3 +21,19 @@ test_that("Confidence interval for binomial proportions", { expect_snapshot(confidence_multinomial(x)) expect_snapshot(confidence_multinomial(x, corrected = TRUE)) }) +test_that("Bootstrap", { + bootstrap_summary <- with_seed({ + bootstrap(rnorm(20), n = 100, do = mean) + }, seed = 12345) + expect_snapshot(bootstrap_summary) + + bootstrap_values <- with_seed({ + bootstrap(rnorm(20), n = 100, do = mean, f = function(x) { x }) + }, seed = 12345) + expect_snapshot(bootstrap_values) +}) +test_that("Jackknife", { + x <- with_seed(rnorm(20), seed = 12345) + expect_snapshot(jackknife(x, do = mean)) + expect_snapshot(jackknife(x, do = mean, f = function(x) { x })) +})