I build browser-based teaching simulations that run entirely in the browser (no installs), so students can experiment with the models directly. Explore them in my Econ Lab.

#| standalone: true
#| viewerHeight: 360

library(shiny)

# ---------------------------------------------------------------------------
# Helper: draw a single random sample from the chosen distribution
# ---------------------------------------------------------------------------
draw_sample <- function(n, dist) {
  switch(dist,
    "Uniform(0, 1)"      = runif(n),
    "Exponential(1)"     = rexp(n, rate = 1),
    "Right-skewed"       = rchisq(n, df = 3),
    "Bimodal"            = {
      k <- rbinom(n, 1, 0.5)
      k * rnorm(n, mean = -2, sd = 0.6) + (1 - k) * rnorm(n, mean = 2, sd = 0.6)
    },
    "Bernoulli(0.3)"     = rbinom(n, size = 1, prob = 0.3),
    runif(n)
  )
}

# Theoretical mean & sd of each population distribution
pop_params <- list(
  "Uniform(0, 1)"  = list(mu = 0.5, sigma = sqrt(1 / 12)),
  "Exponential(1)" = list(mu = 1,   sigma = 1),
  "Right-skewed"   = list(mu = 3,   sigma = sqrt(6)),
  "Bimodal"        = list(mu = 0,   sigma = sqrt(0.6^2 + 4)),
  "Bernoulli(0.3)" = list(mu = 0.3, sigma = sqrt(0.3 * 0.7))
)

# ---------------------------------------------------------------------------
# UI
# ---------------------------------------------------------------------------
ui <- fluidPage(
  tags$head(tags$style(HTML("
    .stats-box {
      background: #eaf2f8; border-radius: 6px; padding: 14px;
      margin-top: 12px; font-size: 14px; line-height: 1.8;
    }
    .stats-box b { color: #2c3e50; }
  "))),

  sidebarLayout(
    sidebarPanel(
      width = 3,

      selectInput("dist", "Population distribution:",
                  choices = names(pop_params)),

      sliderInput("n", "Sample size (n):",
                  min = 1, max = 200, value = 5, step = 1),

      sliderInput("reps", "Number of samples:",
                  min = 100, max = 3000, value = 1000, step = 100),

      actionButton("resample", "Draw new samples",
                   class = "btn-primary", width = "100%"),

      uiOutput("stats_box")
    ),

    mainPanel(
      width = 9,
      fluidRow(
        column(6, plotOutput("parent_plot", height = "230px")),
        column(6, plotOutput("sampling_plot", height = "230px"))
      )
    )
  )
)

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

  sim <- reactive({
    input$resample
    n    <- input$n
    reps <- input$reps
    dist <- input$dist

    means <- replicate(reps, mean(draw_sample(n, dist)))

    params <- pop_params[[dist]]
    theo_mu <- params$mu
    theo_se <- params$sigma / sqrt(n)

    list(means = means, dist = dist, n = n, reps = reps,
         theo_mu = theo_mu, theo_se = theo_se)
  })

  output$parent_plot <- renderPlot({
    dist <- input$dist
    big  <- draw_sample(10000, dist)

    par(mar = c(4.5, 4, 3, 1))
    hist(big, breaks = 60, probability = TRUE,
         col = "#d5e8d4", border = "#82b366",
         main = paste("True Population:", dist),
         xlab = "x", ylab = "Density")
  })

  output$sampling_plot <- renderPlot({
    s <- sim()

    par(mar = c(4.5, 4, 3, 1))
    hist(s$means, breaks = 40, probability = TRUE,
         col = "#dae8fc", border = "#6c8ebf",
         main = paste0("Sampling Distribution of the Mean (n = ", s$n, ")"),
         xlab = "Sample mean", ylab = "Density")

    x_seq <- seq(min(s$means), max(s$means), length.out = 300)
    lines(x_seq, dnorm(x_seq, mean = s$theo_mu, sd = s$theo_se),
          col = "#e74c3c", lwd = 2.5)

    abline(v = s$theo_mu, lty = 2, lwd = 2, col = "#2c3e50")

    legend("topright",
           legend = c("Normal approximation", "Theoretical mean"),
           col    = c("#e74c3c", "#2c3e50"),
           lwd    = c(2.5, 2),
           lty    = c(1, 2),
           bty    = "n", cex = 0.9)
  })

  output$stats_box <- renderUI({
    s <- sim()
    tags$div(class = "stats-box",
      HTML(paste0(
        "<b>Theoretical mean:</b> ",   round(s$theo_mu, 4), "<br>",
        "<b>Observed mean:</b> ",      round(mean(s$means), 4), "<br>",
        "<b>Theoretical SE:</b> ",     round(s$theo_se, 4), "<br>",
        "<b>Observed SD:</b> ",        round(sd(s$means), 4)
      ))
    )
  })
}

shinyApp(ui, server)

#| standalone: true
#| viewerHeight: 360

library(shiny)

ui <- fluidPage(
  titlePanel("Cobb-Douglas Consumer"),

  sidebarLayout(
    sidebarPanel(
      width = 4,
      sliderInput("a", "Cobb-Douglas exponent a:",
                  min = 0.5, max = 4, value = 1, step = 0.25),
      sliderInput("b", "Cobb-Douglas exponent b:",
                  min = 0.5, max = 4, value = 1, step = 0.25),
      sliderInput("m",  "Income (m):",            min =  20, max = 200, value = 100, step = 10),
      sliderInput("p1", "Price of good 1 (p₁):",  min =   1, max =  20, value =   5, step = 1),
      sliderInput("p2", "Price of good 2 (p₂):",  min =   1, max =  20, value =   5, step = 1),
      hr(),
      checkboxInput("show_other_ics",
                    "Show a ladder of indifference curves",
                    value = TRUE),
      htmlOutput("readout")
    ),
    mainPanel(
      width = 8,
      plotOutput("choice_plot", height = "300px")
    )
  )
)

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

  bundle <- reactive({
    a <- input$a; b <- input$b
    m <- input$m; p1 <- input$p1; p2 <- input$p2
    x1 <- a * m / ((a + b) * p1)
    x2 <- b * m / ((a + b) * p2)
    u_opt <- x1^a * x2^b
    list(a = a, b = b, m = m, p1 = p1, p2 = p2,
         x1 = x1, x2 = x2, u_opt = u_opt)
  })

  output$choice_plot <- renderPlot({
    s <- bundle()
    x1_int <- s$m / s$p1
    x2_int <- s$m / s$p2

    xmax <- max(40, x1_int) * 1.05
    ymax <- max(40, x2_int) * 1.05

    par(mar = c(4.2, 4.5, 1, 1))
    plot(NA, xlim = c(0, xmax), ylim = c(0, ymax),
         xlab = expression(x[1]), ylab = expression(x[2]), main = "")

    # Indifference curves: x_1^a x_2^b = U => x_2 = (U / x_1^a)^(1/b)
    ic_x <- seq(0.1, xmax * 1.5, length.out = 600)
    ic_y_opt <- (s$u_opt / ic_x^s$a)^(1 / s$b)

    if (isTRUE(input$show_other_ics)) {
      for (frac in c(0.5, 0.75, 1.25)) {
        u <- s$u_opt * frac
        y <- (u / ic_x^s$a)^(1 / s$b)
        lines(ic_x, y, col = adjustcolor("#7f8c8d", 0.5), lwd = 1.2, lty = 3)
      }
    }

    lines(ic_x, ic_y_opt, col = "#27ae60", lwd = 2.4)

    polygon(c(0, x1_int, 0), c(0, 0, x2_int),
            col = adjustcolor("#3498db", 0.12), border = NA)
    segments(0, x2_int, x1_int, 0, col = "#185FA5", lwd = 3)

    points(s$x1, s$x2, pch = 19, col = "#c0392b", cex = 1.8)
    text(s$x1, s$x2, sprintf("  (%.1f, %.1f)", s$x1, s$x2),
         pos = 4, col = "#c0392b", cex = 1.0)

    legend("topright",
           legend = c("Budget line", "Indifference curve through optimum",
                      if (isTRUE(input$show_other_ics)) "Other indifference curves" else NULL,
                      "Optimal bundle (x₁*, x₂*)"),
           col    = c("#185FA5", "#27ae60",
                      if (isTRUE(input$show_other_ics)) "#7f8c8d" else NULL,
                      "#c0392b"),
           lwd    = c(3, 2.4, if (isTRUE(input$show_other_ics)) 1.2 else NULL, NA),
           lty    = c(1, 1, if (isTRUE(input$show_other_ics)) 3 else NULL, NA),
           pch    = c(NA, NA, if (isTRUE(input$show_other_ics)) NA else NULL, 19),
           bty = "n", cex = 0.9)
  })

  output$readout <- renderUI({
    s <- bundle()
    share1 <- s$a / (s$a + s$b)
    share2 <- s$b / (s$a + s$b)
    HTML(sprintf(paste(
      "<div style='margin-top:10px;font-size:13px;line-height:1.7;'>",
      "<b>Ordinary demand:</b> x₁* = %.2f, x₂* = %.2f<br>",
      "<b>Spending on good 1:</b> p₁·x₁* = %.1f (%.0f%% of income)<br>",
      "<b>Spending on good 2:</b> p₂·x₂* = %.1f (%.0f%% of income)<br>",
      "<b>MRS at optimum:</b> %.2f &nbsp; <b>p₁/p₂:</b> %.2f<br>",
      "<b>Utility achieved:</b> %.2f",
      "</div>"
    ),
    s$x1, s$x2,
    s$p1 * s$x1, 100 * share1,
    s$p2 * s$x2, 100 * share2,
    (s$a / s$b) * (s$x2 / s$x1),
    s$p1 / s$p2, s$u_opt))
  })
}

shinyApp(ui, server)