Rejection Sampling from a Gamma Density with a Cauchy Envelope

#| '!! shinylive warning !!': |
#|   shinylive does not work in self-contained HTML documents.
#|   Please set `embed-resources: false` in your metadata.
#| standalone: true
#| viewerHeight: 610
library(shiny)

log_g_gamma <- function(x, alpha) {
    (alpha - 1) * (log(alpha - 1) - 1) - log(1 + (x - (alpha - 1))^2 / (2 * alpha - 1))
}
g_gamma <- function(x, alpha) exp(log_g_gamma(x, alpha))

ui <- fluidPage(
    titlePanel("Shinylive App for Rejection Sampling from a Gamma Density with a Cauchy Envelope"),
    sidebarLayout(
        sidebarPanel(
            width = 3,
            sliderInput("alpha", "shape alpha (target: Gamma(alpha, 1))",
                        min = 1.1, max = 8, value = 3, step = 0.05),
            helpText("The envelope is only valid (g >= f everywhere) for ",
                     "alpha > 2. Below alpha = 2 the target pokes through ",
                     "the envelope near the mode, and the histogram falls ",
                     "slightly short of the true density there."),
            actionButton("resample", "draw new sample")
        ),
        mainPanel(
            width = 9,
            plotOutput("panels", height = "420px")
        )
    )
)

server <- function(input, output, session) {

    draws <- eventReactive(list(input$alpha, input$resample), {
        set.seed(input$resample + 1)
        a <- input$alpha; n <- 2000
        out <- numeric(0); tries <- 0; max_tries <- 2e6
        while (length(out) < n && tries < max_tries) {
            m <- n - length(out)
            x <- rcauchy(m) * sqrt(2 * a - 1) + (a - 1)
            tries <- tries + m
            keep <- log(runif(m)) < dgamma(x, a, log = TRUE) - log_g_gamma(x, a)
            out <- c(out, x[keep])
        }
        list(x = out[seq_len(min(n, length(out)))], accept = length(out) / tries)
    })

    output$panels <- renderPlot({
        d <- draws(); a <- input$alpha
        hi <- max(15, qgamma(0.999, a) * 1.3)
        xs <- seq(1e-3, hi, length.out = 500)
        ymax <- max(dgamma(xs, a), na.rm = TRUE) * 1.6   # scale to the target's
        # peak, not the envelope's: g touches f at the mode but can be much
        # looser far from it (especially for alpha far above 2), which would
        # otherwise blow up the y-axis and hide the interesting region

        par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))

        plot(xs, g_gamma(xs, a), type = "l", lwd = 3, col = "steelblue",
             ylim = c(0, ymax), xlab = "x", ylab = "density",
             main = sprintf("target vs. envelope (alpha = %.2f)", a))
        lines(xs, dgamma(xs, a), lwd = 3, col = "firebrick")
        bad <- dgamma(xs, a) > g_gamma(xs, a)
        if (any(bad)) lines(xs[bad], dgamma(xs, a)[bad], col = "black", lwd = 5)
        legend("topright", bty = "n", lwd = c(3, 3, 5), col = c("firebrick", "steelblue", "black"),
               legend = c("target f", "envelope g", "f > g: invalid"))

        hist(d$x, breaks = 30, freq = FALSE, col = "steelblue", border = "white",
             xlim = c(0, hi), xlab = "x",
             main = sprintf("%d draws (accept rate %.1f%%)", length(d$x), 100 * d$accept))
        curve(dgamma(x, a), add = TRUE, lwd = 3, col = "firebrick")
        legend("topright", bty = "n", lwd = 3, col = "firebrick", legend = "true Gamma density")
    })
}

shinyApp(ui, server)

About the app

Shows how the efficiency of rejection sampling for the Gamma distribution changes with the shape parameter while the draws remain exact. The app above shows the target \(f(x)=\text{Gamma}(\alpha,1)\) density against the Cauchy-based envelope \(g(x)=\exp(\texttt{log\_g\_gamma}(x,\alpha))\) as \(\alpha\) is dragged (left), and a histogram of 2000 fresh rejection-sampled draws against the true density (right) — an exact match whenever the envelope is valid.

Drag \(\alpha\) below \(2\) and the target pokes through the envelope, as in Table 10.2; the crossing is a modest one for this particular envelope, so the histogram only dips slightly short of the curve in the affected region rather than breaking outright, but the sampler is no longer drawing from the exact target.

NoteR source for this app

library(shiny)

log_g_gamma <- function(x, alpha) {
    (alpha - 1) * (log(alpha - 1) - 1) - log(1 + (x - (alpha - 1))^2 / (2 * alpha - 1))
}
g_gamma <- function(x, alpha) exp(log_g_gamma(x, alpha))

ui <- fluidPage(
    titlePanel("Shinylive App for Rejection Sampling from a Gamma Density with a Cauchy Envelope"),
    sidebarLayout(
        sidebarPanel(
            width = 3,
            sliderInput("alpha", "shape alpha (target: Gamma(alpha, 1))",
                        min = 1.1, max = 8, value = 3, step = 0.05),
            helpText("The envelope is only valid (g >= f everywhere) for ",
                     "alpha > 2. Below alpha = 2 the target pokes through ",
                     "the envelope near the mode, and the histogram falls ",
                     "slightly short of the true density there."),
            actionButton("resample", "draw new sample")
        ),
        mainPanel(
            width = 9,
            plotOutput("panels", height = "420px")
        )
    )
)

server <- function(input, output, session) {

    draws <- eventReactive(list(input$alpha, input$resample), {
        set.seed(input$resample + 1)
        a <- input$alpha; n <- 2000
        out <- numeric(0); tries <- 0; max_tries <- 2e6
        while (length(out) < n && tries < max_tries) {
            m <- n - length(out)
            x <- rcauchy(m) * sqrt(2 * a - 1) + (a - 1)
            tries <- tries + m
            keep <- log(runif(m)) < dgamma(x, a, log = TRUE) - log_g_gamma(x, a)
            out <- c(out, x[keep])
        }
        list(x = out[seq_len(min(n, length(out)))], accept = length(out) / tries)
    })

    output$panels <- renderPlot({
        d <- draws(); a <- input$alpha
        hi <- max(15, qgamma(0.999, a) * 1.3)
        xs <- seq(1e-3, hi, length.out = 500)
        ymax <- max(dgamma(xs, a), na.rm = TRUE) * 1.6   # scale to the target's
        # peak, not the envelope's: g touches f at the mode but can be much
        # looser far from it (especially for alpha far above 2), which would
        # otherwise blow up the y-axis and hide the interesting region

        par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))

        plot(xs, g_gamma(xs, a), type = "l", lwd = 3, col = "steelblue",
             ylim = c(0, ymax), xlab = "x", ylab = "density",
             main = sprintf("target vs. envelope (alpha = %.2f)", a))
        lines(xs, dgamma(xs, a), lwd = 3, col = "firebrick")
        bad <- dgamma(xs, a) > g_gamma(xs, a)
        if (any(bad)) lines(xs[bad], dgamma(xs, a)[bad], col = "black", lwd = 5)
        legend("topright", bty = "n", lwd = c(3, 3, 5), col = c("firebrick", "steelblue", "black"),
               legend = c("target f", "envelope g", "f > g: invalid"))

        hist(d$x, breaks = 30, freq = FALSE, col = "steelblue", border = "white",
             xlim = c(0, hi), xlab = "x",
             main = sprintf("%d draws (accept rate %.1f%%)", length(d$x), 100 * d$accept))
        curve(dgamma(x, a), add = TRUE, lwd = 3, col = "firebrick")
        legend("topright", bty = "n", lwd = 3, col = "firebrick", legend = "true Gamma density")
    })
}

shinyApp(ui, server)

This app accompanies Rejection Sampling in the book.