Sorted Importance Weights
#| '!! shinylive warning !!': |
#| shinylive does not work in self-contained HTML documents.
#| Please set `embed-resources: false` in your metadata.
#| standalone: true
#| viewerHeight: 820
library(shiny)
log_sum_exp <- function(lx) { m <- max(lx); m + log(sum(exp(lx - m))) }
log_lik <- function(x, mu, w) sum(dnorm(x, mu, exp(w / 2), log = TRUE))
log_prior <- function(mu, w, mu_0, sigma_mu, w_0, sigma_w) {
dnorm(mu, mu_0, sigma_mu, log = TRUE) + dnorm(w, w_0, sigma_w, log = TRUE)
}
neg_log_post <- function(theta, x, mu_0, sigma_mu, w_0, sigma_w) {
-log_lik(x, theta[1], theta[2]) -
log_prior(theta[1], theta[2], mu_0, sigma_mu, w_0, sigma_w)
}
log_dmvnorm_batch <- function(Theta, mu, A) {
Tc <- Theta - mu
quad <- colSums((A %*% Tc) * Tc)
0.5 * (-nrow(Theta) * log(2 * pi) + sum(log(svd(A)$d)) - quad)
}
## same grid-quadrature "truth" as log_marlik_grid() above
marlik_truth <- function(x, mu_0, sigma_mu, w_0, sigma_w, n_grid = 400, width = 8) {
n <- length(x)
fit <- nlm(neg_log_post, c(mean(x), log(stats::var(x))), hessian = TRUE,
x = x, mu_0 = mu_0, sigma_mu = sigma_mu, w_0 = w_0, sigma_w = sigma_w)
se <- sqrt(diag(solve(fit$hessian)))
mu_grid <- fit$estimate[1] + seq(-width, width, length.out = n_grid) * se[1]
w_grid <- fit$estimate[2] + seq(-width, width, length.out = n_grid) * se[2]
SS <- sum(x^2) - 2 * mu_grid * sum(x) + n * mu_grid^2
log_lik_grid <- outer(SS, w_grid,
function(ss, w) -0.5 * n * log(2 * pi) - n * w / 2 - 0.5 * exp(-w) * ss)
log_prior_grid <- outer(dnorm(mu_grid, mu_0, sigma_mu, log = TRUE),
dnorm(w_grid, w_0, sigma_w, log = TRUE), "+")
h_mu <- diff(mu_grid[1:2]); h_w <- diff(w_grid[1:2])
log_sum_exp(as.vector(log_lik_grid + log_prior_grid)) + log(h_mu) + log(h_w)
}
## Laplace-shaped proposal: draw, weight, and return the log estimate, its
## standard error (delta method), and the sorted normalized weights
sim_laplace <- function(x, mu_0, sigma_mu, w_0, sigma_w, S) {
fit <- nlm(neg_log_post, c(mean(x), log(stats::var(x))), hessian = TRUE,
x = x, mu_0 = mu_0, sigma_mu = sigma_mu, w_0 = w_0, sigma_w = sigma_w)
mu_hat <- fit$estimate; H <- fit$hessian
L <- t(chol(solve(H)))
Theta <- L %*% matrix(rnorm(2 * S), 2, S) + mu_hat
log_g <- log_dmvnorm_batch(Theta, mu_hat, H)
log_f <- -apply(Theta, 2, neg_log_post, x = x, mu_0 = mu_0, sigma_mu = sigma_mu,
w_0 = w_0, sigma_w = sigma_w)
lw <- log_f - log_g
m <- max(lw); v <- exp(lw - m); wn <- v / sum(v)
list(est = log_sum_exp(lw) - log(S), se = sd(v) / (sqrt(S) * mean(v)),
ess = 1 / sum(wn^2), wn = sort(wn, decreasing = TRUE))
}
## prior proposal (g = pi): same bookkeeping, unfavourable g
sim_prior <- function(x, mu_0, sigma_mu, w_0, sigma_w, S) {
mus <- rnorm(S, mu_0, sigma_mu); ws <- rnorm(S, w_0, sigma_w)
lw <- vapply(seq_len(S), function(i) log_lik(x, mus[i], ws[i]), numeric(1))
m <- max(lw); v <- exp(lw - m); wn <- v / sum(v)
list(est = log_sum_exp(lw) - log(S), se = sd(v) / (sqrt(S) * mean(v)),
ess = 1 / sum(wn^2), wn = sort(wn, decreasing = TRUE))
}
ui <- fluidPage(
titlePanel("Shinylive App for Sorted Importance Weights"),
sidebarLayout(
sidebarPanel(
width = 3,
sliderInput("n", "data size n", min = 5, max = 200, value = 20, step = 5),
sliderInput("sig_mu", HTML("prior sd on μ"), min = 1, max = 10, value = 5, step = 0.5),
sliderInput("sig_w", HTML("prior sd on w"), min = 1, max = 10, value = 5, step = 0.5),
selectInput("S", "importance sample size S",
choices = c(500, 2000, 5000, 20000), selected = 2000),
actionButton("new", "New sample"),
checkboxInput("logw", "log scale for weights", TRUE)
),
mainPanel(width = 9, plotOutput("panel", height = "720px"))
)
)
server <- function(input, output) {
sim <- reactive({
n <- input$n; S <- as.integer(input$S)
sig_mu <- input$sig_mu; sig_w <- input$sig_w
set.seed(input$new + 1)
x <- rnorm(n)
list(S = S,
truth = marlik_truth(x, 0, sig_mu, 0, sig_w),
lap = sim_laplace(x, 0, sig_mu, 0, sig_w, S),
pri = sim_prior(x, 0, sig_mu, 0, sig_w, S))
})
output$panel <- renderPlot({
s <- sim(); S <- s$S
rank_frac <- seq_len(S) / S
uselog <- input$logw
lwr_lap <- s$lap$est - 1.96 * s$lap$se; upr_lap <- s$lap$est + 1.96 * s$lap$se
lwr_pri <- s$pri$est - 1.96 * s$pri$se; upr_pri <- s$pri$est + 1.96 * s$pri$se
cov_lap <- if (lwr_lap <= s$truth && s$truth <= upr_lap) "darkgreen" else "firebrick"
cov_pri <- if (lwr_pri <= s$truth && s$truth <= upr_pri) "darkgreen" else "firebrick"
draw_panel <- function(wn, col, title, xlab, est, lwr, upr, ess_pct, cov_col) {
if (uselog) {
floor <- max(min(wn[wn > 0]), 1e-300)
wn <- pmax(wn, floor)
ylim <- range(wn)
ylab <- "normalized weight (log scale)"
} else {
ylim <- c(0, max(wn))
ylab <- "normalized weight"
}
plot(rank_frac, wn, type = "n", log = if (uselog) "y" else "",
ylim = ylim, xlab = xlab, ylab = ylab, main = title)
lines(rank_frac, wn, lwd = 3, col = col)
abline(h = 1 / S, lty = 2, col = "darkgreen", lwd = 2)
legend("topright", inset = c(0.02, 0.05), bty = "n", cex = 1.05,
text.col = c("black", cov_col),
legend = c(
sprintf("truth: log P(y) = %.3f", s$truth),
sprintf("est=%.3f CI=(%.3f,%.3f) ESS=%.2f%%", est, lwr, upr, ess_pct)))
}
layout(matrix(1:2, 2))
par(mar = c(2.5, 5, 3, 1))
draw_panel(s$lap$wn, "steelblue", "Laplace-shaped proposal", "",
s$lap$est, lwr_lap, upr_lap, 100 * s$lap$ess / S, cov_lap)
par(mar = c(4.5, 5, 2.5, 1))
draw_panel(s$pri$wn, "firebrick", "prior proposal",
"rank / S (draws sorted by weight, largest first)",
s$pri$est, lwr_pri, upr_pri, 100 * s$pri$ess / S, cov_pri)
})
}
shinyApp(ui, server)About the app
Shows how the distribution of the importance weights reveals the quality of a proposal for estimating a marginal likelihood, comparing a Laplace-shaped proposal with the prior. Sliders control the simulated data size \(n\) and the prior spreads on \(\mu\) and \(w\); the top panel plots the sorted, normalized importance weights from the Laplace-shaped proposal and the bottom panel the same for the prior — each on its own axis, since one proposal’s weights routinely span far more orders of magnitude than the other’s — with a checkbox toggling both between log and linear scale, and the dashed line at \(1/S\) marking perfectly balanced weights in both.
A single number like \(S_{\text{eff}}\) summarizes the weight distribution; looking at the weights themselves is more direct. The app above draws \(S\) importance samples from each proposal and plots their normalized weights, sorted from largest to smallest against rank\(/S\), in two separate panels, because the two proposals’ weights live on such different scales that a shared axis would flatten the well-behaved one to a barely visible line. The weight axis defaults to a log scale — a proposal shaped like the target gives a nearly flat curve close to \(1/S\) (every draw contributes about equally), while a poorly matched proposal gives a curve that falls off a cliff, with a handful of draws near rank \(0\) carrying nearly all the weight; a checkbox switches both panels to a linear scale instead, where that same cliff collapses to a single spike against an otherwise invisible floor — a useful reminder that the log scale is doing real work, not just spreading the picture out. Each panel’s legend reports the grid-quadrature truth and that proposal’s estimate and 95% CI, colored green when the interval covers the truth and red when it does not.
Raise \(n\) with the sliders at their defaults and watch the prior’s curve rotate from a gentle slope into a cliff at the very first few ranks, while the Laplace curve barely moves — the same collapse as in Table 11.2, now visible in the shape of the curve rather than only in a single summary number.
NoteR source for this app
library(shiny)
log_sum_exp <- function(lx) { m <- max(lx); m + log(sum(exp(lx - m))) }
log_lik <- function(x, mu, w) sum(dnorm(x, mu, exp(w / 2), log = TRUE))
log_prior <- function(mu, w, mu_0, sigma_mu, w_0, sigma_w) {
dnorm(mu, mu_0, sigma_mu, log = TRUE) + dnorm(w, w_0, sigma_w, log = TRUE)
}
neg_log_post <- function(theta, x, mu_0, sigma_mu, w_0, sigma_w) {
-log_lik(x, theta[1], theta[2]) -
log_prior(theta[1], theta[2], mu_0, sigma_mu, w_0, sigma_w)
}
log_dmvnorm_batch <- function(Theta, mu, A) {
Tc <- Theta - mu
quad <- colSums((A %*% Tc) * Tc)
0.5 * (-nrow(Theta) * log(2 * pi) + sum(log(svd(A)$d)) - quad)
}
## same grid-quadrature "truth" as log_marlik_grid() above
marlik_truth <- function(x, mu_0, sigma_mu, w_0, sigma_w, n_grid = 400, width = 8) {
n <- length(x)
fit <- nlm(neg_log_post, c(mean(x), log(stats::var(x))), hessian = TRUE,
x = x, mu_0 = mu_0, sigma_mu = sigma_mu, w_0 = w_0, sigma_w = sigma_w)
se <- sqrt(diag(solve(fit$hessian)))
mu_grid <- fit$estimate[1] + seq(-width, width, length.out = n_grid) * se[1]
w_grid <- fit$estimate[2] + seq(-width, width, length.out = n_grid) * se[2]
SS <- sum(x^2) - 2 * mu_grid * sum(x) + n * mu_grid^2
log_lik_grid <- outer(SS, w_grid,
function(ss, w) -0.5 * n * log(2 * pi) - n * w / 2 - 0.5 * exp(-w) * ss)
log_prior_grid <- outer(dnorm(mu_grid, mu_0, sigma_mu, log = TRUE),
dnorm(w_grid, w_0, sigma_w, log = TRUE), "+")
h_mu <- diff(mu_grid[1:2]); h_w <- diff(w_grid[1:2])
log_sum_exp(as.vector(log_lik_grid + log_prior_grid)) + log(h_mu) + log(h_w)
}
## Laplace-shaped proposal: draw, weight, and return the log estimate, its
## standard error (delta method), and the sorted normalized weights
sim_laplace <- function(x, mu_0, sigma_mu, w_0, sigma_w, S) {
fit <- nlm(neg_log_post, c(mean(x), log(stats::var(x))), hessian = TRUE,
x = x, mu_0 = mu_0, sigma_mu = sigma_mu, w_0 = w_0, sigma_w = sigma_w)
mu_hat <- fit$estimate; H <- fit$hessian
L <- t(chol(solve(H)))
Theta <- L %*% matrix(rnorm(2 * S), 2, S) + mu_hat
log_g <- log_dmvnorm_batch(Theta, mu_hat, H)
log_f <- -apply(Theta, 2, neg_log_post, x = x, mu_0 = mu_0, sigma_mu = sigma_mu,
w_0 = w_0, sigma_w = sigma_w)
lw <- log_f - log_g
m <- max(lw); v <- exp(lw - m); wn <- v / sum(v)
list(est = log_sum_exp(lw) - log(S), se = sd(v) / (sqrt(S) * mean(v)),
ess = 1 / sum(wn^2), wn = sort(wn, decreasing = TRUE))
}
## prior proposal (g = pi): same bookkeeping, unfavourable g
sim_prior <- function(x, mu_0, sigma_mu, w_0, sigma_w, S) {
mus <- rnorm(S, mu_0, sigma_mu); ws <- rnorm(S, w_0, sigma_w)
lw <- vapply(seq_len(S), function(i) log_lik(x, mus[i], ws[i]), numeric(1))
m <- max(lw); v <- exp(lw - m); wn <- v / sum(v)
list(est = log_sum_exp(lw) - log(S), se = sd(v) / (sqrt(S) * mean(v)),
ess = 1 / sum(wn^2), wn = sort(wn, decreasing = TRUE))
}
ui <- fluidPage(
titlePanel("Shinylive App for Sorted Importance Weights"),
sidebarLayout(
sidebarPanel(
width = 3,
sliderInput("n", "data size n", min = 5, max = 200, value = 20, step = 5),
sliderInput("sig_mu", HTML("prior sd on μ"), min = 1, max = 10, value = 5, step = 0.5),
sliderInput("sig_w", HTML("prior sd on w"), min = 1, max = 10, value = 5, step = 0.5),
selectInput("S", "importance sample size S",
choices = c(500, 2000, 5000, 20000), selected = 2000),
actionButton("new", "New sample"),
checkboxInput("logw", "log scale for weights", TRUE)
),
mainPanel(width = 9, plotOutput("panel", height = "720px"))
)
)
server <- function(input, output) {
sim <- reactive({
n <- input$n; S <- as.integer(input$S)
sig_mu <- input$sig_mu; sig_w <- input$sig_w
set.seed(input$new + 1)
x <- rnorm(n)
list(S = S,
truth = marlik_truth(x, 0, sig_mu, 0, sig_w),
lap = sim_laplace(x, 0, sig_mu, 0, sig_w, S),
pri = sim_prior(x, 0, sig_mu, 0, sig_w, S))
})
output$panel <- renderPlot({
s <- sim(); S <- s$S
rank_frac <- seq_len(S) / S
uselog <- input$logw
lwr_lap <- s$lap$est - 1.96 * s$lap$se; upr_lap <- s$lap$est + 1.96 * s$lap$se
lwr_pri <- s$pri$est - 1.96 * s$pri$se; upr_pri <- s$pri$est + 1.96 * s$pri$se
cov_lap <- if (lwr_lap <= s$truth && s$truth <= upr_lap) "darkgreen" else "firebrick"
cov_pri <- if (lwr_pri <= s$truth && s$truth <= upr_pri) "darkgreen" else "firebrick"
draw_panel <- function(wn, col, title, xlab, est, lwr, upr, ess_pct, cov_col) {
if (uselog) {
floor <- max(min(wn[wn > 0]), 1e-300)
wn <- pmax(wn, floor)
ylim <- range(wn)
ylab <- "normalized weight (log scale)"
} else {
ylim <- c(0, max(wn))
ylab <- "normalized weight"
}
plot(rank_frac, wn, type = "n", log = if (uselog) "y" else "",
ylim = ylim, xlab = xlab, ylab = ylab, main = title)
lines(rank_frac, wn, lwd = 3, col = col)
abline(h = 1 / S, lty = 2, col = "darkgreen", lwd = 2)
legend("topright", inset = c(0.02, 0.05), bty = "n", cex = 1.05,
text.col = c("black", cov_col),
legend = c(
sprintf("truth: log P(y) = %.3f", s$truth),
sprintf("est=%.3f CI=(%.3f,%.3f) ESS=%.2f%%", est, lwr, upr, ess_pct)))
}
layout(matrix(1:2, 2))
par(mar = c(2.5, 5, 3, 1))
draw_panel(s$lap$wn, "steelblue", "Laplace-shaped proposal", "",
s$lap$est, lwr_lap, upr_lap, 100 * s$lap$ess / S, cov_lap)
par(mar = c(4.5, 5, 2.5, 1))
draw_panel(s$pri$wn, "firebrick", "prior proposal",
"rank / S (draws sorted by weight, largest first)",
s$pri$est, lwr_pri, upr_pri, 100 * s$pri$ess / S, cov_pri)
})
}
shinyApp(ui, server)This app accompanies Importance Sampling in the book.