A Quadrature Grid Missing the Posterior Mode

#| '!! 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)

ui <- fluidPage(
    titlePanel("Shinylive App for a Quadrature Grid Missing the Posterior Mode"),
    sidebarLayout(
        sidebarPanel(
            width = 4,
            numericInput("mu_hat", "mu (data location)", value = 8, step = 1),
            sliderInput("n_pts", "number of grid points in u (over 0 to 1)",
                        min = 2, max = 25, value = 9, step = 1),
            helpText("The grid is a midpoint rule, uniform in u over (0, 1); ",
                     "each point's mu is recovered as logit(u). Try mu = 0 ",
                     "first (the grid sits right under the mass, and the two ",
                     "legend numbers agree closely), then move mu up to 8 or ",
                     "20: the quadrature estimate drifts far from the true ",
                     "value however many grid points you add.")
        ),
        mainPanel(
            width = 8,
            plotOutput("squeeze", height = "560px")
        )
    )
)

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

    output$squeeze <- renderPlot({
        n_eff <- 3                                  # likelihood as if 3 obs, sigma = 1
        L     <- function(mu) dnorm(input$mu_hat, mu, 1 / sqrt(n_eff))
        prior <- function(mu) dnorm(mu, 0, 5)
        z     <- function(mu) L(mu) * prior(mu)

        log_sum_exp <- function(lx) { m <- max(lx); m + log(sum(exp(lx - m))) }

        ### this z(mu) is a normal-normal conjugate kernel, so the true log
        ### marginal likelihood is available directly by numerical integration
        log_marlik <- log(integrate(z, -Inf, Inf)$value)

        lo <- min(-8, input$mu_hat - 6)
        hi <- max(14, input$mu_hat + 6)
        mu_grid <- seq(lo, hi, length.out = 500)

        u_pts  <- (seq_len(input$n_pts) - 0.5) / input$n_pts   # midpoint rule, uniform on (0, 1)
        mu_pts <- log(u_pts / (1 - u_pts))          # mu = logit(u)

        ### the midpoint-rule quadrature estimate actually built from this grid:
        ### integral of z(mu) dmu = integral of z(logit(u)) / (u(1-u)) du
        log_z    <- function(mu) dnorm(input$mu_hat, mu, 1 / sqrt(n_eff), log = TRUE) +
            dnorm(mu, 0, 5, log = TRUE)
        h        <- 1 / input$n_pts
        log_marlik_quad <- log_sum_exp(log_z(mu_pts) - log(u_pts * (1 - u_pts))) + log(h)

        par(mar = c(4, 4, 2, 4))
        plot(mu_grid, z(mu_grid), type = "l", lwd = 3, col = "firebrick",
             xlab = expression(mu), ylab = expression(z(mu) == pi(mu) * L(mu)))
        polygon(c(mu_grid, rev(mu_grid)), c(z(mu_grid), rep(0, length(mu_grid))),
                col = adjustcolor("firebrick", 0.15), border = NA)
        segments(mu_pts, 0, mu_pts, z(mu_pts), col = "gray30", lwd = 1.5)
        points(mu_pts, z(mu_pts), pch = 19, col = "gray30")

        par(new = TRUE)
        plot(mu_grid, plogis(mu_grid), type = "l", lwd = 3, col = "steelblue",
             axes = FALSE, xlab = "", ylab = "", ylim = c(0, 1))
        axis(4, col.axis = "steelblue", col.ticks = "steelblue")
        mtext(expression(u == 1 / (1 + e^-mu)), side = 4, line = 2.5, col = "steelblue")

        guide_col <- adjustcolor("gray50", 0.6)
        segments(mu_pts, 0, mu_pts, u_pts, col = guide_col, lty = 3)
        segments(mu_pts, u_pts, max(mu_grid), u_pts, col = guide_col, lty = 3)
        points(mu_pts, u_pts, pch = 19, col = "gray30")

        legend("topleft", bty = "n", lwd = c(3, 3, NA, NA), col = c("firebrick", "steelblue", NA, NA),
               legend = c(expression(z(mu) == pi(mu) * L(mu)),
                          expression(u == 1 / (1 + e^-mu)),
                          as.expression(bquote("true"~log~P(y) == .(round(log_marlik, 3)))),
                          as.expression(bquote("quadrature, this grid"~log~P(y) == .(round(log_marlik_quad, 3))))))
    })
}

shinyApp(ui, server)

About the app

Shows how a fixed quadrature grid can miss the posterior mode entirely, so that the computed marginal likelihood is badly wrong. The app models a single unknown mean \(\mu\).

NoteR source for this app

library(shiny)

ui <- fluidPage(
    titlePanel("Shinylive App for a Quadrature Grid Missing the Posterior Mode"),
    sidebarLayout(
        sidebarPanel(
            width = 4,
            numericInput("mu_hat", "mu (data location)", value = 8, step = 1),
            sliderInput("n_pts", "number of grid points in u (over 0 to 1)",
                        min = 2, max = 25, value = 9, step = 1),
            helpText("The grid is a midpoint rule, uniform in u over (0, 1); ",
                     "each point's mu is recovered as logit(u). Try mu = 0 ",
                     "first (the grid sits right under the mass, and the two ",
                     "legend numbers agree closely), then move mu up to 8 or ",
                     "20: the quadrature estimate drifts far from the true ",
                     "value however many grid points you add.")
        ),
        mainPanel(
            width = 8,
            plotOutput("squeeze", height = "560px")
        )
    )
)

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

    output$squeeze <- renderPlot({
        n_eff <- 3                                  # likelihood as if 3 obs, sigma = 1
        L     <- function(mu) dnorm(input$mu_hat, mu, 1 / sqrt(n_eff))
        prior <- function(mu) dnorm(mu, 0, 5)
        z     <- function(mu) L(mu) * prior(mu)

        log_sum_exp <- function(lx) { m <- max(lx); m + log(sum(exp(lx - m))) }

        ### this z(mu) is a normal-normal conjugate kernel, so the true log
        ### marginal likelihood is available directly by numerical integration
        log_marlik <- log(integrate(z, -Inf, Inf)$value)

        lo <- min(-8, input$mu_hat - 6)
        hi <- max(14, input$mu_hat + 6)
        mu_grid <- seq(lo, hi, length.out = 500)

        u_pts  <- (seq_len(input$n_pts) - 0.5) / input$n_pts   # midpoint rule, uniform on (0, 1)
        mu_pts <- log(u_pts / (1 - u_pts))          # mu = logit(u)

        ### the midpoint-rule quadrature estimate actually built from this grid:
        ### integral of z(mu) dmu = integral of z(logit(u)) / (u(1-u)) du
        log_z    <- function(mu) dnorm(input$mu_hat, mu, 1 / sqrt(n_eff), log = TRUE) +
            dnorm(mu, 0, 5, log = TRUE)
        h        <- 1 / input$n_pts
        log_marlik_quad <- log_sum_exp(log_z(mu_pts) - log(u_pts * (1 - u_pts))) + log(h)

        par(mar = c(4, 4, 2, 4))
        plot(mu_grid, z(mu_grid), type = "l", lwd = 3, col = "firebrick",
             xlab = expression(mu), ylab = expression(z(mu) == pi(mu) * L(mu)))
        polygon(c(mu_grid, rev(mu_grid)), c(z(mu_grid), rep(0, length(mu_grid))),
                col = adjustcolor("firebrick", 0.15), border = NA)
        segments(mu_pts, 0, mu_pts, z(mu_pts), col = "gray30", lwd = 1.5)
        points(mu_pts, z(mu_pts), pch = 19, col = "gray30")

        par(new = TRUE)
        plot(mu_grid, plogis(mu_grid), type = "l", lwd = 3, col = "steelblue",
             axes = FALSE, xlab = "", ylab = "", ylim = c(0, 1))
        axis(4, col.axis = "steelblue", col.ticks = "steelblue")
        mtext(expression(u == 1 / (1 + e^-mu)), side = 4, line = 2.5, col = "steelblue")

        guide_col <- adjustcolor("gray50", 0.6)
        segments(mu_pts, 0, mu_pts, u_pts, col = guide_col, lty = 3)
        segments(mu_pts, u_pts, max(mu_grid), u_pts, col = guide_col, lty = 3)
        points(mu_pts, u_pts, pch = 19, col = "gray30")

        legend("topleft", bty = "n", lwd = c(3, 3, NA, NA), col = c("firebrick", "steelblue", NA, NA),
               legend = c(expression(z(mu) == pi(mu) * L(mu)),
                          expression(u == 1 / (1 + e^-mu)),
                          as.expression(bquote("true"~log~P(y) == .(round(log_marlik, 3)))),
                          as.expression(bquote("quadrature, this grid"~log~P(y) == .(round(log_marlik_quad, 3))))))
    })
}

shinyApp(ui, server)

This app accompanies Numerical Quadrature in the book.