From c776c9ac196c1dfd2c826d0d50437ac404bdd9bb Mon Sep 17 00:00:00 2001
From: Joshua Nathaniel Pritikin <jpritikin@pobox.com>
Date: Fri, 17 May 2019 07:53:23 -0400
Subject: [PATCH 1/7] Add sbc_hist

https://github.com/stan-dev/bayesplot/issues/146
---
 NAMESPACE |   1 +
 R/sbc.R   | 106 ++++++++++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 107 insertions(+)
 create mode 100644 R/sbc.R

diff --git a/NAMESPACE b/NAMESPACE
index 0436c4cb..7d557ce7 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -126,6 +126,7 @@ export(ppc_stat_freqpoly_grouped)
 export(ppc_stat_grouped)
 export(ppc_violin_grouped)
 export(rhat)
+export(sbc_hist)
 export(scatter_style_np)
 export(theme_default)
 export(trace_style_np)
diff --git a/R/sbc.R b/R/sbc.R
new file mode 100644
index 00000000..d5412268
--- /dev/null
+++ b/R/sbc.R
@@ -0,0 +1,106 @@
+klUniform <- function(v, numPriorDraws, samplesPerPrior) {
+  # D_{KL}(v || uniform)
+  # https://en.wikipedia.org/wiki/Kullback%E2%80%93Leibler_divergence
+  expectedPr <- 1.0/samplesPerPrior
+  observedPr <- table(v) / numPriorDraws
+  sum(observedPr * log(observedPr/expectedPr))
+}
+
+#' Histograms for Simulation Based Calibration
+#'
+#' @param ranks A list of sampling realizations.
+#' @param thin An integer vector of length one indicating the thinning interval
+#' when plotting
+#' @param perBin Number of histogram entries to combine into a single bar.
+#' @param worst If NA, plots all parameters. Otherwise how many parameters to show.
+#' Parameters are ordered by the degree of non-uniformity.
+#' @param alpha Uncertainty interval probability for a false positive (alpha level).
+#' @param hideAxes Whether to hide the plot axes.
+#'
+#' Each list element of \code{ranks} should be a matrix of rank
+#' comparison results (encoded as 0 or 1) associated with a single
+#' draw from the prior distribution. Each draw from the posterior is
+#' in the row and parameters are in columns. The matrix should have
+#' column names to correctly label the parameters.
+#'
+#' So that the histograms consist of independent realizations,
+#' draws from the posterior should be thinned to remove
+#' autocorrelation. Set \code{thin} such that the number of
+#' draws approximately matches the effective sample size.
+#'
+#' For best results, one plus the number of draws from the posterior
+#' should be evenly divisible by the number of histogram bins after
+#' thinning. For example, 511 draws after thinning results in 128
+#' draws. If perBin is set to 4 then 32 histogram bars are drawn.
+#' 
+#' @template return-ggplot
+#'
+#' @references
+#'   Talts, S., Betancourt, M., Simpson, D., Vehtari, A., and Gelman, A. (2018).
+#'   Validating Bayesian Inference Algorithms with Simulation-Based Calibration.
+#'   arXiv preprint arXiv:1804.06788. \url{https://arxiv.org/abs/1804.06788}
+#' @seealso
+#' \link[rstan]{sbc}
+#' @examples
+#' pars <- paste0('parameter',1:2)
+#' samplesPerPrior <- 511
+#' ranks <- list()
+#' for (px in 1:500) {
+#'   r1 <- matrix(0, nrow=samplesPerPrior, ncol=length(pars),
+#'          dimnames=list(NULL, pars))
+#'   for (p1 in 1:length(pars)) {
+#'     r1[sample.int(samplesPerPrior,
+#'                   floor(runif(1, 0, samplesPerPrior))), p1] <- 1
+#'   }
+#'   ranks[[px]] <- r1
+#' }
+#' sbc_hist(ranks)
+#' @export
+
+sbc_hist <- function(ranks, thin = 4, perBin=4, worst=16, ...,
+                     alpha = 0.01, hideAxes=TRUE) {
+  numPriorDraws <- length(ranks)
+  thinner <- seq(from = 1, to = nrow(ranks[[1]]), by = thin)
+  samplesPerPrior <- length(thinner)
+  u <- t(sapply(ranks, FUN = function(r) 1 + colSums(r[thinner, , drop = FALSE])))
+  if (ncol(ranks[[1]]) == 1) {
+    u <- t(u)
+    dimnames(u) <- list(NULL, colnames(ranks[[1]]))
+  }
+  
+  if (!is.na(worst)) {
+    kl <- apply(u, 2, function(v) klUniform(v, numPriorDraws, samplesPerPrior))
+    filter <- order(-kl)[1:min(worst,ncol(u))]
+#    print(filter)
+    u <- u[, filter, drop=FALSE ]
+  }
+  
+  parameter <- ordered(rep(colnames(u), each = nrow(u)),
+                       levels=colnames(u))
+  d <- data.frame(u = c(u), parameter)
+  if (samplesPerPrior %% perBin != 0) {
+    warning(paste("perBin (", perBin, ") does not evenly divide the",
+      "number of samples per prior (",samplesPerPrior,")"))
+  }
+  numBins <- samplesPerPrior/perBin
+  CI <- qbinom(c(alpha/2,0.5,1-alpha/2), numPriorDraws, numBins^-1) + c(-.5,0,.5)
+  offset <- perBin*2
+  pl <- ggplot(d, aes(x = u)) + 
+    geom_polygon(data=data.frame(x=c(-offset,0,-offset,samplesPerPrior + offset,
+      samplesPerPrior, samplesPerPrior + offset,-offset),
+      y=c(CI[1],CI[2],CI[3],CI[3],CI[2],CI[1],CI[1])),
+      aes(x=x,y=y),fill="grey45",color="grey25",alpha=0.5) +
+    geom_histogram(bins=numBins, na.rm=TRUE) +
+#    xlim(1,samplesPerPrior) +
+    # https://github.com/tidyverse/ggplot2/issues/3332
+    facet_wrap("parameter") +
+    geom_hline(yintercept=CI[1], color='green', linetype="dotted", alpha=.5) +
+  geom_hline(yintercept=CI[3], color='green', linetype="dotted", alpha=.5)
+  if (hideAxes) {
+    pl <- pl + theme(axis.text.x=element_blank(),
+      axis.text.y=element_blank(),axis.ticks=element_blank(),
+      axis.title.x=element_blank(),
+      axis.title.y=element_blank())
+    }
+  pl
+}

From 79c4d5acef3c5c762369c203135d5300b2212f13 Mon Sep 17 00:00:00 2001
From: jgabry <jgabry@gmail.com>
Date: Fri, 17 May 2019 18:21:11 -0400
Subject: [PATCH 2/7] edits to sbc_hist to fit in with bayesplot conventions

---
 R/sbc.R         | 233 +++++++++++++++++++++++++++++++-----------------
 man/sbc_hist.Rd |  90 +++++++++++++++++++
 2 files changed, 241 insertions(+), 82 deletions(-)
 create mode 100644 man/sbc_hist.Rd

diff --git a/R/sbc.R b/R/sbc.R
index d5412268..e596560f 100644
--- a/R/sbc.R
+++ b/R/sbc.R
@@ -1,106 +1,175 @@
-klUniform <- function(v, numPriorDraws, samplesPerPrior) {
-  # D_{KL}(v || uniform)
-  # https://en.wikipedia.org/wiki/Kullback%E2%80%93Leibler_divergence
-  expectedPr <- 1.0/samplesPerPrior
-  observedPr <- table(v) / numPriorDraws
-  sum(observedPr * log(observedPr/expectedPr))
-}
-
-#' Histograms for Simulation Based Calibration
+#' Diagnostics for Simulation Based Calibration
 #'
-#' @param ranks A list of sampling realizations.
-#' @param thin An integer vector of length one indicating the thinning interval
-#' when plotting
-#' @param perBin Number of histogram entries to combine into a single bar.
-#' @param worst If NA, plots all parameters. Otherwise how many parameters to show.
-#' Parameters are ordered by the degree of non-uniformity.
-#' @param alpha Uncertainty interval probability for a false positive (alpha level).
-#' @param hideAxes Whether to hide the plot axes.
+#' @export
+#' @param ranks A list of matrices. See **Details**.
+#' @param thin An integer indicating the thinning interval to use when plotting
+#'   so that the histograms consist of (close to) independent realizations. Set
+#'   the `thin` argument such that the resulting number of draws approximately
+#'   matches the effective sample size.
+#' @param per_bin An integer indicating the number of histogram entries to
+#'   combine into a single bar. For best results, one plus the number of draws
+#'   from the posterior should be evenly divisible by the number of histogram
+#'   bins after thinning. For example, if there are `511` posterior draws (per
+#'   matrix in `ranks`) and `thin=4`, then after thinning there will be `128`
+#'   draws. If `per_bin=4`, then `128/4=32` histogram bars will be drawn.
+#' @param worst How many parameters to show, or `NA` to plot all parameters. If
+#'   `worst` is not `NA` then parameters are ordered by the degree of
+#'   non-uniformity so, for example, `worst = 10` means to plot only the `10`
+#'   worst parameters.
+#' @param prob The size of the interval plotted to show the expected behavior
+#'   under uniformity. The default is `prob=0.99`.
+#' @template args-facet_args
 #'
-#' Each list element of \code{ranks} should be a matrix of rank
-#' comparison results (encoded as 0 or 1) associated with a single
-#' draw from the prior distribution. Each draw from the posterior is
-#' in the row and parameters are in columns. The matrix should have
-#' column names to correctly label the parameters.
+#' @details
+#' Each element of the list `ranks` should be a matrix of rank comparison
+#' results (encoded as `0` or `1`) associated with a single draw from the prior
+#' distribution. (These are not actually "ranks" but can be used afterwards to
+#' reconstruct (thinned) ranks.) The columns of each matrix are model parameters
+#' and the rows are posterior draws from the model fit to the data generated
+#' from the corresponding prior draw. An element is `0` or `1` depending on
+#' whether the posterior draw is greater or less than the corresponding "true"
+#' realization (from the prior). The matrices should have column names to
+#' correctly label the parameters.
 #'
-#' So that the histograms consist of independent realizations,
-#' draws from the posterior should be thinned to remove
-#' autocorrelation. Set \code{thin} such that the number of
-#' draws approximately matches the effective sample size.
+#' @section Plot Descriptions:
+#' \describe{
+#'   \item{`sbc_hist()`}{
+#'   SBC histogram from Talts et al. (2018). A separate plot (facet) is created
+#'   for each parameter.
+#'   }
+#' }
 #'
-#' For best results, one plus the number of draws from the posterior
-#' should be evenly divisible by the number of histogram bins after
-#' thinning. For example, 511 draws after thinning results in 128
-#' draws. If perBin is set to 4 then 32 histogram bars are drawn.
-#' 
 #' @template return-ggplot
 #'
 #' @references
-#'   Talts, S., Betancourt, M., Simpson, D., Vehtari, A., and Gelman, A. (2018).
-#'   Validating Bayesian Inference Algorithms with Simulation-Based Calibration.
-#'   arXiv preprint arXiv:1804.06788. \url{https://arxiv.org/abs/1804.06788}
-#' @seealso
-#' \link[rstan]{sbc}
+#' Talts, S., Betancourt, M., Simpson, D., Vehtari, A., and Gelman, A. (2018).
+#' Validating Bayesian Inference Algorithms with Simulation-Based Calibration.
+#' arXiv preprint arXiv:1804.06788. \url{https://arxiv.org/abs/1804.06788}
+#'
+#' @seealso [rstan::sbc()]
+#'
 #' @examples
-#' pars <- paste0('parameter',1:2)
-#' samplesPerPrior <- 511
+#' # create some fake inputs to use for sbc_hist()
+#' set.seed(19)
+#' pars <- paste0("beta[", 1:4, "]")
+#' samples_per_prior <- 511
+#' n_replications <- 500
 #' ranks <- list()
-#' for (px in 1:500) {
-#'   r1 <- matrix(0, nrow=samplesPerPrior, ncol=length(pars),
-#'          dimnames=list(NULL, pars))
+#' for (n in 1:n_replications) {
+#'   r1 <- matrix(0, nrow=samples_per_prior, ncol=length(pars),
+#'                dimnames=list(NULL, pars))
 #'   for (p1 in 1:length(pars)) {
-#'     r1[sample.int(samplesPerPrior,
-#'                   floor(runif(1, 0, samplesPerPrior))), p1] <- 1
+#'     r1[sample.int(samples_per_prior, floor(runif(1, 0, samples_per_prior))), p1] <- 1
 #'   }
-#'   ranks[[px]] <- r1
+#'   ranks[[n]] <- r1
 #' }
+#'
+#' color_scheme_set("purple")
 #' sbc_hist(ranks)
-#' @export
+#'
+sbc_hist <- function(ranks, ...,
+                     thin = 4, per_bin = 4, worst = 16, prob = 0.99,
+                     facet_args = list()) {
+  check_ignored_arguments(...)
+  stopifnot(is.list(ranks), all(sapply(ranks, is.matrix)))
+  rows <- sapply(ranks, nrow)
+  cols <- sapply(ranks, ncol)
+  if (any(rows != rows[1]) || any(cols != cols[1])) {
+    stop("Not all matrices in 'ranks' have the same dimensions.")
+  }
 
-sbc_hist <- function(ranks, thin = 4, perBin=4, worst=16, ...,
-                     alpha = 0.01, hideAxes=TRUE) {
-  numPriorDraws <- length(ranks)
+  num_prior_draws <- length(ranks)
   thinner <- seq(from = 1, to = nrow(ranks[[1]]), by = thin)
-  samplesPerPrior <- length(thinner)
+  samples_per_prior <- length(thinner)
   u <- t(sapply(ranks, FUN = function(r) 1 + colSums(r[thinner, , drop = FALSE])))
   if (ncol(ranks[[1]]) == 1) {
     u <- t(u)
     dimnames(u) <- list(NULL, colnames(ranks[[1]]))
   }
-  
+
   if (!is.na(worst)) {
-    kl <- apply(u, 2, function(v) klUniform(v, numPriorDraws, samplesPerPrior))
-    filter <- order(-kl)[1:min(worst,ncol(u))]
-#    print(filter)
-    u <- u[, filter, drop=FALSE ]
+    # order starting with worst (least uniform)
+    kl <- apply(u, 2, function(v) kl_uniform(v, num_prior_draws, samples_per_prior))
+    filter <- order(-kl)[1:min(worst, ncol(u))]
+    u <- u[, filter, drop = FALSE]
   }
-  
-  parameter <- ordered(rep(colnames(u), each = nrow(u)),
-                       levels=colnames(u))
-  d <- data.frame(u = c(u), parameter)
-  if (samplesPerPrior %% perBin != 0) {
-    warning(paste("perBin (", perBin, ") does not evenly divide the",
-      "number of samples per prior (",samplesPerPrior,")"))
+
+  pars <- parameter_names(u)
+  data <- data.frame(
+    Parameter = ordered(rep(pars, each = nrow(u)), levels = pars),
+    u = c(u)
+  )
+
+  num_bins <- samples_per_prior / per_bin
+  if (samples_per_prior %% per_bin != 0) {
+    warning(paste0("per_bin (", per_bin, ") does not evenly divide the",
+                   "number of samples per prior (", samples_per_prior ,")."))
   }
-  numBins <- samplesPerPrior/perBin
-  CI <- qbinom(c(alpha/2,0.5,1-alpha/2), numPriorDraws, numBins^-1) + c(-.5,0,.5)
-  offset <- perBin*2
-  pl <- ggplot(d, aes(x = u)) + 
-    geom_polygon(data=data.frame(x=c(-offset,0,-offset,samplesPerPrior + offset,
-      samplesPerPrior, samplesPerPrior + offset,-offset),
-      y=c(CI[1],CI[2],CI[3],CI[3],CI[2],CI[1],CI[1])),
-      aes(x=x,y=y),fill="grey45",color="grey25",alpha=0.5) +
-    geom_histogram(bins=numBins, na.rm=TRUE) +
-#    xlim(1,samplesPerPrior) +
-    # https://github.com/tidyverse/ggplot2/issues/3332
-    facet_wrap("parameter") +
-    geom_hline(yintercept=CI[1], color='green', linetype="dotted", alpha=.5) +
-  geom_hline(yintercept=CI[3], color='green', linetype="dotted", alpha=.5)
-  if (hideAxes) {
-    pl <- pl + theme(axis.text.x=element_blank(),
-      axis.text.y=element_blank(),axis.ticks=element_blank(),
-      axis.title.x=element_blank(),
-      axis.title.y=element_blank())
-    }
-  pl
+
+  # data for polygon showing expected behavior under uniformity
+  alpha <- 1 - prob
+  CI <- qbinom(
+    p = c(alpha / 2, 0.5, 1 - alpha / 2),
+    size = num_prior_draws,
+    prob = 1 / num_bins
+  )
+  CI <- CI + c(-.5, 0, .5)
+  offset <- 2 * per_bin
+  polygon_data <- data.frame(
+    x= c(-offset, 0, -offset, samples_per_prior + offset,
+         samples_per_prior, samples_per_prior + offset, -offset),
+    y = c(CI[1], CI[2], CI[3], CI[3], CI[2], CI[1], CI[1])
+  )
+
+  graph <- ggplot(data, aes_(x = ~ u)) +
+    geom_polygon(
+      aes_(x = ~ x, y = ~ y),
+      data = polygon_data,
+      fill = "lightgray",
+      color = NA,
+      alpha = 1
+    ) +
+    geom_segment(
+      x = 0,
+      xend = samples_per_prior,
+      y = CI[2],
+      yend = CI[2],
+      color = "darkgray",
+      alpha = 0.5,
+      size = 0.1
+    ) +
+    geom_histogram(
+      bins = num_bins,
+      fill = get_color("mid"),
+      color = get_color("mid_highlight"),
+      size = .25,
+      na.rm = TRUE
+    ) +
+    xlab("Rank statistic") +
+    coord_cartesian(
+      # xlim = c(-per_bin, samples_per_prior + per_bin),
+      expand = FALSE
+    )
+
+  facet_args[["facets"]] <- ~ Parameter
+  graph +
+    do.call("facet_wrap", facet_args) +
+    bayesplot_theme_get() +
+    xaxis_text(FALSE) +
+    yaxis_text(FALSE) +
+    xaxis_title(TRUE) +
+    yaxis_title(FALSE) +
+    xaxis_ticks(FALSE) +
+    yaxis_ticks(FALSE)
+}
+
+
+
+# internal ----------------------------------------------------------------
+kl_uniform <- function(v, num_prior_draws, samples_per_prior) {
+  # D_{KL}(v || uniform)
+  # https://en.wikipedia.org/wiki/Kullback%E2%80%93Leibler_divergence
+  expected_pr <- 1.0/samples_per_prior
+  observed_pr <- table(v) / num_prior_draws
+  sum(observed_pr * log(observed_pr/expected_pr))
 }
diff --git a/man/sbc_hist.Rd b/man/sbc_hist.Rd
new file mode 100644
index 00000000..567c72a6
--- /dev/null
+++ b/man/sbc_hist.Rd
@@ -0,0 +1,90 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/sbc.R
+\name{sbc_hist}
+\alias{sbc_hist}
+\title{Diagnostics for Simulation Based Calibration}
+\usage{
+sbc_hist(ranks, ..., thin = 4, per_bin = 4, worst = 16,
+  prob = 0.99, facet_args = list())
+}
+\arguments{
+\item{ranks}{A list of matrices. See \strong{Details}.}
+
+\item{thin}{An integer indicating the thinning interval to use when plotting
+so that the histograms consist of (close to) independent realizations. Set
+the \code{thin} argument such that the resulting number of draws approximately
+matches the effective sample size.}
+
+\item{per_bin}{An integer indicating the number of histogram entries to
+combine into a single bar. For best results, one plus the number of draws
+from the posterior should be evenly divisible by the number of histogram
+bins after thinning. For example, if there are \code{511} posterior draws (per
+matrix in \code{ranks}) and \code{thin=4}, then after thinning there will be \code{128}
+draws. If \code{per_bin=4}, then \code{128/4=32} histogram bars will be drawn.}
+
+\item{worst}{How many parameters to show, or \code{NA} to plot all parameters. If
+\code{worst} is not \code{NA} then parameters are ordered by the degree of
+non-uniformity so, for example, \code{worst = 10} means to plot only the \code{10}
+worst parameters.}
+
+\item{prob}{The size of the interval plotted to show the expected behavior
+under uniformity. The default is \code{prob=0.99}.}
+
+\item{facet_args}{A named list of arguments (other than \code{facets}) passed
+to \code{\link[ggplot2:facet_wrap]{ggplot2::facet_wrap()}} or \code{\link[ggplot2:facet_grid]{ggplot2::facet_grid()}}
+to control faceting.}
+}
+\value{
+A ggplot object that can be further customized using the \strong{ggplot2} package.
+}
+\description{
+Diagnostics for Simulation Based Calibration
+}
+\details{
+Each element of the list \code{ranks} should be a matrix of rank comparison
+results (encoded as \code{0} or \code{1}) associated with a single draw from the prior
+distribution. The columns of each matrix are model parameters and the rows
+are posterior draws from the model fit to the data generated from the
+corresponding prior draw. An element is \code{0} or \code{1} depending on whether the
+posterior draw is greater or less than the corresponding "true" realization
+(from the prior). The matrices should have column names to correctly label
+the parameters.
+}
+\section{Plot Descriptions}{
+
+\describe{
+\item{\code{sbc_hist()}}{
+SBC histogram from Talts et al. (2018). A separate plot (facet) is created
+for each parameter.
+}
+}
+}
+
+\examples{
+# create some fake inputs to use for sbc_hist()
+set.seed(19)
+pars <- paste0("beta[", 1:4, "]")
+samples_per_prior <- 511
+n_replications <- 500
+ranks <- list()
+for (n in 1:n_replications) {
+  r1 <- matrix(0, nrow=samples_per_prior, ncol=length(pars),
+               dimnames=list(NULL, pars))
+  for (p1 in 1:length(pars)) {
+    r1[sample.int(samples_per_prior, floor(runif(1, 0, samples_per_prior))), p1] <- 1
+  }
+  ranks[[n]] <- r1
+}
+
+color_scheme_set("purple")
+sbc_hist(ranks)
+
+}
+\references{
+Talts, S., Betancourt, M., Simpson, D., Vehtari, A., and Gelman, A. (2018).
+Validating Bayesian Inference Algorithms with Simulation-Based Calibration.
+arXiv preprint arXiv:1804.06788. \url{https://arxiv.org/abs/1804.06788}
+}
+\seealso{
+\code{\link[rstan:sbc]{rstan::sbc()}}
+}

From 8c8d910e07c9f66833d3f2c46826b0ce4b4c58a9 Mon Sep 17 00:00:00 2001
From: jgabry <jgabry@gmail.com>
Date: Fri, 17 May 2019 18:35:09 -0400
Subject: [PATCH 3/7] fixes for r cmd check

---
 R/sbc.R         |  3 ++-
 man/sbc_hist.Rd | 17 ++++++++++-------
 2 files changed, 12 insertions(+), 8 deletions(-)

diff --git a/R/sbc.R b/R/sbc.R
index e596560f..a33c4372 100644
--- a/R/sbc.R
+++ b/R/sbc.R
@@ -18,6 +18,7 @@
 #'   worst parameters.
 #' @param prob The size of the interval plotted to show the expected behavior
 #'   under uniformity. The default is `prob=0.99`.
+#' @param ... Currently ignored.
 #' @template args-facet_args
 #'
 #' @details
@@ -46,7 +47,7 @@
 #' Validating Bayesian Inference Algorithms with Simulation-Based Calibration.
 #' arXiv preprint arXiv:1804.06788. \url{https://arxiv.org/abs/1804.06788}
 #'
-#' @seealso [rstan::sbc()]
+#' @seealso `rstan::sbc()`
 #'
 #' @examples
 #' # create some fake inputs to use for sbc_hist()
diff --git a/man/sbc_hist.Rd b/man/sbc_hist.Rd
index 567c72a6..9e6a2515 100644
--- a/man/sbc_hist.Rd
+++ b/man/sbc_hist.Rd
@@ -10,6 +10,8 @@ sbc_hist(ranks, ..., thin = 4, per_bin = 4, worst = 16,
 \arguments{
 \item{ranks}{A list of matrices. See \strong{Details}.}
 
+\item{...}{Currently ignored.}
+
 \item{thin}{An integer indicating the thinning interval to use when plotting
 so that the histograms consist of (close to) independent realizations. Set
 the \code{thin} argument such that the resulting number of draws approximately
@@ -43,12 +45,13 @@ Diagnostics for Simulation Based Calibration
 \details{
 Each element of the list \code{ranks} should be a matrix of rank comparison
 results (encoded as \code{0} or \code{1}) associated with a single draw from the prior
-distribution. The columns of each matrix are model parameters and the rows
-are posterior draws from the model fit to the data generated from the
-corresponding prior draw. An element is \code{0} or \code{1} depending on whether the
-posterior draw is greater or less than the corresponding "true" realization
-(from the prior). The matrices should have column names to correctly label
-the parameters.
+distribution. (These are not actually "ranks" but can be used afterwards to
+reconstruct (thinned) ranks.) The columns of each matrix are model parameters
+and the rows are posterior draws from the model fit to the data generated
+from the corresponding prior draw. An element is \code{0} or \code{1} depending on
+whether the posterior draw is greater or less than the corresponding "true"
+realization (from the prior). The matrices should have column names to
+correctly label the parameters.
 }
 \section{Plot Descriptions}{
 
@@ -86,5 +89,5 @@ Validating Bayesian Inference Algorithms with Simulation-Based Calibration.
 arXiv preprint arXiv:1804.06788. \url{https://arxiv.org/abs/1804.06788}
 }
 \seealso{
-\code{\link[rstan:sbc]{rstan::sbc()}}
+\code{rstan::sbc()}
 }

From 45fc425c3b3e003778e5dd1f895b4912ae2599e1 Mon Sep 17 00:00:00 2001
From: jgabry <jgabry@gmail.com>
Date: Fri, 17 May 2019 18:40:55 -0400
Subject: [PATCH 4/7] missing whitespace

---
 R/sbc.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/R/sbc.R b/R/sbc.R
index a33c4372..1e5a4a1d 100644
--- a/R/sbc.R
+++ b/R/sbc.R
@@ -103,7 +103,7 @@ sbc_hist <- function(ranks, ...,
 
   num_bins <- samples_per_prior / per_bin
   if (samples_per_prior %% per_bin != 0) {
-    warning(paste0("per_bin (", per_bin, ") does not evenly divide the",
+    warning(paste0("per_bin (", per_bin, ") does not evenly divide the ",
                    "number of samples per prior (", samples_per_prior ,")."))
   }
 

From 577fff0f6b4067e38535f9e40ca3cba3e97b79a4 Mon Sep 17 00:00:00 2001
From: jgabry <jgabry@gmail.com>
Date: Fri, 17 May 2019 20:07:37 -0400
Subject: [PATCH 5/7] fixes to sbc_hist()

* remove segment in center of CI
* add segments at extremes of CI (lower one is in front of histogram, upper is behind, I think it looks nice)
* add x-axis ticks and labels at 0, samples_per_prior/2,samples_per_prior (this does convey some information)
---
 R/sbc.R | 35 +++++++++++++++++++++--------------
 1 file changed, 21 insertions(+), 14 deletions(-)

diff --git a/R/sbc.R b/R/sbc.R
index 1e5a4a1d..a5aacc81 100644
--- a/R/sbc.R
+++ b/R/sbc.R
@@ -69,7 +69,8 @@
 #' sbc_hist(ranks)
 #'
 sbc_hist <- function(ranks, ...,
-                     thin = 4, per_bin = 4, worst = 16, prob = 0.99,
+                     thin = 4, per_bin = 4, worst = 16,
+                     prob = 0.99,
                      facet_args = list()) {
   check_ignored_arguments(...)
   stopifnot(is.list(ranks), all(sapply(ranks, is.matrix)))
@@ -131,12 +132,12 @@ sbc_hist <- function(ranks, ...,
       alpha = 1
     ) +
     geom_segment(
-      x = 0,
-      xend = samples_per_prior,
-      y = CI[2],
-      yend = CI[2],
+      x = -offset,
+      xend = samples_per_prior + offset,
+      y = CI[3],
+      yend = CI[3],
       color = "darkgray",
-      alpha = 0.5,
+      alpha = 0.25,
       size = 0.1
     ) +
     geom_histogram(
@@ -146,21 +147,27 @@ sbc_hist <- function(ranks, ...,
       size = .25,
       na.rm = TRUE
     ) +
-    xlab("Rank statistic") +
-    coord_cartesian(
-      # xlim = c(-per_bin, samples_per_prior + per_bin),
-      expand = FALSE
-    )
+    geom_segment(
+      x = -offset,
+      xend = samples_per_prior + offset,
+      y = CI[1],
+      yend = CI[1],
+      color = "darkgray",
+      alpha = 0.1,
+      size = 0.1
+    ) +
+    scale_x_continuous(
+      name = "Rank statistic",
+      breaks = c(0, round(samples_per_prior / 2), samples_per_prior)
+    ) +
+    coord_cartesian(expand = FALSE) # xlim = c(-per_bin, samples_per_prior + per_bin)
 
   facet_args[["facets"]] <- ~ Parameter
   graph +
     do.call("facet_wrap", facet_args) +
     bayesplot_theme_get() +
-    xaxis_text(FALSE) +
     yaxis_text(FALSE) +
-    xaxis_title(TRUE) +
     yaxis_title(FALSE) +
-    xaxis_ticks(FALSE) +
     yaxis_ticks(FALSE)
 }
 

From be4ecdfc1f0895c850da1ffa30b401d2f9da1b68 Mon Sep 17 00:00:00 2001
From: jgabry <jgabry@gmail.com>
Date: Fri, 17 May 2019 20:25:06 -0400
Subject: [PATCH 6/7] make CI line colors depend on current color scheme

also rename samples_per_prior to thinned_sample_size in internal code
---
 R/sbc.R | 35 +++++++++++++++++++----------------
 1 file changed, 19 insertions(+), 16 deletions(-)

diff --git a/R/sbc.R b/R/sbc.R
index a5aacc81..9e150638 100644
--- a/R/sbc.R
+++ b/R/sbc.R
@@ -68,6 +68,9 @@
 #' color_scheme_set("purple")
 #' sbc_hist(ranks)
 #'
+#' color_scheme_set("blue")
+#' sbc_hist(ranks, worst = 3, facet_args = list(labeller = ggplot2::label_parsed))
+#'
 sbc_hist <- function(ranks, ...,
                      thin = 4, per_bin = 4, worst = 16,
                      prob = 0.99,
@@ -82,7 +85,7 @@ sbc_hist <- function(ranks, ...,
 
   num_prior_draws <- length(ranks)
   thinner <- seq(from = 1, to = nrow(ranks[[1]]), by = thin)
-  samples_per_prior <- length(thinner)
+  thinned_sample_size <- length(thinner)
   u <- t(sapply(ranks, FUN = function(r) 1 + colSums(r[thinner, , drop = FALSE])))
   if (ncol(ranks[[1]]) == 1) {
     u <- t(u)
@@ -91,7 +94,7 @@ sbc_hist <- function(ranks, ...,
 
   if (!is.na(worst)) {
     # order starting with worst (least uniform)
-    kl <- apply(u, 2, function(v) kl_uniform(v, num_prior_draws, samples_per_prior))
+    kl <- apply(u, 2, function(v) kl_uniform(v, num_prior_draws, thinned_sample_size))
     filter <- order(-kl)[1:min(worst, ncol(u))]
     u <- u[, filter, drop = FALSE]
   }
@@ -102,10 +105,10 @@ sbc_hist <- function(ranks, ...,
     u = c(u)
   )
 
-  num_bins <- samples_per_prior / per_bin
-  if (samples_per_prior %% per_bin != 0) {
+  num_bins <- thinned_sample_size / per_bin
+  if (thinned_sample_size %% per_bin != 0) {
     warning(paste0("per_bin (", per_bin, ") does not evenly divide the ",
-                   "number of samples per prior (", samples_per_prior ,")."))
+                   "number of samples per prior (", thinned_sample_size ,")."))
   }
 
   # data for polygon showing expected behavior under uniformity
@@ -118,8 +121,8 @@ sbc_hist <- function(ranks, ...,
   CI <- CI + c(-.5, 0, .5)
   offset <- 2 * per_bin
   polygon_data <- data.frame(
-    x= c(-offset, 0, -offset, samples_per_prior + offset,
-         samples_per_prior, samples_per_prior + offset, -offset),
+    x= c(-offset, 0, -offset, thinned_sample_size + offset,
+         thinned_sample_size, thinned_sample_size + offset, -offset),
     y = c(CI[1], CI[2], CI[3], CI[3], CI[2], CI[1], CI[1])
   )
 
@@ -133,12 +136,12 @@ sbc_hist <- function(ranks, ...,
     ) +
     geom_segment(
       x = -offset,
-      xend = samples_per_prior + offset,
+      xend = thinned_sample_size + offset,
       y = CI[3],
       yend = CI[3],
-      color = "darkgray",
-      alpha = 0.25,
-      size = 0.1
+      color = get_color("mid_highlight"),
+      alpha = 0.5,
+      size = 0.2
     ) +
     geom_histogram(
       bins = num_bins,
@@ -149,18 +152,18 @@ sbc_hist <- function(ranks, ...,
     ) +
     geom_segment(
       x = -offset,
-      xend = samples_per_prior + offset,
+      xend = thinned_sample_size + offset,
       y = CI[1],
       yend = CI[1],
-      color = "darkgray",
+      color = get_color("mid_highlight"),
       alpha = 0.1,
-      size = 0.1
+      size = 0.2
     ) +
     scale_x_continuous(
       name = "Rank statistic",
-      breaks = c(0, round(samples_per_prior / 2), samples_per_prior)
+      breaks = c(0, round(thinned_sample_size / 2), thinned_sample_size)
     ) +
-    coord_cartesian(expand = FALSE) # xlim = c(-per_bin, samples_per_prior + per_bin)
+    coord_cartesian(expand = FALSE) # xlim = c(-per_bin, thinned_sample_size + per_bin)
 
   facet_args[["facets"]] <- ~ Parameter
   graph +

From 633bd0b9002be733242061a37da86f5b8ff9f870 Mon Sep 17 00:00:00 2001
From: jgabry <jgabry@gmail.com>
Date: Sat, 18 May 2019 14:29:20 -0400
Subject: [PATCH 7/7] a few tests for sbc_hist

---
 R/sbc.R                   | 20 ++++++++++--------
 man/sbc_hist.Rd           |  2 ++
 tests/testthat/test-sbc.R | 44 +++++++++++++++++++++++++++++++++++++++
 3 files changed, 57 insertions(+), 9 deletions(-)
 create mode 100644 tests/testthat/test-sbc.R

diff --git a/R/sbc.R b/R/sbc.R
index 9e150638..65558848 100644
--- a/R/sbc.R
+++ b/R/sbc.R
@@ -67,12 +67,14 @@
 #'
 #' color_scheme_set("purple")
 #' sbc_hist(ranks)
-#'
-#' color_scheme_set("blue")
+#' sbc_hist(ranks, worst = NA) # uses original parameter ordering
 #' sbc_hist(ranks, worst = 3, facet_args = list(labeller = ggplot2::label_parsed))
 #'
-sbc_hist <- function(ranks, ...,
-                     thin = 4, per_bin = 4, worst = 16,
+sbc_hist <- function(ranks,
+                     ...,
+                     thin = 4,
+                     per_bin = 4,
+                     worst = 16,
                      prob = 0.99,
                      facet_args = list()) {
   check_ignored_arguments(...)
@@ -118,7 +120,7 @@ sbc_hist <- function(ranks, ...,
     size = num_prior_draws,
     prob = 1 / num_bins
   )
-  CI <- CI + c(-.5, 0, .5)
+  CI <- CI + c(-0.5, 0, 0.5)
   offset <- 2 * per_bin
   polygon_data <- data.frame(
     x= c(-offset, 0, -offset, thinned_sample_size + offset,
@@ -147,7 +149,7 @@ sbc_hist <- function(ranks, ...,
       bins = num_bins,
       fill = get_color("mid"),
       color = get_color("mid_highlight"),
-      size = .25,
+      size = 0.25,
       na.rm = TRUE
     ) +
     geom_segment(
@@ -163,7 +165,7 @@ sbc_hist <- function(ranks, ...,
       name = "Rank statistic",
       breaks = c(0, round(thinned_sample_size / 2), thinned_sample_size)
     ) +
-    coord_cartesian(expand = FALSE) # xlim = c(-per_bin, thinned_sample_size + per_bin)
+    coord_cartesian(expand = FALSE)
 
   facet_args[["facets"]] <- ~ Parameter
   graph +
@@ -180,7 +182,7 @@ sbc_hist <- function(ranks, ...,
 kl_uniform <- function(v, num_prior_draws, samples_per_prior) {
   # D_{KL}(v || uniform)
   # https://en.wikipedia.org/wiki/Kullback%E2%80%93Leibler_divergence
-  expected_pr <- 1.0/samples_per_prior
+  expected_pr <- 1.0 / samples_per_prior
   observed_pr <- table(v) / num_prior_draws
-  sum(observed_pr * log(observed_pr/expected_pr))
+  sum(observed_pr * log(observed_pr / expected_pr))
 }
diff --git a/man/sbc_hist.Rd b/man/sbc_hist.Rd
index 9e6a2515..424534e2 100644
--- a/man/sbc_hist.Rd
+++ b/man/sbc_hist.Rd
@@ -81,6 +81,8 @@ for (n in 1:n_replications) {
 
 color_scheme_set("purple")
 sbc_hist(ranks)
+sbc_hist(ranks, worst = NA) # uses original parameter ordering
+sbc_hist(ranks, worst = 3, facet_args = list(labeller = ggplot2::label_parsed))
 
 }
 \references{
diff --git a/tests/testthat/test-sbc.R b/tests/testthat/test-sbc.R
new file mode 100644
index 00000000..ccd27965
--- /dev/null
+++ b/tests/testthat/test-sbc.R
@@ -0,0 +1,44 @@
+library(bayesplot)
+context("SBC")
+
+set.seed(19)
+pars <- paste0("beta[", 1:2, "]")
+samples_per_prior <- 511
+n_replications <- 50
+ranks <- list()
+for (n in 1:n_replications) {
+  r1 <- matrix(0, nrow=samples_per_prior, ncol=length(pars), dimnames=list(NULL, pars))
+  for (p1 in 1:length(pars)) {
+    r1[sample.int(samples_per_prior, floor(runif(1, 0, samples_per_prior))), p1] <- 1
+  }
+  ranks[[n]] <- r1
+}
+set.seed(NULL)
+
+test_that("sbc_hist returns a ggplot object", {
+  expect_gg(sbc_hist(ranks, worst = NA))
+
+  g <- sbc_hist(ranks)
+  expect_gg(g)
+  expect_equal(nlevels(g$data$Parameter), ncol(ranks[[1]]))
+
+  g2 <- sbc_hist(ranks, worst = 1)
+  expect_gg(g2)
+  expect_equal(nlevels(g2$data$Parameter), 1)
+})
+
+test_that("sbc_hist throws correct warnings and errors", {
+  expect_warning(sbc_hist(ranks, per_bin = 27),
+                 "does not evenly divide the number of samples per prior")
+  expect_warning(sbc_hist(ranks, thin = 3),
+                 "does not evenly divide the number of samples per prior")
+
+  ranks2 <- ranks
+  ranks2[[1]] <- ranks[[1]][, 1, drop=FALSE]
+  expect_error(sbc_hist(ranks2), "Not all matrices in 'ranks' have the same dimensions")
+
+  ranks2[[1]] <- ranks[[1]][, 1, drop=TRUE]
+  expect_error(sbc_hist(ranks2), "all(sapply(ranks, is.matrix)) is not TRUE", fixed = TRUE)
+
+  expect_error(sbc_hist(ranks[[1]]), "is.list(ranks) is not TRUE", fixed=TRUE)
+})