Merge pull request #316 from stan-dev/discrete_PIT_in_ppc_pit_ecdf

Discrete pit in ppc pit ecdf

Author: jgabry

DESCRIPTION

@@ -55,7 +55,7 @@ Suggests:
     survival,
     testthat (>= 2.0.0),
     vdiffr (>= 1.0.2)
-RoxygenNote: 7.2.3
+RoxygenNote: 7.3.0
 VignetteBuilder: knitr
 Encoding: UTF-8
 Roxygen: list(markdown = TRUE)

NEWS.md

@@ -2,6 +2,7 @@
 
 # bayesplot 1.10.0.9000
 
+* `ppc_pit_ecdf()` and `ppc_pit_ecdf_grouped()` now support discrete variables, and their default method for selecting the number of ECDF evaluation points has been updated.
 * Items for next release here
 
 # bayesplot 1.10.0

R/helpers-ppc.R

@@ -416,12 +416,12 @@ adjust_gamma_simulate <- function(N, L, K, prob, M) {
 #' simultaneous coverage of the ECDF traces.
 #' @noRd
 interpolate_gamma <- function(N, K, prob, L) {
-  # Find the precomputed values ueful for the interpolation task.
+  # Find the precomputed values useful for the interpolation task.
   vals <- get_interpolation_values(N, K, L, prob)
   # Largest lower bound and smalles upper bound for N among precomputed values.
   N_lb <- max(vals[vals$N <= N, ]$N)
   N_ub <- min(vals[vals$N >= N, ]$N)
-  # Approximate largest lower bound and smalles upper bound for gamma.
+  # Approximate largest lower bound and smallest upper bound for gamma.
   log_gamma_lb <- approx(
     x = log(vals[vals$N == N_lb, ]$K),
     y = log(vals[vals$N == N_lb, ]$val),

R/ppc-distributions.R

@@ -606,13 +606,13 @@ ppc_pit_ecdf <- function(y,
   if (is.null(pit)) {
     pit <- ppc_data(y, yrep) %>%
       group_by(.data$y_id) %>%
-      dplyr::group_map(~ mean(.x$value[.x$is_y] >= .x$value[!.x$is_y])) %>%
+      dplyr::group_map(
+        ~ mean(.x$value[.x$is_y] > .x$value[!.x$is_y]) +
+        runif(1, max = mean(.x$value[.x$is_y] == .x$value[!.x$is_y]))
+        ) %>%
       unlist()
     if (is.null(K)) {
-      K <- min(
-        length(unique(ppc_data(y, yrep)$rep_id)) + 1,
-        length(pit)
-      )
+      K <- nrow(yrep) + 1
     }
   } else {
     inform("'pit' specified so ignoring 'y', and 'yrep' if specified.")
@@ -632,20 +632,22 @@ ppc_pit_ecdf <- function(y,
   ggplot() +
     aes(
       x = 1:K / K,
-      y = ecdf(pit)(seq(0, 1, length.out = K)) - (plot_diff == TRUE) * 1:K / K,
+      y = ecdf(pit)(seq(0, 1, length.out = K)) -
+          (plot_diff == TRUE) * seq(0, 1, length.out = K),
       color = "y"
     ) +
     geom_step(show.legend = FALSE) +
     geom_step(aes(
-      y = lims$upper[-1] / N - (plot_diff == TRUE) * 1:K / K,
+      y = lims$upper[-1] / N - (plot_diff == TRUE) * seq(0, 1, length.out = K),
       color = "yrep"
-    ), show.legend = FALSE) +
+    ),
+    linetype = 2, show.legend = FALSE) +
     geom_step(aes(
-      y = lims$lower[-1] / N - (plot_diff == TRUE) * 1:K / K,
+      y = lims$lower[-1] / N - (plot_diff == TRUE) * seq(0, 1, length.out = K),
       color = "yrep"
-    ), show.legend = FALSE) +
-    yaxis_title(FALSE) +
-    xaxis_title(FALSE) +
+    ),
+    linetype = 2, show.legend = FALSE) +
+    labs(y = ifelse(plot_diff,"ECDF - difference","ECDF"), x = "PIT") +
     yaxis_ticks(FALSE) +
     scale_color_ppc() +
     bayesplot_theme_get()
@@ -671,10 +673,13 @@ ppc_pit_ecdf_grouped <-
     if (is.null(pit)) {
       pit <- ppc_data(y, yrep, group) %>%
         group_by(.data$y_id) %>%
-        dplyr::group_map(~ mean(.x$value[.x$is_y] >= .x$value[!.x$is_y])) %>%
+        dplyr::group_map(
+          ~ mean(.x$value[.x$is_y] > .x$value[!.x$is_y]) +
+          runif(1, max = mean(.x$value[.x$is_y] == .x$value[!.x$is_y]))
+          ) %>%
         unlist()
       if (is.null(K)) {
-        K <- length(unique(ppc_data(y, yrep)$rep_id)) + 1
+        K <- nrow(yrep) + 1
       }
     } else {
       inform("'pit' specified so ignoring 'y' and 'yrep' if specified.")
@@ -723,13 +728,14 @@ ppc_pit_ecdf_grouped <-
       geom_step(aes(
         y = .data$lims_upper - (plot_diff == TRUE) * .data$x,
         color = "yrep"
-      ), show.legend = FALSE) +
+      ),
+      linetype = 2, show.legend = FALSE) +
       geom_step(aes(
         y = .data$lims_lower - (plot_diff == TRUE) * .data$x,
         color = "yrep"
-      ), show.legend = FALSE) +
-      xaxis_title(FALSE) +
-      yaxis_title(FALSE) +
+      ),
+      linetype = 2, show.legend = FALSE) +
+      labs(y = ifelse(plot_diff,"ECDF - difference","ECDF"), x = "PIT") +
       yaxis_ticks(FALSE) +
       bayesplot_theme_get() +
       facet_wrap("group") +

R/sysdata.rda

@@ -1,7 +1,9 @@
 #' @param K An optional integer defining the number of equally spaced evaluation
-#'   points for the ECDF. If the submitted PIT values are known to be discrete,
-#'   this should be the number of the discrete cases. Defaults to the smaller of
-#'   `length(y)` and `ncol(yrep)` when applicable.
+#'   points for the ECDF. Reducing K when using `interpolate_adj = FALSE` makes
+#' computing the confidence bands faster. For `ppc_pit_ecdf` and
+#' `ppc_pit_ecdf_grouped`, defaults to `ncol(yrep) + 1`, or `length(pit)` if PIT
+#'  values are supplied. For `mcmc_rank_ecdf` defaults to the number of
+#'  iterations per chain in `x`.
 #' @param prob The desired simultaneous coverage level of the bands around the
 #'   ECDF. A value in (0,1).
 #' @param plot_diff A boolean defining whether to plot the difference between