Continuous density example

In this notebook, we create the first three visualizations of a sample from a smooth unbounded distribution.

library(ggplot2)
library(ggdist)
library(bayesplot)

source("code/R/helpers.R")
source("code/R/helpers_qdotplot.R")
source("code/R/helpers_kde.R")
source("code/R/helpers_histogram.R")

theme_set(ppc_paper_theme())

set.seed(7656754)

n <- 1e3
nq <- 100
bw_dots <- sqrt(1 / nq)

# Smooth
sd <- .5
x <- rnorm(n, 0, sd)

We start by generating 1000 normally distributed draws, and make a plot of the true density function of the sampling distribution:

xlim <- max(abs(x))
true_density <- geom_line(aes(
  x = seq(-xlim, xlim, length.out = 1e3),
  y = dnorm(seq(
    -xlim, xlim,
    length.out = 1e3
  ), 0, sd)
),
linewidth = 1)

ggplot() + true_density + labs(x = "", y = "")

1 Quantile dotplot

We visualize the sample with a quantile dot plot using 100 quantiles.

qdots <- ggplot() +
  stat_dots(
    aes(x = x),
    quantiles = nq,
    binwidth = bw_dots,
    overflow = "keep",
    height = 1.1 * (max(
      bin_dots(
        x = quantile(x, stats::ppoints(nq, a = 1 / 2), type = 5),
        y = 0,
        binwidth = bw_dots
      )$y + bw_dots
    )),
    fill = paper_colors["mid"],
    color = paper_colors["mid"]
  ) +
  true_density +
  scale_y_continuous(breaks = c(0,0.5,1)) +
  coord_equal() +
  labs(x = "", y = "")

qdots

PIT ECDF with simultaneous 95% confidence bands

To assess the goodness-of-fit of the visualization to the data, we use a graphical goodness-of-fit test with 95% simultaneous confidence bands for the density visualization matching the data.

n_eval <- nq
ecdf_difference_limits <-
  bayesplot:::ecdf_intervals(
    gamma = bayesplot:::adjust_gamma(N = n, K = n_eval, prob = .95),
    N = n,
    K = n_eval
  )

x0 <- 1:n_eval / n_eval

ecdf_qdot <- ggplot(mapping = aes(x = x0)) +
  geom_step(aes(y = ecdf_difference_limits$lower[-1] / n - x0)) +
  geom_step(aes(y = ecdf_difference_limits$upper[-1] / n - x0)) +
  geom_step(aes(
    y = ecdf(pit_qdotplot(x, nq))(x0) - x0
  ), color = paper_colors["mid"], linewidth = 1) +
  scale_y_continuous(breaks = c(-0.05, 0, 0.05)) +
  scale_x_continuous(breaks = c(0, 0.5, 1)) +
  coord_fixed(ratio = 6.5) +
  labs(x = "", y = "")

ecdf_qdot

2 Histogram

To visualize the sample using a histogram, we first need to determine the bin width. We use the

bw_hist <- 2 * IQR(x) / length(x)^(1 / 3)

hist_plot <- ggplot() +
  geom_histogram(
    aes(x = x, y = after_stat(density)),
    color = paper_colors["dark_highlight"],
    fill = paper_colors["mid"],
    binwidth = bw_hist,
    # center = mean(range(x)),
    #fill = "gray",
  ) +
  true_density +
  scale_y_continuous(breaks = c(0,0.5,1)) +
  coord_equal(
    xlim = layer_scales(qdots)$x$range$range,
    ylim = layer_scales(qdots)$y$range$range,
  ) +
  labs(x = "", y = "")

hist_plot

Goodness-of-fit evaluation shows that the visualization seems to agree with the data.

n_eval <- n
ecdf_difference_limits <-
  bayesplot:::ecdf_intervals(
    gamma = bayesplot:::adjust_gamma(N = n, K = n_eval, prob = .95),
    N = n,
    K = n_eval
  )

x0 <- 1:n_eval / n_eval

ecdf_hist <- ggplot(mapping = aes(x = x0)) +
  geom_step(aes(y = ecdf_difference_limits$lower[-1] / n - x0)) +
  geom_step(aes(y = ecdf_difference_limits$upper[-1] / n - x0)) +
  geom_step(aes(
    y = ecdf(pit_from_hist(hist_plot, x, bw = bw_hist, 1))(x0) - x0
  ), color = paper_colors["mid"], linewidth = 1) +
  scale_y_continuous(breaks = c(-0.05, 0, 0.05)) +
  scale_x_continuous(breaks = c(0, 0.5, 1)) +
  coord_fixed(ratio = 6.5) +
  labs(x = "", y = "")

ecdf_hist

3 Density plot

For the kernel density plots, we use the Sheather Jones method for bandwidth selection.

kde_plot <- ggplot() +
  true_density +
  stat_slab(
    aes(x = x, colour = "Unbounded"),
    density = "unbounded",
    fill = NA,
    normalize = "none",
    scale = 1,
    color = paper_colors["mid"],
    linewidth = 1
  ) +
  labs(
    x = "",
    y = "",
    colour = "Method"
  )

kde_plot +
  scale_y_continuous(breaks = c(0,0.5,1)) +
  ppc_paper_theme(26)

Again, the visualization seems to agree with the data.

ecdf_kde_reg <- ggplot(mapping = aes(x = x0)) +
  geom_step(aes(y = (ecdf_difference_limits$lower[-1]) / n)) +
  geom_step(aes(y = (ecdf_difference_limits$upper[-1]) / n)) +
  geom_step(aes(
    y = ecdf(pit_from_densityplot(kde_plot, 2, x, T))(x0)
  ), color = paper_colors["mid"], linewidth = 1) +
  labs(
    x = "PIT",
    y = "ECDF",
    colour = "Method"
  )

ecdf_kde_reg +
  scale_y_continuous(breaks = c(0,0.5,1)) +
  ppc_paper_theme(26)

ecdf_kde <- ggplot(mapping = aes(x = x0)) +
  geom_step(aes(y = (ecdf_difference_limits$lower[-1] - 1:n) / n)) +
  geom_step(aes(y = (ecdf_difference_limits$upper[-1] - 1:n) / n)) +
  geom_step(aes(
    y = ecdf(pit_from_densityplot(kde_plot, 2, x, T))(x0) - x0
  ), color = paper_colors["mid"], linewidth = 1) +
  labs(
    x = "PIT",
    y = "ECDF difference",
    colour = "Method"
  )

ecdf_kde +
  scale_y_continuous(breaks = c(-0.05,0,0.05)) +
  ppc_paper_theme(26)