Envelope Tightness in Rejection Sampling

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

f_demo <- function(x) exp(-x^2 / 2)
g_demo <- function(x, c) c * exp(0.5 - abs(x))
rlaplace <- function(n) { u <- runif(n); ifelse(u < 0.5, log(2 * u), -log(2 * (1 - u))) }

ui <- fluidPage(
    titlePanel("Shinylive App for Envelope Tightness in Rejection Sampling"),
    sidebarLayout(
        sidebarPanel(
            width = 4,
            sliderInput("c", "envelope scale c (g = c * e^(1/2 - |x|))",
                        min = 0.6, max = 3, value = 1.2, step = 0.05),
            helpText("c >= 1 keeps g above f everywhere: a valid, if looser, ",
                     "envelope. c < 1 breaks the envelope condition near ",
                     "|x| = 1, where the target actually exceeds it."),
            actionButton("resample", "draw new points"),
            tableOutput("tab")
        ),
        mainPanel(
            width = 8,
            plotOutput("scatter", height = "440px")
        )
    )
)

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

    pts <- eventReactive(list(input$c, input$resample), {
        set.seed(input$resample + 1)
        N <- 500
        x <- rlaplace(N)
        v <- runif(N) * g_demo(x, input$c)
        list(x = x, v = v, acc = v < f_demo(x))
    })

    output$scatter <- renderPlot({
        p <- pts()
        xs <- seq(-6, 6, length.out = 500)
        par(mar = c(4, 4, 2, 1))
        plot(xs, g_demo(xs, input$c), type = "l", lwd = 3, col = "steelblue",
             ylim = c(0, max(2.2, 1.05 * input$c * exp(0.5))),
             xlab = "x", ylab = "v", main = sprintf("c = %.2f", input$c))
        lines(xs, f_demo(xs), lwd = 3, col = "firebrick")
        bad <- f_demo(xs) > g_demo(xs, input$c)
        if (any(bad)) lines(xs[bad], f_demo(xs)[bad], col = "black", lwd = 5)
        points(p$x[!p$acc], p$v[!p$acc], pch = 4, col = adjustcolor("gray40", 0.6), cex = 0.7)
        points(p$x[p$acc], p$v[p$acc], pch = 19, col = adjustcolor("firebrick", 0.7), cex = 0.7)
        legend("topright", bty = "n", lwd = c(3, 3, 5), col = c("firebrick", "steelblue", "black"),
               legend = c("target f", "envelope g", "f exceeds g here: envelope invalid"))
    })

    output$tab <- renderTable({
        p <- pts()
        data.frame(quantity = c("theoretical acceptance", "observed acceptance", "envelope valid?"),
                   value = c(sprintf("%.3f", 0.7602 / input$c), sprintf("%.3f", mean(p$acc)),
                             ifelse(input$c >= 1, "yes", "NO (c < 1)")))
    }, colnames = FALSE)
}

shinyApp(ui, server)

About the app

Shows how the tightness of the envelope determines the efficiency of rejection sampling: a looser envelope still gives exact draws but rejects more of them. The app above scales the envelope by a factor \(c \ge 1\), \(g_c(x) = c\, e^{1/2-|x|}\): larger \(c\) keeps the envelope valid but wastes more draws (acceptance rate \(= 0.760/c\)).

Dragging \(c\) below \(1\) shows what an invalid envelope looks like — the target pokes through it, and rejection sampling would silently under-represent that region.

NoteReproducing this figure

The chunk above uses the shinylive extension, compiled to WebAssembly so it runs in the reader’s browser with no Shiny server (quarto add quarto-ext/shinylive, then filters: [shinylive] in the document or project YAML). Only base R and shiny are available inside the app.

NoteR source for this app

library(shiny)

f_demo <- function(x) exp(-x^2 / 2)
g_demo <- function(x, c) c * exp(0.5 - abs(x))
rlaplace <- function(n) { u <- runif(n); ifelse(u < 0.5, log(2 * u), -log(2 * (1 - u))) }

ui <- fluidPage(
    titlePanel("Shinylive App for Envelope Tightness in Rejection Sampling"),
    sidebarLayout(
        sidebarPanel(
            width = 4,
            sliderInput("c", "envelope scale c (g = c * e^(1/2 - |x|))",
                        min = 0.6, max = 3, value = 1.2, step = 0.05),
            helpText("c >= 1 keeps g above f everywhere: a valid, if looser, ",
                     "envelope. c < 1 breaks the envelope condition near ",
                     "|x| = 1, where the target actually exceeds it."),
            actionButton("resample", "draw new points"),
            tableOutput("tab")
        ),
        mainPanel(
            width = 8,
            plotOutput("scatter", height = "440px")
        )
    )
)

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

    pts <- eventReactive(list(input$c, input$resample), {
        set.seed(input$resample + 1)
        N <- 500
        x <- rlaplace(N)
        v <- runif(N) * g_demo(x, input$c)
        list(x = x, v = v, acc = v < f_demo(x))
    })

    output$scatter <- renderPlot({
        p <- pts()
        xs <- seq(-6, 6, length.out = 500)
        par(mar = c(4, 4, 2, 1))
        plot(xs, g_demo(xs, input$c), type = "l", lwd = 3, col = "steelblue",
             ylim = c(0, max(2.2, 1.05 * input$c * exp(0.5))),
             xlab = "x", ylab = "v", main = sprintf("c = %.2f", input$c))
        lines(xs, f_demo(xs), lwd = 3, col = "firebrick")
        bad <- f_demo(xs) > g_demo(xs, input$c)
        if (any(bad)) lines(xs[bad], f_demo(xs)[bad], col = "black", lwd = 5)
        points(p$x[!p$acc], p$v[!p$acc], pch = 4, col = adjustcolor("gray40", 0.6), cex = 0.7)
        points(p$x[p$acc], p$v[p$acc], pch = 19, col = adjustcolor("firebrick", 0.7), cex = 0.7)
        legend("topright", bty = "n", lwd = c(3, 3, 5), col = c("firebrick", "steelblue", "black"),
               legend = c("target f", "envelope g", "f exceeds g here: envelope invalid"))
    })

    output$tab <- renderTable({
        p <- pts()
        data.frame(quantity = c("theoretical acceptance", "observed acceptance", "envelope valid?"),
                   value = c(sprintf("%.3f", 0.7602 / input$c), sprintf("%.3f", mean(p$acc)),
                             ifelse(input$c >= 1, "yes", "NO (c < 1)")))
    }, colnames = FALSE)
}

shinyApp(ui, server)

This app accompanies Rejection Sampling in the book.