Perfect. I spent some time after dinner to put together a `R`

demonstration for anyone interested.

```
c_dark <- c("#8F2727")
c_dark_highlight <- c("#7C0000")
m <- 0
s <- 0.5
s2 <- s * s
sigma <- 1
sigma2 <- sigma * sigma
N <- 100000 # Number of replications
L <- 25 # Number of posterior samples
B <- 25 # Number of rank bins
# Correct SBC
ranks <- rep(0, 100)
for (n in 1:N) {
theta <- rnorm(1, m, s)
y <- rnorm(1, theta, sigma)
post_m <- (y * s2 + m * sigma2) / (s2 + sigma2)
post_s <- sqrt( s2 * sigma2 / (s2 + sigma2) )
post_theta <- rnorm(L, post_m, post_s)
ranks[n] <- sum(post_theta < theta)
}
sbc_hist <- hist(ranks, breaks=seq(0, L + 1, 1) - 0.5, plot=FALSE)
plot(sbc_hist, main="SBC (Correct)", xlab="Prior Rank", yaxt='n', ylab="")
low <- qbinom(0.005, N, 1 / (L + 1))
mid <- qbinom(0.5, N, 1 / (L + 1))
high <- qbinom(0.995, N, 1 / (L + 1))
bar_x <- c(-1, L + 1, L + 0.5, L + 1, -1, -0.5, -1)
bar_y <- c(high, high, mid, low, low, mid, high)
polygon(bar_x, bar_y, col=c("#DDDDDD"), border=NA)
segments(x0=-0.5, x1=(L + 0.5), y0=mid, y1=mid, col=c("#999999"), lwd=2)
plot(sbc_hist, col=c_dark, border=c_dark_highlight, add=T)
# Incorrect SBC
ranks <- rep(0, 100)
for (n in 1:N) {
theta <- rnorm(1, m, s)
y <- rnorm(L, theta, sigma)
post_m <- (y * s2 + m * sigma2) / (s2 + sigma2)
post_s <- sqrt( s2 * sigma2 / (s2 + sigma2) )
post_theta <- rnorm(L, post_m, post_s)
ranks[n] <- sum(post_theta < theta)
}
sbc_hist <- hist(ranks, breaks=seq(0, L + 1, 1) - 0.5, plot=FALSE)
plot(sbc_hist, main="SBC (Incorrect)", xlab="Prior Rank", yaxt='n', ylab="")
low <- qbinom(0.005, N, 1 / (L + 1))
mid <- qbinom(0.5, N, 1 / (L + 1))
high <- qbinom(0.995, N, 1 / (L + 1))
bar_x <- c(-1, L + 1, L + 0.5, L + 1, -1, -0.5, -1)
bar_y <- c(high, high, mid, low, low, mid, high)
polygon(bar_x, bar_y, col=c("#DDDDDD"), border=NA)
segments(x0=-0.5, x1=(L + 0.5), y0=mid, y1=mid, col=c("#999999"), lwd=2)
plot(sbc_hist, col=c_dark, border=c_dark_highlight, add=T)
```