The Failure of Importance Sampling

Author

Longhai Li

Published

October 6, 2026

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

d <- 5

log_f <- function(th, d) {
    a <- dnorm(th, -d, log = TRUE); b <- dnorm(th, d, log = TRUE)
    m <- pmax(a, b)
    m + log(exp(a - m) + exp(b - m)) - log(2)
}

log_Ew2 <- function(d, sg) {
    A <- 2 - 1 / sg^2
    if (A <= 0) return(Inf)
    lterm <- function(a, b) {
        m <- (a + b) / 2
        B <- 2 * m + d / sg^2
        C <- 2 * m^2 - d^2 / sg^2
        dnorm(a, b, sqrt(2), log = TRUE) + log(sg) + 0.5 * log(2) +
            0.5 * log(2 * pi / A) - 0.5 * (C - B^2 / A)
    }
    lv <- c(lterm(-d, -d), log(2) + lterm(-d, d), lterm(d, d))
    m  <- max(lv)
    m + log(sum(exp(lv - m))) - log(4)
}

ui <- fluidPage(
    titlePanel("Shinylive App for the Failure of Importance Sampling"),
    sidebarLayout(
        sidebarPanel(
            width = 3,
            sliderInput("sg", HTML("proposal sd &sigma;"), min = 0.25, max = 8,
                        value = 1, step = 0.05),
            selectInput("S", "sample size  S",
                        choices = c(500, 2000, 10000), selected = 2000),
            numericInput("seed", "seed", value = 11, min = 1, step = 1),
            actionButton("new", "New sample"),
            tags$hr(style = "margin:8px 0"),
            htmlOutput("readout")
        ),
        mainPanel(width = 9, plotOutput("panels", height = "660px"))
    )
)

server <- function(input, output) {

    sim <- reactive({
        S  <- as.integer(input$S); sg <- input$sg
        set.seed(input$seed + 1000 * input$new)

        th   <- rnorm(S, -d, sg)
        logw <- log_f(th, d) - dnorm(th, -d, sg, log = TRUE)
        mx   <- max(logw); v <- exp(logw - mx)
        wn   <- v / sum(v)

        list(th = th, wn = wn, S = S, sg = sg,
             Chat = mean(v) * exp(mx),
             se   = sd(v) / sqrt(S) * exp(mx),
             ess  = 1 / sum(wn^2),
             Ehat = sum(wn * th),
             lsum = mx + log(sum(v)),
             ess_x = 100 * exp(-log_Ew2(d, sg)))
    })

    output$panels <- renderPlot({
        s  <- sim(); th <- s$th; sg <- s$sg
        xr <- range(c(th, -d - 4, d + 4))
        xs <- seq(xr[1], xr[2], length.out = 1000)
        fx <- exp(log_f(xs, d)); gx <- dnorm(xs, -d, sg)
        lwc <- log_f(xs, d) - dnorm(xs, -d, sg, log = TRUE) - s$lsum

        layout(matrix(1:2, 2), heights = c(3, 2))

        par(mar = c(0.4, 5, 3, 1))
        yl <- c(1e-10, max(fx, gx) * 2)
        plot(xs, fx, type = "n", log = "y", ylim = yl, xlim = xr, xaxt = "n",
             xlab = "", ylab = "density (log scale)",
             main = "bimodal target and a proposal centred on one mode")
        lines(xs, fx, lwd = 4, col = "red")
        lines(xs, gx, lwd = 3, col = "grey60")
        rug(th, col = adjustcolor("grey45", 0.35))
        legend("bottomright", bty = "n", cex = 1.05, text.col = "navy",
               lwd = c(4, 3, NA, NA, NA),
               col = c("grey25", "grey60", NA, NA, NA),
               legend = TeX(c(
                   r"(target $f$)",
                   sprintf(r"(proposal $g$, sd $= %.2f$)", sg),
                   sprintf(r"($C_f/C_g$: true 1, est. %.3f)", s$Chat),
                   sprintf(r"(95%% CI $(%.3f, %.3f)$)", s$Chat - 1.96 * s$se, s$Chat + 1.96 * s$se),
                   sprintf(r"($E_f(\theta)$: true 0, est. %.2f)", s$Ehat))))

        par(mar = c(4.5, 5, 0.4, 1))
        wn <- s$wn
        wl <- c(min(wn) / 5, max(wn) * 5)
        plot(th, wn, type = "n", log = "y", ylim = wl, xlim = xr,
             xlab = expression(theta), ylab = "normalized weight (log scale)")
        lines(xs, exp(pmin(lwc, 700)), lwd = 2, col = "grey25")
        segments(th, wl[1], th, wn, col = adjustcolor("grey45", 0.5))
        abline(h = 1 / s$S, lty = 2, col = "grey55", lwd = 2)
        legend("bottomright", bty = "n", cex = 1.05, text.col = "navy",
               lty = c(2, 1, NA, NA), lwd = c(2, 2, NA, NA),
               col = c("grey55", "grey25", NA, NA),
               legend = TeX(c(
                   r"($1/S$)",
                   r"(normalized weight function)",
                   sprintf(r"(observed $S_{eff}$ = %.1f%% of S)", 100 * s$ess / s$S),
                   sprintf(r"(exact $S_{eff}$ = %.3g%% of S)", s$ess_x))))
    })

    output$readout <- renderUI({
        s <- sim()
        HTML(sprintf(
            "<b>C<sub>f</sub>/C<sub>g</sub></b><br>true 1, est. %.4f<br>
             <b>95%% CI</b> (%.4f, %.4f)<br>
             <b>E<sub>f</sub>(theta)</b><br>true 0, est. %.3f<br>
             <b>observed ESS</b> %.1f%%<br>
             <b>exact ESS</b> %.3g%%<br>
             <b>Var<sub>f</sub>(w)</b> %s",
            s$Chat, s$Chat - 1.96 * s$se, s$Chat + 1.96 * s$se, s$Ehat,
            100 * s$ess / s$S, s$ess_x,
            if (s$sg <= 1 / sqrt(2)) "infinite" else "finite"))
    })
}

shinyApp(ui, server)

About the app

Illustrates how importance sampling can fail without warning when the proposal is too narrow to cover the target. The proposal’s standard deviation \(\sigma\) is under slider control and both panels are redrawn for each setting. For \(\sigma \le 1/\sqrt{2}\) the weight variance is infinite; near \(\sigma = 1\) the failure is silent, with the observed effective sample size close to \(S\) while the estimate of \(C_f/C_g\) sits at one half; past \(\sigma \approx 4\) the proposal reaches the second mode, the observed effective sample size falls, and the estimate moves onto the truth.

the app above varies \(\sigma\) with everything else held fixed. Three regimes appear. Below \(\sigma=1/\sqrt{2}\) the variance is infinite and the exact effective fraction is reported as zero. Between there and roughly \(\sigma=3\) the failure is silent in the sense above: observed \(S_{\text{eff}}\) near \(S\), estimate near \(\tfrac12\). Past \(\sigma\approx4\) the proposal is over-dispersed enough to place draws near \(+d\); the observed \(S_{\text{eff}}\) then drops, which looks like deterioration and is in fact the diagnostic finally working, and the estimate moves onto the truth. The cure for a missed mode is a proposal that is too diffuse rather than too well matched.

library(shiny)
library(latex2exp)

d <- 5

log_f <- function(th, d) {
    a <- dnorm(th, -d, log = TRUE); b <- dnorm(th, d, log = TRUE)
    m <- pmax(a, b)
    m + log(exp(a - m) + exp(b - m)) - log(2)
}

log_Ew2 <- function(d, sg) {
    A <- 2 - 1 / sg^2
    if (A <= 0) return(Inf)
    lterm <- function(a, b) {
        m <- (a + b) / 2
        B <- 2 * m + d / sg^2
        C <- 2 * m^2 - d^2 / sg^2
        dnorm(a, b, sqrt(2), log = TRUE) + log(sg) + 0.5 * log(2) +
            0.5 * log(2 * pi / A) - 0.5 * (C - B^2 / A)
    }
    lv <- c(lterm(-d, -d), log(2) + lterm(-d, d), lterm(d, d))
    m  <- max(lv)
    m + log(sum(exp(lv - m))) - log(4)
}

ui <- fluidPage(
    titlePanel("Shinylive App for the Failure of Importance Sampling"),
    sidebarLayout(
        sidebarPanel(
            width = 3,
            sliderInput("sg", HTML("proposal sd &sigma;"), min = 0.25, max = 8,
                        value = 1, step = 0.05),
            selectInput("S", "sample size  S",
                        choices = c(500, 2000, 10000), selected = 2000),
            numericInput("seed", "seed", value = 11, min = 1, step = 1),
            actionButton("new", "New sample"),
            tags$hr(style = "margin:8px 0"),
            htmlOutput("readout")
        ),
        mainPanel(width = 9, plotOutput("panels", height = "660px"))
    )
)

server <- function(input, output) {

    sim <- reactive({
        S  <- as.integer(input$S); sg <- input$sg
        set.seed(input$seed + 1000 * input$new)

        th   <- rnorm(S, -d, sg)
        logw <- log_f(th, d) - dnorm(th, -d, sg, log = TRUE)
        mx   <- max(logw); v <- exp(logw - mx)
        wn   <- v / sum(v)

        list(th = th, wn = wn, S = S, sg = sg,
             Chat = mean(v) * exp(mx),
             se   = sd(v) / sqrt(S) * exp(mx),
             ess  = 1 / sum(wn^2),
             Ehat = sum(wn * th),
             lsum = mx + log(sum(v)),
             ess_x = 100 * exp(-log_Ew2(d, sg)))
    })

    output$panels <- renderPlot({
        s  <- sim(); th <- s$th; sg <- s$sg
        xr <- range(c(th, -d - 4, d + 4))
        xs <- seq(xr[1], xr[2], length.out = 1000)
        fx <- exp(log_f(xs, d)); gx <- dnorm(xs, -d, sg)
        lwc <- log_f(xs, d) - dnorm(xs, -d, sg, log = TRUE) - s$lsum

        layout(matrix(1:2, 2), heights = c(3, 2))

        par(mar = c(0.4, 5, 3, 1))
        yl <- c(1e-10, max(fx, gx) * 2)
        plot(xs, fx, type = "n", log = "y", ylim = yl, xlim = xr, xaxt = "n",
             xlab = "", ylab = "density (log scale)",
             main = "bimodal target and a proposal centred on one mode")
        lines(xs, fx, lwd = 4, col = "red")
        lines(xs, gx, lwd = 3, col = "grey60")
        rug(th, col = adjustcolor("grey45", 0.35))
        legend("bottomright", bty = "n", cex = 1.05, text.col = "navy",
               lwd = c(4, 3, NA, NA, NA),
               col = c("grey25", "grey60", NA, NA, NA),
               legend = TeX(c(
                   r"(target $f$)",
                   sprintf(r"(proposal $g$, sd $= %.2f$)", sg),
                   sprintf(r"($C_f/C_g$: true 1, est. %.3f)", s$Chat),
                   sprintf(r"(95%% CI $(%.3f, %.3f)$)", s$Chat - 1.96 * s$se, s$Chat + 1.96 * s$se),
                   sprintf(r"($E_f(\theta)$: true 0, est. %.2f)", s$Ehat))))

        par(mar = c(4.5, 5, 0.4, 1))
        wn <- s$wn
        wl <- c(min(wn) / 5, max(wn) * 5)
        plot(th, wn, type = "n", log = "y", ylim = wl, xlim = xr,
             xlab = expression(theta), ylab = "normalized weight (log scale)")
        lines(xs, exp(pmin(lwc, 700)), lwd = 2, col = "grey25")
        segments(th, wl[1], th, wn, col = adjustcolor("grey45", 0.5))
        abline(h = 1 / s$S, lty = 2, col = "grey55", lwd = 2)
        legend("bottomright", bty = "n", cex = 1.05, text.col = "navy",
               lty = c(2, 1, NA, NA), lwd = c(2, 2, NA, NA),
               col = c("grey55", "grey25", NA, NA),
               legend = TeX(c(
                   r"($1/S$)",
                   r"(normalized weight function)",
                   sprintf(r"(observed $S_{eff}$ = %.1f%% of S)", 100 * s$ess / s$S),
                   sprintf(r"(exact $S_{eff}$ = %.3g%% of S)", s$ess_x))))
    })

    output$readout <- renderUI({
        s <- sim()
        HTML(sprintf(
            "<b>C<sub>f</sub>/C<sub>g</sub></b><br>true 1, est. %.4f<br>
             <b>95%% CI</b> (%.4f, %.4f)<br>
             <b>E<sub>f</sub>(theta)</b><br>true 0, est. %.3f<br>
             <b>observed ESS</b> %.1f%%<br>
             <b>exact ESS</b> %.3g%%<br>
             <b>Var<sub>f</sub>(w)</b> %s",
            s$Chat, s$Chat - 1.96 * s$se, s$Chat + 1.96 * s$se, s$Ehat,
            100 * s$ess / s$S, s$ess_x,
            if (s$sg <= 1 / sqrt(2)) "infinite" else "finite"))
    })
}

shinyApp(ui, server)

This app accompanies Importance Sampling in the book.