Shinylive app for simple random sampling
Simple random sampling
#| '!! shinylive warning !!': |
#| shinylive does not work in self-contained HTML documents.
#| Please set `embed-resources: false` in your metadata.
#| standalone: true
#| viewerHeight: 780
## file: app.R
library(shiny)
data_url <- "https://raw.githubusercontent.com/longhaiSK/sampling/main/acres92.csv"
download.file(data_url, "acres92.csv")
pop <- read.csv("acres92.csv")$acres92
N <- length(pop)
ybarU <- mean(pop)
S <- sd(pop)
## ---- broken x-axis for the population panel --------------------------------
## The drawn region is [0, 2 * ybarU]: the visible scale [0, cut] plus a
## compressed lane holding the right tail. That total width puts the population
## mean exactly at the midpoint, matching the histogram panel below.
lane_frac <- 0.10
cut <- 2 * ybarU / (1 + lane_frac)
lane_w <- lane_frac * cut
xmax_top <- cut + lane_w # == 2 * ybarU
tail_id <- which(pop > cut)
n_tail <- length(tail_id)
px <- pop # plotting positions
if (n_tail > 0) {
r <- rank(pop[tail_id], ties.method = "first") / (n_tail + 1)
px[tail_id] <- cut + lane_w * (0.10 + 0.80 * r)
}
# fixed jitter band for the rug of raw population values
set.seed(42)
rug_y <- runif(N, -0.42, -0.08)
dens <- density(pop, from = 0, to = cut)
dens_y <- dens$y / max(dens$y)
## ---- fixed window for the sampling distribution ----------------------------
se100 <- sqrt(1 - 100 / N) * S / sqrt(100)
half <- 5 * se100
xlim_b <- ybarU + c(-1, 1) * half
MAR <- c(4.5, 4.5, 3, 1.5) # identical in both panels, so the scales align
## ---- inline SVG button icons (no font or icon-library dependency) ----------
svg_icon <- function(paths)
HTML(sprintf(paste0('<svg width="16" height="16" viewBox="0 0 16 16" ',
'fill="currentColor" style="vertical-align:-2px">%s</svg>'), paths))
icon_play <- svg_icon('<path d="M3 2v12l10-6z"/>')
icon_pause <- svg_icon('<rect x="3" y="2" width="4" height="12"/><rect x="9" y="2" width="4" height="12"/>')
icon_skip <- svg_icon('<path d="M1 2v12l7-6zM8 2v12l6-6z"/><rect x="14" y="2" width="2" height="12"/>')
ui <- fluidPage(
titlePanel("Shinylive App for Simple Random Sampling"),
tags$style(HTML(
".well {padding: 10px 12px; margin-bottom: 10px;}
.form-group {margin-bottom: 6px;}
.irs {margin-bottom: 0;}
.btn {margin-bottom: 6px;}"
)),
wellPanel(
fluidRow(
column(
4,
sliderInput("n", "Sample size (n)", min = 10, max = 1000, value = 300, step = 10),
sliderInput("conf", "Confidence level", min = 0.80, max = 0.99, value = 0.95, step = 0.01)
),
column(
4,
sliderInput("speed", "Samples per second", min = 1, max = 10, value = 3, step = 1),
numericInput("nrep", "Stop after this many samples",
value = 500, min = 10, max = 20000, step = 50)
),
column(
4,
actionButton("draw", "Draw new sample", class = "btn-primary", width = "100%"),
br(), br(),
fluidRow(
column(4, actionButton("play", icon_play, class = "btn-success", width = "100%",
title = "Play")),
column(4, actionButton("pause", icon_pause, class = "btn-warning", width = "100%",
title = "Pause")),
column(4, actionButton("skip", icon_skip, class = "btn-info", width = "100%",
title = "Skip to the end"))
),
br(),
actionButton("reset", "Reset", width = "100%")
)
)
),
plotOutput("mainplot", height = "480px")
)
server <- function(input, output, session) {
rv <- reactiveValues(
idx = NULL, mean = NULL, ci = NULL, covered = NULL,
means = numeric(0), hits = logical(0)
)
running <- reactiveVal(FALSE)
# guard against a blank or absurd entry in the numeric box
nrep <- reactive({
v <- input$nrep
if (is.null(v) || is.na(v) || v < 1) 1L else min(as.integer(v), 20000L)
})
# draw k samples; the last one becomes the "current" sample, and the reactive
# state is written once so a large batch triggers a single redraw
draw_many <- function(k) {
n <- isolate(input$n)
conf <- isolate(input$conf)
tq <- qt(1 - (1 - conf) / 2, df = n - 1)
means <- numeric(k); hits <- logical(k)
for (i in seq_len(k)) {
idx <- sample(N, n)
s <- pop[idx]
m <- mean(s)
me <- tq * sqrt(1 - n / N) * sd(s) / sqrt(n)
means[i] <- m
hits[i] <- (m - me) <= ybarU && ybarU <= (m + me)
}
rv$idx <- idx
rv$mean <- m
rv$ci <- c(m - me, m + me)
rv$covered <- hits[k]
rv$means <- c(rv$means, means)
rv$hits <- c(rv$hits, hits)
}
draw_one <- function() draw_many(1)
# one sample per click; Play repeats until Pause or until the cap is reached
observeEvent(input$draw, { running(FALSE); draw_one() })
observeEvent(input$play, running(TRUE))
observeEvent(input$pause, running(FALSE))
# Skip: draw every remaining sample up to the cap in one go
observeEvent(input$skip, {
running(FALSE)
k <- nrep() - length(rv$means)
if (k > 0) draw_many(k)
})
observe({
if (!isTRUE(running())) return()
if (length(isolate(rv$means)) >= isolate(nrep())) {
running(FALSE)
return()
}
invalidateLater(1000 / isolate(input$speed), session)
isolate(draw_one())
})
# Reset clears the accumulated means and halts any run in progress
reset_history <- function() {
running(FALSE)
rv$idx <- NULL; rv$mean <- NULL; rv$ci <- NULL; rv$covered <- NULL
rv$means <- numeric(0); rv$hits <- logical(0)
}
observeEvent(input$reset, reset_history())
observeEvent(input$n, reset_history())
observeEvent(input$conf, reset_history())
fmt <- function(x) format(round(x), big.mark = ",")
output$mainplot <- renderPlot({
layout(matrix(1:2, nrow = 1), widths = c(1, 1))
## ---- left: population, with the current SRS sample (acres92 on the y axis) ----
par(mar = MAR, yaxs = "i")
plot(dens_y, dens$x, type = "l", col = "grey40", lwd = 2,
xlim = c(-0.62, 1.15), ylim = c(0, xmax_top),
xlab = "", ylab = "acres92 (thousands)", xaxt = "n", yaxt = "n", bty = "n",
main = "Population distribution with the current SRS sample")
# raw values, jittered; tail units sit in the compressed lane
points(rug_y, px, pch = 16, cex = 0.5, col = adjustcolor("grey60", 0.5))
if (!is.null(rv$idx))
points(rug_y[rv$idx], px[rv$idx], pch = 16, cex = 0.6,
col = adjustcolor("blue", 0.85))
# axis up to the cut, then the break marker and the lane
at <- pretty(c(0, cut), 6); at <- at[at <= cut]
axis(2, at = at, labels = format(at / 1000, big.mark = ",", trim = TRUE), las = 1)
segments(-0.62, cut, 1.15, cut, col = "grey75", lty = 3)
text(1.10, cut + lane_w / 2, sprintf("compressed: > %s (%d counties)", fmt(cut / 1000), n_tail),
adj = c(1, 0), col = "grey40", cex = 0.7)
# current sample mean and interval
if (!is.null(rv$mean)) {
segments(-0.55, rv$ci[1], -0.55, rv$ci[2], col = "blue", lwd = 3)
points(-0.55, rv$mean, pch = 19, cex = 1.3, col = "blue")
}
abline(h = ybarU, col = "black", lwd = 2, lty = 2)
legend("topright", bty = "o", bg = "white", box.col = NA, cex = 0.8, inset = c(0.01, 0.12),
legend = c("Population", "Current sample", "Sample mean and CI", "True mean"),
col = c("grey60", "blue", "blue", "black"),
pch = c(16, 16, 19, NA), lty = c(NA, NA, NA, 2))
## ---- right: sampling distribution of the sample mean (same y axis units) ----
par(mar = MAR, yaxs = "i")
if (length(rv$means) == 0) {
plot(1, type = "n", xlim = c(0, 1), ylim = xlim_b, bty = "n", xaxt = "n", yaxt = "n",
xlab = "", ylab = "Sample mean of acres92 (thousands)",
main = "Sampling distribution of the sample mean")
at_r <- pretty(xlim_b, 6); at_r <- at_r[at_r >= xlim_b[1] & at_r <= xlim_b[2]]
axis(2, at = at_r, labels = format(at_r / 1000, big.mark = ",", trim = TRUE), las = 1)
abline(h = ybarU, col = "black", lwd = 2, lty = 2)
text(0.5, ybarU, "Click 'Draw new sample' or Play to begin", col = "grey50")
return(invisible())
}
# fixed bin width; breaks extended to cover means outside the fixed window
bw <- 2 * half / 40
lo <- min(xlim_b[1], min(rv$means)) - bw
hi <- max(xlim_b[2], max(rv$means)) + bw
h <- hist(rv$means, breaks = seq(lo, hi, by = bw), plot = FALSE)
# theoretical sampling distribution of the mean at the current sample size
n <- input$n
se_theory <- sqrt(1 - n / N) * S / sqrt(n)
curve_pk <- dnorm(ybarU, mean = ybarU, sd = se_theory)
# empirical density of the accumulated means, once there are enough of them
dens_emp <- if (length(rv$means) >= 10) density(rv$means) else NULL
vis <- h$density[h$mids >= xlim_b[1] & h$mids <= xlim_b[2]]
top <- max(c(vis, curve_pk, if (!is.null(dens_emp)) dens_emp$y else NULL, 1e-6))
# histogram bars drawn horizontally: density along x, sample mean along y
plot(NA, xlim = c(0, top * 1.30), ylim = xlim_b, yaxt = "n", bty = "n",
xlab = "Density", ylab = "Sample mean of acres92 (thousands)",
main = "Sampling distribution of the sample mean")
at_r <- pretty(xlim_b, 6); at_r <- at_r[at_r >= xlim_b[1] & at_r <= xlim_b[2]]
axis(2, at = at_r, labels = format(at_r / 1000, big.mark = ",", trim = TRUE), las = 1)
rect(0, h$breaks[-length(h$breaks)], h$density, h$breaks[-1],
col = adjustcolor("steelblue", 0.55), border = "white")
if (!is.null(dens_emp)) lines(dens_emp$y, dens_emp$x, col = "blue", lwd = 2)
yg <- seq(xlim_b[1], xlim_b[2], length.out = 300)
lines(dnorm(yg, mean = ybarU, sd = se_theory), yg, col = "red", lwd = 2)
abline(h = ybarU, col = "black", lwd = 2, lty = 2)
# transient interval, for the current sample only (t-based, from its own sd)
ci_col <- if (rv$covered) "forestgreen" else "red"
x <- top * 0.90
segments(x, rv$ci[1], x, rv$ci[2], col = ci_col, lwd = 3)
segments(x - top * 0.045, rv$ci, x + top * 0.045, rv$ci, col = ci_col, lwd = 3)
points(x, rv$mean, pch = 19, cex = 1.3, col = ci_col)
# theoretical margin of deviation from the normal approximation, centred on
# mu and fixed for the current n and confidence level (not the sample CI);
# drawn close to the y-axis, at the base of the histogram bars
z_crit <- qnorm(1 - (1 - input$conf) / 2)
me_normal <- z_crit * se_theory
x2 <- top * 0.06
segments(x2, ybarU - me_normal, x2, ybarU + me_normal, col = "red", lwd = 2, lty = 2)
segments(x2 - top * 0.03, ybarU + c(-1, 1) * me_normal,
x2 + top * 0.03, ybarU + c(-1, 1) * me_normal, col = "red", lwd = 2)
outside <- sum(rv$means < xlim_b[1] | rv$means > xlim_b[2])
leg <- c(
sprintf("True mean = %sK", fmt(ybarU / 1000)),
sprintf("Theoretical SE (n = %d) = %sK", n, fmt(se_theory / 1000)),
sprintf("Normal-approx %.0f%% margin = ±%sK", 100 * input$conf, fmt(me_normal / 1000)),
sprintf("Samples drawn = %d of %d", length(rv$means), nrep()),
sprintf("Empirical coverage = %.1f%% (nominal %.0f%%)",
100 * mean(rv$hits), 100 * input$conf),
sprintf("Misses = %d of %d", sum(!rv$hits), length(rv$hits))
)
if (outside > 0)
leg <- c(leg, sprintf("%d mean(s) beyond the axis", outside))
legend("topright", bty = "n", cex = 0.8,
legend = c("Theoretical Normal", "Empirical density", "Normal-approx margin", leg),
lty = c(1, 1, 2, rep(NA, length(leg))),
lwd = c(2, 2, 2, rep(NA, length(leg))),
col = c("red", "blue", "red", rep(NA, length(leg))))
})
}
shinyApp(ui, server)
About this app
This app draws SRS samples one at a time, or continuously with play (▶) until pause (❚❚) or until the requested number of samples has been reached. The skip button (▶▶❘) draws all remaining samples at once, leaving the full histogram and the last sample on display.
The two panels sit side by side and share acres92 on the y-axis, so the two dashed lines at \(\mu\) line up horizontally across both. The left panel shows the population: the grey curve is a sideways kernel density estimate and the grey points to its left are the raw values, jittered horizontally. The axis is broken at twice the population mean; counties above the cut are compressed, in rank order, into the narrow lane at the top, so that a tail unit drawn into the sample is still visible. Points selected into the current sample turn blue, matching the sample mean and CI bar drawn beside them, and the sample means and density accumulated in the right panel.
The right panel accumulates the sample means as a blue density histogram (bars now run horizontally, density on the x-axis), with a blue curve tracing their empirical density once enough draws have accumulated, and a red curve for the theoretical sampling distribution \(N(\mu, \mathrm{SE}(\bar y))\) at the current sample size. Its window is fixed at five theoretical standard errors for \(n = 100\) on either side of \(\mu\), so the spread is comparable across sample sizes; for \(n < 100\) some means fall outside the window and are counted in the legend. The interval to the right belongs to the current sample, drawn green when it covers \(\mu\) and red when it does not, and replaced at every draw; the red dashed bar at the base of the histogram, close to the y-axis, is the theoretical margin of deviation \(\pm z_{\alpha/2}\,\mathrm{SE}(\bar y)\) from the normal approximation, centred on \(\mu\) and fixed for the current \(n\) and confidence level.
This app accompanies Shinylive App for Illustrating SRS Theory in Elements of Sampling Survey.