add bins argument

[ci skip]

Author: jgabry

R/ppc-errors.R

@@ -311,18 +311,19 @@ ppc_error_scatter_avg_vs_x <-
 
 #' @rdname PPC-errors
 #' @export
-ppc_error_binned <- function(y, yrep, ..., size = 1, alpha = 0.25) {
+#' @param bins For \code{ppc_error_binned}, the number of bins to use (approximately).
+ppc_error_binned <- function(y, yrep, ..., bins = NULL, size = 1, alpha = 0.25) {
   check_ignored_arguments(...)
 
   y <- validate_y(y)
   yrep <- validate_yrep(yrep, y)
-  binned <- binned_error_data(y, yrep)
+  binned <- binned_error_data(y, yrep, bins = bins)
 
   mixed_scheme <- is_mixed_scheme(color_scheme_get())
   point_fill <- get_color(ifelse(mixed_scheme, "m", "d"))
   point_color <- get_color(ifelse(mixed_scheme, "mh", "dh"))
   graph <-
-    ggplot(binned, aes_(x = ~ xbar)) +
+    ggplot(binned, aes_(x = ~ ey_bar)) +
     geom_hline(
       yintercept = 0,
       linetype = 2,
@@ -345,7 +346,7 @@ ppc_error_binned <- function(y, yrep, ..., size = 1, alpha = 0.25) {
       size = size
     ) +
     geom_point(
-      mapping = aes_(y = ~ ybar),
+      mapping = aes_(y = ~ err_bar),
       shape = 21,
       fill = point_fill,
       color = point_color
@@ -388,21 +389,48 @@ grouped_error_data <- function(y, yrep, group) {
   }
   dat <- dplyr::bind_rows(errs)
   dat$y_id <- NULL
-  dat
+  return(dat)
 }
 
-bin_errors <- function(rep_id, ey, r, nbins) {
+binned_error_data <- function(y, yrep, bins = NULL) {
+  if (is.null(bins)) {
+    bins <- n_bins(length(y))
+  }
+
+  errors <- compute_errors(y, yrep)
+  binned_errs <- list()
+  for (s in 1:nrow(errors)) {
+    binned_errs[[s]] <- bin_errors(ey = yrep[s,], r = errors[s,], bins = bins,
+                                   rep_id = s)
+  }
+  dat <- dplyr::bind_rows(binned_errs)
+  return(dat)
+}
+
+# calculate number of bins binned_error_data()
+# @parmam N Number of data points, length(y)
+n_bins <- function(N) {
+  if (N <= 10) {
+    return(floor(N / 2))
+  } else if (N > 10 && N < 100) {
+    return(10)
+  } else { # N >= 100
+    return(floor(sqrt(N)))
+  }
+}
+
+bin_errors <- function(ey, r, bins, rep_id = NULL) {
   N <- length(ey)
-  break_ids <- floor(N * (1:(nbins - 1)) / nbins)
+  break_ids <- floor(N * (1:(bins - 1)) / bins)
   if (any(break_ids == 0)) {
-    nbins <- 1
+    bins <- 1
   }
-  if (nbins == 1) {
+  if (bins == 1) {
     breaks <- c(-Inf, sum(range(ey)) / 2, Inf)
   } else {
     ey_sort <- sort(ey)
     breaks <- -Inf
-    for (i in 1:(nbins - 1)) {
+    for (i in 1:(bins - 1)) {
       break_i <- break_ids[i]
       ey_range <- ey_sort[c(break_i, break_i + 1)]
       if (diff(ey_range) == 0) {
@@ -417,43 +445,21 @@ bin_errors <- function(rep_id, ey, r, nbins) {
     breaks <- unique(c(breaks, Inf))
   }
 
-  nbins <- length(breaks) - 1
   ey_binned <- as.numeric(cut(ey, breaks))
-
-  out <- matrix(NA, nrow = nbins, ncol = 3)
-  for (i in 1:nbins) {
+  bins <- length(breaks) - 1
+  out <- matrix(NA, nrow = bins, ncol = 4)
+  colnames(out) <- c("ey_bar", "err_bar", "se2", "bin")
+  for (i in 1:bins) {
     mark <- which(ey_binned == i)
     ey_bar <- mean(ey[mark])
     r_bar <- mean(r[mark])
     s <- if (length(r[mark]) > 1) sd(r[mark]) else 0
-    out[i, ] <- c(ey_bar, r_bar, 2 * s / sqrt(length(mark)))
+    out[i, ] <- c(ey_bar, r_bar, 2 * s / sqrt(length(mark)), i)
   }
   out <- as.data.frame(out)
-  colnames(out) <- c("xbar", "ybar", "se2")
-  out$rep_id <- as.integer(rep_id)
-  return(out)
-}
-
-
-binned_error_data <- function(y, yrep) {
-  N <- length(y)
-  if (N >= 100) {
-    nbins <- floor(sqrt(N))
-  } else if (N > 10 && N < 100) {
-    nbins <- 10
-  } else {
-    # if (N <= 10)
-    nbins <- floor(N / 2)
+  if (!is.null(rep_id)) {
+    out$rep_id <- as.integer(rep_id)
   }
-
-  errors <- compute_errors(y, yrep)
-  binned_df <- bin_errors(rep_id = 1, ey = yrep[1, ], r = errors[1, ], nbins = nbins)
-  if (nrow(errors) > 1) {
-    for (s in 2:nrow(errors)) {
-      binned_s <- bin_errors(rep_id = s, ey = yrep[s,], r = errors[s,], nbins = nbins)
-      binned_df <- rbind(binned_df, binned_s)
-    }
-  }
-  return(binned_df)
+  return(out)
 }