Need help on sbc uniform proof

General
11

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)
3 Likes