From 954d4f66266cb7d53eddb6be7b9821638220c678 Mon Sep 17 00:00:00 2001 From: Matthew Kay Date: Mon, 27 Dec 2021 21:22:43 -0600 Subject: [PATCH 1/6] fix minor doc issues --- R/hilo.R | 1 - R/plot.R | 2 +- man/autoplot.distribution.Rd | 1 + 3 files changed, 2 insertions(+), 2 deletions(-) diff --git a/R/hilo.R b/R/hilo.R index 61add7dd..43fd4f61 100644 --- a/R/hilo.R +++ b/R/hilo.R @@ -109,7 +109,6 @@ vec_math.hilo <- function(.fn, .x, ...){ vec_restore(out, .x) } -#' @rdname vctrs-compat #' @method vec_arith hilo #' @export vec_arith.hilo <- function(op, x, y, ...){ diff --git a/R/plot.R b/R/plot.R index b91b9932..448be04b 100644 --- a/R/plot.R +++ b/R/plot.R @@ -10,7 +10,7 @@ #' @param x The distribution(s) to plot. #' @param ... Unused. #' -#' @keyword internal +#' @keywords internal #' #' @export autoplot.distribution <- function(x, ...){ diff --git a/man/autoplot.distribution.Rd b/man/autoplot.distribution.Rd index 1121f295..5cf2cdc8 100644 --- a/man/autoplot.distribution.Rd +++ b/man/autoplot.distribution.Rd @@ -20,3 +20,4 @@ the {ggdist} package to produce your own distribution plots. You can learn more about how this plot can be produced using {ggdist} here: https://mjskay.github.io/ggdist/articles/slabinterval.html } +\keyword{internal} From 77dace998d04ef2501dcb65615b8468a7eeccf1e Mon Sep 17 00:00:00 2001 From: Matthew Kay Date: Mon, 27 Dec 2021 22:15:49 -0600 Subject: [PATCH 2/6] basic support for inverses of some known functions --- R/default.R | 29 ++++++++++++++++++++++++++++- 1 file changed, 28 insertions(+), 1 deletion(-) diff --git a/R/default.R b/R/default.R index 60ab3126..05ee442d 100644 --- a/R/default.R +++ b/R/default.R @@ -159,12 +159,39 @@ dim.dist_default <- function(x){ invert_fail <- function(...) stop("Inverting transformations for distributions is not yet supported.") +inverse_functions <- exprs( + sqrt = function(x) x^2, + exp = log, + log = function(x, base = exp(1)) base ^ x, + log2 = function(x) 2^x, + log10 = function(x) 10^x, + expm1 = log1p, + log1p = expm1, + cos = acos, + sin = asin, + tan = atan, + acos = cos, + asin = sin, + atan = tan, + cosh = acosh, + sinh = asinh, + tanh = atanh, + acosh = cosh, + asinh = sinh, + atanh = tanh +) + #' @method Math dist_default #' @export Math.dist_default <- function(x, ...) { if(dim(x) > 1) stop("Transformations of multivariate distributions are not yet supported.") + trans <- new_function(exprs(x = ), body = expr((!!sym(.Generic))(x, !!!dots_list(...)))) - vec_data(dist_transformed(wrap_dist(list(x)), trans, invert_fail))[[1]] + + inverse_fun <- inverse_functions[[.Generic]] %||% invert_fail + inverse <- new_function(exprs(x = ), body = expr((!!inverse_fun)(x, !!!dots_list(...)))) + + vec_data(dist_transformed(wrap_dist(list(x)), trans, inverse))[[1]] } #' @method Ops dist_default From e63165be474116420ee3b21a976861ada4dbed25 Mon Sep 17 00:00:00 2001 From: Matthew Kay Date: Mon, 27 Dec 2021 22:16:30 -0600 Subject: [PATCH 3/6] support difference bases for log(dist_lognormal()) --- R/dist_lognormal.R | 17 ++++++++++++++--- 1 file changed, 14 insertions(+), 3 deletions(-) diff --git a/R/dist_lognormal.R b/R/dist_lognormal.R index 2d86fab9..142414fe 100644 --- a/R/dist_lognormal.R +++ b/R/dist_lognormal.R @@ -145,10 +145,21 @@ kurtosis.dist_lognormal <- function(x, ...) { exp(4*s2) + 2*exp(3*s2) + 3*exp(2*s2) - 6 } +# make a normal distribution from a lognormal distribution using the +# specified base +normal_dist_with_base <- function(x, base = exp(1)) { + vec_data(dist_normal(x[["mu"]], x[["sigma"]]) / log(base))[[1]] +} + #' @method Math dist_lognormal #' @export Math.dist_lognormal <- function(x, ...) { - # Shortcut to get Normal distribution from log-normal. - if(.Generic == "log") return(vec_data(dist_normal(x[["mu"]], x[["sigma"]]))[[1]]) - NextMethod() + switch(.Generic, + # Shortcuts to get Normal distribution from log-normal. + log = normal_dist_with_base(x, ...), + log2 = normal_dist_with_base(x, 2), + log10 = normal_dist_with_base(x, 10), + + NextMethod() + ) } From a41b710aaadf1cbdd76a4fd669b0575961c81ffb Mon Sep 17 00:00:00 2001 From: Matthew Kay Date: Mon, 27 Dec 2021 22:36:56 -0600 Subject: [PATCH 4/6] support inverse functions in Math.dist_transformed --- R/default.R | 10 ++++++++-- R/transformed.R | 6 +++++- 2 files changed, 13 insertions(+), 3 deletions(-) diff --git a/R/default.R b/R/default.R index 05ee442d..c7496db4 100644 --- a/R/default.R +++ b/R/default.R @@ -158,7 +158,6 @@ dim.dist_default <- function(x){ } invert_fail <- function(...) stop("Inverting transformations for distributions is not yet supported.") - inverse_functions <- exprs( sqrt = function(x) x^2, exp = log, @@ -180,6 +179,13 @@ inverse_functions <- exprs( asinh = sinh, atanh = tanh ) +#' Attempt to get the inverse of a known function by name. Returns invert_fail +#' (a function that raises an error if called) if there is no known inverse. +#' @param f string. Name of a function. +#' @noRd +get_inverse_function <- function(f) { + inverse_functions[[f]] %||% invert_fail +} #' @method Math dist_default #' @export @@ -188,7 +194,7 @@ Math.dist_default <- function(x, ...) { trans <- new_function(exprs(x = ), body = expr((!!sym(.Generic))(x, !!!dots_list(...)))) - inverse_fun <- inverse_functions[[.Generic]] %||% invert_fail + inverse_fun <- get_inverse_function(.Generic) inverse <- new_function(exprs(x = ), body = expr((!!inverse_fun)(x, !!!dots_list(...)))) vec_data(dist_transformed(wrap_dist(list(x)), trans, inverse))[[1]] diff --git a/R/transformed.R b/R/transformed.R index 0ce7ca06..12517faa 100644 --- a/R/transformed.R +++ b/R/transformed.R @@ -83,7 +83,11 @@ covariance.dist_transformed <- function(x, ...){ #' @export Math.dist_transformed <- function(x, ...) { trans <- new_function(exprs(x = ), body = expr((!!sym(.Generic))((!!x$transform)(x), !!!dots_list(...)))) - vec_data(dist_transformed(wrap_dist(list(x[["dist"]])), trans, invert_fail))[[1]] + + inverse_fun <- get_inverse_function(.Generic) + inverse <- new_function(exprs(x = ), body = expr((!!x$inverse)((!!inverse_fun)(x, !!!dots_list(...))))) + + vec_data(dist_transformed(wrap_dist(list(x[["dist"]])), trans, inverse))[[1]] } #' @method Ops dist_transformed From fbbbd648f07e1f0e1a91ca3feb401f3ae7e2ad3e Mon Sep 17 00:00:00 2001 From: Matthew Kay Date: Tue, 28 Dec 2021 00:07:08 -0600 Subject: [PATCH 5/6] tests for inverting transformations of distributions --- R/default.R | 104 +++++++++++++++++++------- R/transformed.R | 15 +++- tests/testthat/test-transformations.R | 43 +++++++++++ 3 files changed, 132 insertions(+), 30 deletions(-) mode change 100644 => 100755 R/default.R mode change 100644 => 100755 R/transformed.R diff --git a/R/default.R b/R/default.R old mode 100644 new mode 100755 index c7496db4..4a4ba481 --- a/R/default.R +++ b/R/default.R @@ -158,33 +158,73 @@ dim.dist_default <- function(x){ } invert_fail <- function(...) stop("Inverting transformations for distributions is not yet supported.") -inverse_functions <- exprs( - sqrt = function(x) x^2, - exp = log, - log = function(x, base = exp(1)) base ^ x, - log2 = function(x) 2^x, - log10 = function(x) 10^x, - expm1 = log1p, - log1p = expm1, - cos = acos, - sin = asin, - tan = atan, - acos = cos, - asin = sin, - atan = tan, - cosh = acosh, - sinh = asinh, - tanh = atanh, - acosh = cosh, - asinh = sinh, - atanh = tanh -) -#' Attempt to get the inverse of a known function by name. Returns invert_fail + +#' Attempt to get the inverse of f(x) by name. Returns invert_fail +#' (a function that raises an error if called) if there is no known inverse. +#' @param f string. Name of a function. +#' @noRd +get_unary_inverse <- function(f) { + switch(f, + sqrt = function(x) x^2, + exp = log, + log = function(x, base = exp(1)) base ^ x, + log2 = function(x) 2^x, + log10 = function(x) 10^x, + expm1 = log1p, + log1p = expm1, + cos = acos, + sin = asin, + tan = atan, + acos = cos, + asin = sin, + atan = tan, + cosh = acosh, + sinh = asinh, + tanh = atanh, + acosh = cosh, + asinh = sinh, + atanh = tanh, + + invert_fail + ) +} + +#' Attempt to get the inverse of f(x, constant) by name. Returns invert_fail #' (a function that raises an error if called) if there is no known inverse. #' @param f string. Name of a function. +#' @param constant a constant value #' @noRd -get_inverse_function <- function(f) { - inverse_functions[[f]] %||% invert_fail +get_binary_inverse_1 <- function(f, constant) { + force(constant) + + switch(f, + `+` = function(x) x - constant, + `-` = function(x) x + constant, + `*` = function(x) x / constant, + `/` = function(x) x * constant, + `^` = function(x) x ^ (1/constant), + + invert_fail + ) +} + +#' Attempt to get the inverse of f(constant, x) by name. Returns invert_fail +#' (a function that raises an error if called) if there is no known inverse. +#' @param f string. Name of a function. +#' @param constant a constant value +#' @noRd +get_binary_inverse_2 <- function(f, constant) { + force(constant) + + switch(f, + `+` = function(x) x - constant, + `-` = function(x) constant - x, + `*` = function(x) x / constant, + `/` = function(x) x * constant, + `^` = function(x) log(x, base = constant), + + invert_fail + ) } #' @method Math dist_default @@ -194,7 +234,7 @@ Math.dist_default <- function(x, ...) { trans <- new_function(exprs(x = ), body = expr((!!sym(.Generic))(x, !!!dots_list(...)))) - inverse_fun <- get_inverse_function(.Generic) + inverse_fun <- get_unary_inverse(.Generic) inverse <- new_function(exprs(x = ), body = expr((!!inverse_fun)(x, !!!dots_list(...)))) vec_data(dist_transformed(wrap_dist(list(x)), trans, inverse))[[1]] @@ -215,10 +255,18 @@ Ops.dist_default <- function(e1, e2) { stop(sprintf("The %s operation is not supported for <%s> and <%s>", .Generic, class(e1)[1], class(e2)[1])) } } else if(is_dist[1]){ - new_function(exprs(x = ), body = expr((!!sym(.Generic))((!!e1$transform)(x), !!e2))) + new_function(exprs(x = ), body = expr((!!sym(.Generic))(x, !!e2))) + } else { + new_function(exprs(x = ), body = expr((!!sym(.Generic))(!!e1, x))) + } + + inverse <- if(all(is_dist)) { + invert_fail + } else if(is_dist[1]){ + get_binary_inverse_1(.Generic, e2) } else { - new_function(exprs(x = ), body = expr((!!sym(.Generic))(!!e1, (!!e2$transform)(x)))) + get_binary_inverse_2(.Generic, e1) } - vec_data(dist_transformed(wrap_dist(list(e1,e2)[which(is_dist)]), trans, invert_fail))[[1]] + vec_data(dist_transformed(wrap_dist(list(e1,e2)[which(is_dist)]), trans, inverse))[[1]] } diff --git a/R/transformed.R b/R/transformed.R old mode 100644 new mode 100755 index 12517faa..f8eb7c59 --- a/R/transformed.R +++ b/R/transformed.R @@ -84,7 +84,7 @@ covariance.dist_transformed <- function(x, ...){ Math.dist_transformed <- function(x, ...) { trans <- new_function(exprs(x = ), body = expr((!!sym(.Generic))((!!x$transform)(x), !!!dots_list(...)))) - inverse_fun <- get_inverse_function(.Generic) + inverse_fun <- get_unary_inverse(.Generic) inverse <- new_function(exprs(x = ), body = expr((!!x$inverse)((!!inverse_fun)(x, !!!dots_list(...))))) vec_data(dist_transformed(wrap_dist(list(x[["dist"]])), trans, inverse))[[1]] @@ -105,5 +105,16 @@ Ops.dist_transformed <- function(e1, e2) { } else { new_function(exprs(x = ), body = expr((!!sym(.Generic))(!!e1, (!!e2$transform)(x)))) } - vec_data(dist_transformed(wrap_dist(list(list(e1,e2)[[which(is_dist)[1]]][["dist"]])), trans, invert_fail))[[1]] + + inverse <- if(all(is_dist)) { + invert_fail + } else if(is_dist[1]){ + inverse_fun <- get_binary_inverse_1(.Generic, e2) + new_function(exprs(x = ), body = expr((!!e1$inverse)((!!inverse_fun)(x)))) + } else { + inverse_fun <- get_binary_inverse_2(.Generic, e1) + new_function(exprs(x = ), body = expr((!!e2$inverse)((!!inverse_fun)(x)))) + } + + vec_data(dist_transformed(wrap_dist(list(list(e1,e2)[[which(is_dist)[1]]][["dist"]])), trans, inverse))[[1]] } diff --git a/tests/testthat/test-transformations.R b/tests/testthat/test-transformations.R index 0eb72147..a5b25555 100644 --- a/tests/testthat/test-transformations.R +++ b/tests/testthat/test-transformations.R @@ -43,6 +43,11 @@ test_that("LogNormal distributions", { dist_normal(0, 0.5) ) + # Test log() shortcut with different bases + expect_equal(log(dist_lognormal(0, log(3)), base = 3), dist_normal(0, 1)) + expect_equal(log2(dist_lognormal(0, log(2))), dist_normal(0, 1)) + expect_equal(log10(dist_lognormal(0, log(10))), dist_normal(0, 1)) + # format expect_equal(format(dist), sprintf("t(%s)", format(dist_normal(0, 0.5)))) @@ -71,3 +76,41 @@ test_that("LogNormal distributions", { expect_equal(mean(dist), exp(0.25/2), tolerance = 0.01) expect_equal(variance(dist), (exp(0.25) - 1)*exp(0.25), tolerance = 0.1) }) + +test_that("inverses are applied automatically", { + dist <- dist_gamma(1,1) + log2dist <- log(dist, base = 2) + log2dist_t <- dist_transformed(dist, log2, function(x) 2 ^ x) + + expect_equal(density(log2dist, 0.5), density(log2dist_t, 0.5)) + expect_equal(cdf(log2dist, 0.5), cdf(log2dist_t, 0.5)) + expect_equal(quantile(log2dist, 0.5), quantile(log2dist_t, 0.5)) + + # test multiple transformations that get stacked together by dist_transformed + explogdist <- exp(log(dist)) + expect_equal(density(dist, 0.5), density(explogdist, 0.5)) + expect_equal(cdf(dist, 0.5), cdf(explogdist, 0.5)) + expect_equal(quantile(dist, 0.5), quantile(explogdist, 0.5)) + + # test multiple transformations created by operators (via Ops) + explog2dist <- 2 ^ log2dist + expect_equal(density(dist, 0.5), density(explog2dist, 0.5)) + expect_equal(cdf(dist, 0.5), cdf(explog2dist, 0.5)) + expect_equal(quantile(dist, 0.5), quantile(explog2dist, 0.5)) + + # basic set of inverses + expect_equal(density(sqrt(dist^2), 0.5), density(dist, 0.5)) + expect_equal(density(exp(log(dist)), 0.5), density(dist, 0.5)) + expect_equal(density(10^(log10(dist)), 0.5), density(dist, 0.5)) + expect_equal(density(expm1(log1p(dist)), 0.5), density(dist, 0.5)) + expect_equal(density(cos(acos(dist)), 0.5), density(dist, 0.5)) + expect_equal(density(sin(asin(dist)), 0.5), density(dist, 0.5)) + expect_equal(density(tan(atan(dist)), 0.5), density(dist, 0.5)) + expect_equal(density(cosh(acosh(dist + 1)) - 1, 0.5), density(dist, 0.5)) + expect_equal(density(sinh(asinh(dist)), 0.5), density(dist, 0.5)) + expect_equal(density(tanh(atanh(dist)), 0.5), density(dist, 0.5)) + + expect_equal(density(dist + 1 - 1, 0.5), density(dist, 0.5)) + expect_equal(density(dist * 2 / 2, 0.5), density(dist, 0.5)) + +}) From cf3a2e54ee32589cd073cd28ec02b096a3002044 Mon Sep 17 00:00:00 2001 From: Matthew Kay Date: Tue, 28 Dec 2021 13:26:07 -0600 Subject: [PATCH 6/6] minor fixes for inverses and Jacobian correction to density --- R/default.R | 2 +- R/transformed.R | 2 +- tests/testthat/test-transformations.R | 4 ++++ 3 files changed, 6 insertions(+), 2 deletions(-) diff --git a/R/default.R b/R/default.R index 4a4ba481..1014ecd6 100755 --- a/R/default.R +++ b/R/default.R @@ -220,7 +220,7 @@ get_binary_inverse_2 <- function(f, constant) { `+` = function(x) x - constant, `-` = function(x) constant - x, `*` = function(x) x / constant, - `/` = function(x) x * constant, + `/` = function(x) constant / x, `^` = function(x) log(x, base = constant), invert_fail diff --git a/R/transformed.R b/R/transformed.R index f8eb7c59..fbd06ea8 100755 --- a/R/transformed.R +++ b/R/transformed.R @@ -38,7 +38,7 @@ format.dist_transformed <- function(x, ...){ #' @export density.dist_transformed <- function(x, at, ...){ - density(x[["dist"]], x[["inverse"]](at))*vapply(at, numDeriv::jacobian, numeric(1L), func = x[["inverse"]]) + density(x[["dist"]], x[["inverse"]](at))*abs(vapply(at, numDeriv::jacobian, numeric(1L), func = x[["inverse"]])) } #' @export diff --git a/tests/testthat/test-transformations.R b/tests/testthat/test-transformations.R index a5b25555..cf830e2e 100644 --- a/tests/testthat/test-transformations.R +++ b/tests/testthat/test-transformations.R @@ -113,4 +113,8 @@ test_that("inverses are applied automatically", { expect_equal(density(dist + 1 - 1, 0.5), density(dist, 0.5)) expect_equal(density(dist * 2 / 2, 0.5), density(dist, 0.5)) + # inverting a gamma distribution + expect_equal(density(1/dist_gamma(4, 3), 0.5), density(dist_inverse_gamma(4, 1/3), 0.5)) + expect_equal(density(1/(1/dist_gamma(4, 3)), 0.5), density(dist_gamma(4, 3), 0.5)) + })