add more extensive tests

[venpopov]

While we figure out what the best approach for #101 is, I thought it might be good to extract some of the tests I wrote for it in a separate PR that could be reviewed more easily and implemented. I was thinking about this, because of your comment here. I hadn't realized that the implementation of nested distributions broke some functions such as quantile, generate, mean and covariance, because they use x[["transform"]]. This is because all the tests passed, but there are no tests for the behavior of these functions for nested transformations.

While trying to fix those for the PR, I realized that mean.transformed and covariance.transformed are alrady broken (#102), and I noted a problem with quantile a few days ago (#100). It would be good to have tests for those as well, to make sure that changes like those introduced by #101 or variations of it work with all of the package functionality. So here's my proposal:

Add methods eval_transform(dist, at), eval_inverse(dist, at) and eval_deriv(dist, at)

In the current master branch, these will be just wrappers for x[["transform"]], x[['inverse']] and numDeriv... But in the future, any changes to how the transform, inverse and deriv functions are stored will have to change these methods. This way the change needs to happen only in one place, whereas all functions like density.transformed, quantile.transformed, etc, while be the same. This should make code maintenance easier, and allow for more robust tests as well

add the tests from #101 that inverses and derivatives are calculated directly

there are some tests in the master branch, but again just for simple transformations, and would be good to have the extensive tests I wrote there

add more tests for quantile, generate, mean and covariance

the only test for mean and covariance of transformed distribution is for a longnormal with mean 0 and sd=0.25, and there the results are close enough, but as I noted in #102, the results are quite incorrect for other parameter values. I know these are supposed to be based on numeric approximation, but the error is quite substantial. There also no tests for nested transforms.

Let me know what you think and I can put together a PR. I would really love to be able to use the package with more confidence and I think this will also help catch more problems during development.

[mitchelloharawild]

Yes, evidently the tests for transformed distributions need to be fleshed out. Especially transformations which nest other transformations. It's also worth noting that the transformations should be monotonic over the supported domain (as per the documentation), but there aren't any checks for this (#7).

I'm not sure what you mean with the eval_() family proposed here, but I'm inclined to avoid tests which use the same code as the package itself - otherwise it is effectively testing nothing. Instead I think it is better to compare the package results with ground truth, and if there is slight inaccuracy due to numeric differentiation we can incorporate tolerance in the tests.

And yes, additional tests for other methods, especially at the boundaries/extremes where symbolic differentiation produces more accurate/exact results should be added to prevent future degradation.

[venpopov]

for the eval_() functions I meant something like this:

# Functions for evaluating the transformation of a random variable at a vector of values
eval_transform <- function(dist, at) {
  UseMethod("eval_transform")
}

#' @export
eval_transform.distribution <- function(dist, at) {
  at <- arg_listable(at, .ptype = NULL)
  dist_apply(dist, eval_transform, at = at)
}

#' @export
eval_transform.dist_default <- function(dist, at) {
  at
}

#' @export
eval_transform.dist_transformed <- function(dist, at) {
  dist$transform(at)
}

# Functions for evaluating the inverse transformation of a random variable at a vector of values
eval_inverse <- function(dist, at) {
  UseMethod("eval_inverse")
}

#' @export
eval_inverse.distribution <- function(dist, at) {
  at <- arg_listable(at, .ptype = NULL)
  dist_apply(dist, eval_inverse, at = at)
}

#' @export
eval_inverse.dist_default <- function(dist, at) {
  at
}

#' @export
eval_inverse.dist_transformed <- function(dist, at) {
  dist$inverse(at)
}

# Functions for evaluating the derivative of the inverse transformation of a random variable at a vector of values
eval_deriv <- function(dist, at) {
  UseMethod("eval_deriv")
}

#' @export
eval_deriv.distribution <- function(dist, at) {
  at <- arg_listable(at, .ptype = NULL)
  dist_apply(dist, eval_deriv, at = at)
}

#' @export
eval_deriv.dist_default <- function(dist, at) {
  at
}

#' @export
eval_deriv.dist_transformed <- function(dist, at) {
  inv <- function(at) suppressWarnings(eval_inverse(dist, at))
  vapply(at, numDeriv::jacobian, numeric(1L), func = inv)
}

The point is that this is a bit more convenient to use than vec_data(dist)[[1]]$transform(), and can be used to define tests that shouldn't change even if in the future you decide to reimplement how the inverses, transformations or derivatives are stored/accessed. E.g.

library(distributional)
library(testthat)
d <- dist_uniform(0, 1)
d2 <- -log(-log(d))

# can use the eval functions to evaluate the transformation, inverse, and derivative of the distribution at a given value (instead of vec_data(d2)[[1]]$transform)
at <- seq(0, 1, 0.1)
eval_transform(d2, at)
#> [[1]]
#>  [1]        -Inf -0.83403245 -0.47588500 -0.18562676  0.08742157  0.36651292
#>  [7]  0.67172699  1.03093043  1.49993999  2.25036733         Inf
eval_inverse(d2, -log(-log(at)))
#> [[1]]
#>  [1] 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
eval_deriv(d2, -log(-log(at)))
#> [[1]]
#>  [1]        NaN 0.23025851 0.32188758 0.36119184 0.36651629 0.34657359
#>  [7] 0.30649537 0.24967246 0.17851484 0.09482446        NaN

# can apply to distribution vectors
eval_transform(c(d2,d2,d2), 0.5)
#> [1] 0.3665129 0.3665129 0.3665129
eval_inverse(c(d2,d2,d2), 0.3665129)
#> [1] 0.5 0.5 0.5
eval_deriv(c(d2,d2,d2), 0.3665129)
#> [1] 0.3465736 0.3465736 0.3465736

# can be used to define tests
expect_correct_inverse <- function(dist, value) {
  tvalue <- eval_transform(dist, value)[[1]]
  inv_tvalue <- eval_inverse(dist, tvalue)[[1]]
  expect_equal(inv_tvalue, value)
}

expect_correct_inverse(d2, 0.5)

# also tests for derivatives
expect_correct_derivative <- function(dist, true_fun, value) {
  res <- eval_deriv(dist, value)[[1]]
  ref <- true_fun(value)
  expect_equal(res, ref, tolerance = 1e-5)
}

expect_correct_derivative(d2, function(x) exp(-exp(-x)) * exp(-x), 0.5)

^{Created on 2024-04-05 with reprex v2.1.0}

For example, now density.dist_transformed could be:

density.dist_transformed <- function(x, at, ...){
  inv <- function(at) suppressWarnings(eval_inverse(x, at))
  jacobian <- eval_deriv(x, at)
  d <- density(x[["dist"]], inv(at)) * abs(jacobian)
  limits <- field(support(x), "lim")[[1]]
  closed <- field(support(x), "closed")[[1]]
  if (!any(is.na(limits))) {
    `%less_than%` <- if (closed[1]) `<` else `<=`
    `%greater_than%` <- if (closed[2]) `>` else `>=`
    d[which(at %less_than% limits[1] | at %greater_than% limits[2])] <- 0
  }
  d
}

which will work with the current master branch, but also with the symbolic derivatives PR - no need for changing density.dist_transformed() if the procedure for derivatives changes - it always uses eval_deriv(x, at), which itself could be whatever function is desired (currently the call to numDeriv, soon the symbolic derivatives).

I just think this makes it a bit more future-proof

I can push a draft PR so you can see it in action