Merge pull request #46 from bayesplay/approximation

updated plot tests

DESCRIPTION

@@ -25,7 +25,7 @@ Imports:
     gginnards,
     stats
 Authors@R: 
-    person(given = c("Lincoln", "J"), 
+    person(given = c("Lincoln", "John"), 
             family = "Colling",
             comment = c(ORCID = "0000-0002-3572-7758"),
             email = "lincoln@colling.net.nz",

NAMESPACE

@@ -16,6 +16,7 @@ export(sd_ratio)
 export(visual_compare)
 exportMethods(summary)
 importFrom(methods,new)
+importFrom(methods,slot)
 importFrom(stats,integrate)
 importFrom(stats,pbeta)
 importFrom(stats,pcauchy)

R/bayesplay-package.R

@@ -3,6 +3,7 @@
 
 ## usethis namespace: start
 #' @importFrom methods new
+#' @importFrom methods slot
 #' @importFrom stats integrate
 #' @importFrom stats pbeta
 #' @importFrom stats pcauchy

R/helper_functions.R

@@ -25,10 +25,6 @@ in_range <- function(x, range) {
 
 
 
-`/.product` <- function(e1, e2) {
-  bf <- e1@data[["integral"]] / e2@data[["integral"]]
-  return(bf)
-}
 
 setOldClass("numeric")
 
@@ -64,7 +60,7 @@ integral <- function(obj) {
     stop("obj must be of class product", call. = FALSE)
   }
   rv <- new("auc", obj[["integral"]])
-  attr(rv, "approximate") <- obj[["approximate"]]
+  attr(rv, "approximate") <- slot(obj, "approximation")
   ## prior description
   p <- obj@prior_obj
   p_desc <- p[["parameters"]]
@@ -75,38 +71,29 @@ integral <- function(obj) {
 
 #' @export
 `/.auc` <- function(e1, e2) {
-  # FIXME: This is a hack... attribute should be set earlier
-  is_e1_approx <- attributes(e1)[["approximate"]] %||% FALSE
-  is_e2_approx <- attributes(e2)[["approximate"]] %||% FALSE
+
+  is_e1_approx <- is_approx(e1)
+  is_e2_approx <- is_approx(e2)
 
   has_approximation <- is_e1_approx | is_e2_approx
 
-  is_e1_point <- attributes(e1)[["prior"]][["family"]] == "point"
-  is_e2_point <- attributes(e2)[["prior"]][["family"]] == "point"
+  is_e1_point <- is_point(e1, 0L)
+  is_e2_point <- is_point(e2, 0L)
 
-  is_e1_point <- ifelse(is_e1_point,
-    attributes(e1)[["prior"]][["point"]] == 0L, FALSE
-  )
-  is_e2_point <- ifelse(is_e2_point,
-    attributes(e2)[["prior"]][["point"]] == 0L, FALSE
-  )
 
   one_is_point <- is_e1_point | is_e2_point
 
 
   # if one of the objects is an approximation
   # then one of the them must also be a point null
-  # FIXME: This condition is not correct
+
   if (has_approximation == TRUE && one_is_point == FALSE) {
     stop("Marginal likelihood is a approximation. One prior must ",
       "be a point prior at 0",
       call. = FALSE
     )
   }
 
-
-
-
   new("bf", unclass(e1) / unclass(e2))
 }
 
@@ -231,20 +218,20 @@ integer_breaks <- function(n = 5L, ...) {
 sd_ratio <- function(x, point) {
   is_estimated <- x@approximation %||% FALSE
   if (is_estimated && point != 0L) {
-    stop("point must be 0 if the marginal likelihood is estimated",
+    stop("point must be 0 if the marginal likelihood is an approximation",
       call. = FALSE
     )
   } else if (is_estimated && point == 0L) {
     bf <- x@data[["integral"]] / x@likelihood_obj@func(0L)
     return(new("bf", bf))
   }
-  # TODO: Move approximation it's own slot in class
+
   bf <- x@prior_obj@func(point) / x[["posterior_function"]](point)
   new("bf", bf)
 }
 
 
-### FIXME: This needs to be updated
+
 
 #' @export
 `*.bayesplay` <- function(e1, e2) {
@@ -258,7 +245,7 @@ sd_ratio <- function(x, point) {
 
 
   likelihood_family <- likelihood[["family"]]
-  # prior_family <- prior[["family"]]
+
 
   approx <- check_approximation(likelihood, prior)
   do_approximation <- approx[["approximation"]]
@@ -271,16 +258,17 @@ sd_ratio <- function(x, point) {
 
 
   if (do_approximation == TRUE && likelihood_family == "noncentral_t") {
-    stop("t value is large; approximation needed.
-    Reparametrise using a `noncentral_d` or `noncentral_d2` likelihood.",
+    stop(
+    "t value is large; approximation needed
+     Reparameterize using a `noncentral_d` of `noncentral_d2` likelihood.",
       call. = FALSE
     )
   }
 
 
   if (do_approximation == TRUE && supported_prior == FALSE) {
-    stop("Obseration is large; approximation needed.
-    Approximations are only supported with cauchy priors.", call. = FALSE)
+    warning("Observation is large; approximation needed.")
+    stop("Approximations are only supported with cauchy priors.", call. = FALSE)
   }
 
   if (do_approximation == TRUE) {
@@ -289,14 +277,10 @@ sd_ratio <- function(x, point) {
     df <- approximation_params[["df"]]
     approximation <- estimate_marginal(n, t, df, prior)
     marginal_likelihood_approx <- approximation[["marginal"]]
-    observation_shift <- approximation[["observation_shift"]]
-    if (observation_shift != 0L) {
-      obs <- recompute_observation(likelihood, observation_shift)
-      warning("Observation of d = ", obs[["d"]], " is unstable.\n",
-        "Shifting observation to d = ", obs[["new_d"]],
-        call. = FALSE
-      )
-    }
+
+
+
+
 
     warning("Observation is large; approximation needed.", call. = FALSE)
   }
@@ -403,7 +387,7 @@ reparameterise_d_to_t <- function(likelihood_obj) {
   d <- likelihood_obj[["parameters"]][["d"]]
   t <- d * sqrt(n)
   df <- n - 1L
-  # observation_type <- "d"
+
   list(n = n, t = t, df = df)
 }
 
@@ -486,50 +470,10 @@ check_approximation <- function(likelihood_obj, prior_obj) {
 }
 
 
-shift_observation <- function(likelihood_obj, shift) {
-  family <- likelihood_obj[["family"]]
-  params <- likelihood_obj[["parameters"]]
-  observation <- params[[1L]]
-  # round_observation <- round(observation, 2L)
-  round_observation <- observation
-  # fallback for second pass
-  # if (observation == round_observation) {
-  round_observation <- round_observation + shift
-  # }
-  params[[1L]] <- round_observation
-  params[["family"]] <- family
-  do.call(likelihood, params)
-}
-
-
-total_n <- function(likelihood_obj) {
-  if (likelihood_obj[["family"]] == "noncentral_d") {
-    n <- likelihood_obj[["parameters"]][["n"]]
-  } else if (likelihood_obj[["family"]] == "noncentral_d2") {
-    n1 <- likelihood_obj[["parameters"]][["n1"]]
-    n2 <- likelihood_obj[["parameters"]][["n2"]]
-    n <- n1 * n2 / (n1 + n2)
-  } else {
-    stop("This should not happen")
-  }
-  n
-}
-
-
-recompute_observation <- function(likelihood_obj, observation_shift) {
-  n <- total_n(likelihood_obj)
-  d <- likelihood_obj[["parameters"]][["d"]]
-  t <- d * sqrt(n)
-  new_t <- t + observation_shift
-  new_d <- new_t / sqrt(n)
-  list(d = d, new_d = new_d)
-}
-
-
 
 estimate_marginal <- function(n, t, df, prior) {
   prior_limits <- prior[["parameters"]][["range"]]
-  # prior_family <- prior[["family"]]
+
   upper <- prior_limits[[2L]]
   lower <- prior_limits[[1L]]
   var_delta <- 1L / n
@@ -567,66 +511,45 @@ estimate_marginal <- function(n, t, df, prior) {
   # The numerical integration can silently fail
   # This is extremely rare, but when it does occur we just need to
   # shift the observation by a small amount and try again
-  # TODO: Move up and down so that we can get a better estimate
-
-  # FIXME: Step out in both directions at the same time
-  # Then the one that stops first is the closest one
-  pass_1_likelihood <- new_likelihood
-  pass_2_likelihood <- new_likelihood
-  original_observation <- new_likelihood@observation
-  while (TRUE) {
-    auc_h1_pass1 <- integrate(
-      Vectorize(\(x)  pass_1_likelihood@func(x) * new_prior@func(x)),
-      -Inf, Inf,
-      subdivisions = 1000L, rel.tol = 1e-10, abs.tol = 1e-10
-    )
-    error1 <- auc_h1_pass1[["abs.error"]]
 
-    auc_h1_pass2 <- integrate(
-      Vectorize(\(x)  pass_2_likelihood@func(x) * new_prior@func(x)),
+
+
+
+
+
+
+
+    auc_h1_pass1 <- integrate(
+      Vectorize(\(x)  new_likelihood@func(x) * new_prior@func(x)),
       -Inf, Inf,
-      subdivisions = 1000L, rel.tol = 1e-10, abs.tol = 1e-10
+      subdivisions = 1000L, abs.tol = 1e-14
     )
-    error2 <- auc_h1_pass2[["abs.error"]]
 
 
-    error <- error1 == 0L | error2 == 0L
-    if (error) {
-      pass_1_likelihood <- shift_observation(pass_1_likelihood, -0.01)
-      pass_2_likelihood <- shift_observation(pass_2_likelihood, 0.01)
-    } else {
-      break
-    }
-  }
-
-  # TODO: This should pick the most conservative one
-  observation_shift_1 <- abs(
-    pass_1_likelihood@observation - original_observation
-  )
-  observation_shift_2 <- abs(
-    pass_2_likelihood@observation - original_observation
-  )
-  if (observation_shift_1 < observation_shift_2) {
     auc_h1 <- auc_h1_pass1
-    new_likelihood <- pass_1_likelihood
-  } else {
-    auc_h1 <- auc_h1_pass2
-    new_likelihood <- pass_2_likelihood
-  }
+    new_likelihood <- new_likelihood
 
 
   auc_h1 <- auc_h1[["value"]]
   auc_h0 <- new_likelihood@func(0L)
   log_bf_uncorrected <- log(auc_h1 / auc_h0)
   log_bf <- log_bf_interval + log_bf_uncorrected
   bf <- exp(log_bf)
-  new_observation <- new_likelihood@observation
 
-  if (new_observation != original_observation) {
-    observation_shift <- new_observation - original_observation
-  } else {
-    observation_shift <- 0L
-  }
 
-  list(marginal = bf * auc_h0, bf = bf, observation_shift = observation_shift)
+
+  list(marginal = bf * auc_h0, bf = bf)
+}
+
+
+is_approx <- function(e1) {
+  atr <- attributes(e1)
+  atr[["approximate"]]
+}
+
+is_point <- function(e1, value) {
+  if (attributes(e1)[["prior"]][["family"]] == "point") {
+   return(attributes(e1)[["prior"]][["point"]] == value)
+  }
+  FALSE
 }

R/priors.r

@@ -286,10 +286,10 @@ Adjust the mean or the limits.")
 }
 
 truncate_normalise_cauchy <- function(family, range, location, scale, ...) {
-  # nolint
-  unnormalised <- function(x) get_function(family)(x = x, location, scale) # nolint
-  # ll <- min(range)
-  # ul <- max(range)
+
+  unnormalised <- function(x) get_function(family)(x = x, location, scale)
+
+
 
   truncated_function <- function(x) {
     ifelse(in_range(x, range),
@@ -312,10 +312,6 @@ truncate_normalise_cauchy <- function(family, range, location, scale, ...) {
 }
 
 
-# truncate_normalise_beta <- function(family = family, range = range,
-# beta = beta, alpha = alpha) {
-#   stop("Error")
-# }
 
 
 truncate_normalise_beta <- function(family, range, alpha, beta, ...) {
@@ -332,12 +328,7 @@ truncate_normalise_beta <- function(family, range, alpha, beta, ...) {
   }
 
   k <- range_area_beta(alpha, beta, ll, ul)
-  if (k != 0L) {
-    constant <- 1L / k
-  } else {
-    warning("Could not normalise the truncated prior. Adjust the limits.")
-    constant <- 0L
-  }
+  constant <- 1L / k
 
 
   normalised <- function(x) truncated_function(x) * constant
@@ -403,10 +394,10 @@ make_prior.normal <- function(family, mean, sd, range = NULL) { # nolint
 #' @noRd
 make_prior.point <- function(family, point = 0L) { # nolint
   func <- function(x) get_function(family)(x = x, point = point)
-  # width <- 4L
-  # range <- c(point - width, point + width)
+
+
   params <- list(point = point)
-  # func <- make_distribution("point", list(point = point)) # nolint
+
   desc <- describe_prior(family, params)
 
   data <- make_prior_data(family, params, func)

cran-comments.md

@@ -1,3 +1,22 @@
+## Resubmission after archive
+
+Fixed warning on undocumented s3 methods
+
+
+## Test environments
+
+* R-hub builder macOS R-release
+* R-hub builder Debian Linux, R-release, GCC
+* R-hub builder Fedora Linux, R-devel, clang, gfortran
+* R win-builder R-devel 
+
+## R CMD check results
+
+1 NOTE
+New submission
+
+Package was archived on CRAN because of warning on undocumented S3 method
+
 ## Update
 
 Minor bug-fix to `make_likelihood.noncentral_d2()`