R Bootcamp: Day 3 PM


Probability (sort of) and simulation in R


Eric Manning

September 1, 2026



princeton-ddss.github.io/r-bootcamp

This session


1) Drawing random numbers

  • sample(), set.seed(), and a couple of named distributions
  • The d/p/q/r family


2) Simulation

  • Wax on, wax off


We are not teaching statistics today. We are showing you the machinery, so the concepts seem vaguely familiar later.

Some sources



Additional: Geyer, Distributions in R — a one-page lookup table for every distribution R ships with. Bookmark it.

Drawing random numbers

sample(): draw from a vector


Give it something to draw from and how many draws you want:

sample(1:10, size = 5)
[1] 10  4  8  2  9


With no size, you get a permutation — the whole vector, reshuffled:

sample(1:10)
 [1] 10  5  2  9  3  8  6  4  7  1

With or without replacement


By default, each element can be drawn once. size can’t exceed the vector:

sample(1:5, size = 10)
Error in `sample.int()`:
! cannot take a sample larger than the population when 'replace = FALSE'


replace = TRUE puts each draw back in the bag:

sample(1:5, size = 10, replace = TRUE)
 [1] 2 4 1 5 1 3 2 3 1 1

A coin, ten times


The vector doesn’t have to be numbers:

sample(c("heads", "tails"), size = 10, replace = TRUE)
 [1] "tails" "heads" "heads" "tails" "tails" "heads" "heads" "heads" "heads"
[10] "heads"

An unfair coin: prob


prob weights the draws. It lines up with the vector you’re sampling from:

sample(c("heads", "tails"), size = 10, replace = TRUE, prob = c(0.9, 0.1))
 [1] "heads" "heads" "heads" "heads" "heads" "heads" "heads" "heads" "heads"
[10] "heads"


prob gets normalized, so c(9, 1) and c(0.9, 0.1) behave identically.

Your results won’t match mine


Run the same line twice:

sample(1:100, 4)
sample(1:100, 4)
[1] 32 74  6 51
[1] 44 29 63 14


Results are not reproducible between runs.

set.seed()


R’s randomness is not random — it’s a deterministic sequence that looks random. set.seed() says where to start:

set.seed(2026)
sample(1:100, 4)

set.seed(2026)
sample(1:100, 4)
[1] 93 97 38 45
[1] 93 97 38 45

Using set.seed() in practice


  • Set it once, near the top of the script, with any integer you like
  • Then never touch it again in that script

Two distributions worth naming


Bernoulli

One trial, two outcomes, probability p of a 1/T/‘yes’/etc. A (biased?) coin flip.


Uniform

Every value in a range is equally likely.

Bernoulli draws: rbinom()


size = 1 means one coin per draw; prob is the chance of a 1:

set.seed(2026)
rbinom(n = 20, size = 1, prob = 0.5)
 [1] 1 1 0 0 1 0 0 1 0 1 0 1 0 1 0 0 1 0 0 0


Bump size and each draw becomes “how many heads in size flips?”:

rbinom(n = 10, size = 100, prob = 0.5)
 [1] 54 54 56 57 45 53 50 54 50 54

Uniform draws: runif()


runif(n) gives n values between 0 and 1, none more likely than any other:

set.seed(2026)
runif(5)
[1] 0.6986735 0.5565305 0.1401400 0.2857233 0.5553690


min and max move the range:

runif(5, min = 10, max = 20)
[1] 10.25131 14.66231 18.61011 12.52501 15.80806

What “uniform” looks like


set.seed(2026)
hist(runif(10000), breaks = 40, main = "10,000 draws from runif()")

And now one we’re not going to explain


set.seed(2026)
hist(rnorm(10000), breaks = 40, main = "10,000 draws from rnorm()")

This is the normal distribution. Note the shape and move on.

rnorm() arguments


mean centers it, sd sets its width. The defaults are 0 and 1:

set.seed(2026)
rnorm(5)
rnorm(5, mean = 100, sd = 15)
[1]  0.52058907 -1.07969076  0.13923812 -0.08474878 -0.66663962
[1]  62.25866  88.97280  84.69817 101.70332  92.89314


Every r*() function follows this shape: how many draws, then the knobs for that distribution.

The d / p / q / r family

One distribution, four questions


Every distribution in R comes as four functions sharing one name stem:

Prefix Question it answers Normal
r Give me random draws rnorm()
d How dense/tall is the curve here? dnorm()
p How much of the distribution’s mass lies below here? pnorm()
q Which value has this much below it? qnorm()

Visualized

Jason Mercer, “A brief visualization of R’s distribution functions” (2020).

Reading that figure


Top row works on the density curve, bottom row on the cumulative curve:

  • dnorm() — you hand it a value, it hands back the height of the curve
  • rnorm() — draws values, more often where the curve is tall
  • pnorm() — you hand it a value, it hands back the probability below it
  • qnorm() — the inverse: hand it a probability, get the value


p and q undo each other.

p and q


Half the curve lies below its center:

pnorm(0)
[1] 0.5


Go the other way — which value has 2.5% of the curve below it?

qnorm(0.025)
[1] -1.959964


And back again:

pnorm(qnorm(0.025))
[1] 0.025

d is not a probability


dnorm(0)
[1] 0.3989423


0.399 is a height For a continuous distribution the probability of any exact value is zero.


Probability lives in areas under the curve (pnorm()). dnorm() is for drawing the curve.

Drawing a density with d


stat_function() takes the function itself — no data frame needed:

ggplot() +
    stat_function(fun = dnorm, xlim = c(-4, 4), linewidth = 1) +
    labs(x = NULL, y = "Density") +
    theme_minimal(base_size = 20)

p is the shaded area


ggplot() +
    stat_function(fun = dnorm, xlim = c(-4, 4), linewidth = 1) +
    stat_function(
        fun = dnorm, xlim = c(-4, qnorm(0.025)),
        geom = "area", fill = "steelblue"
    ) +
    labs(x = NULL, y = "Density") +
    theme_minimal(base_size = 20)

The blue strip runs from the far left up to qnorm(0.025) = -1.96, and its area is pnorm(-1.96) = 0.025.

“2.5% of the distribution sits below -1.96.”

Every other distribution


?Distributions lists every one R ships with. Geyer’s lookup table is friendlier.

Simulation

Why


Some probability questions have clean formulas. Many do not.

Recipe


1. Do the random thing once

2. Wrap it in a function that returns one number or TRUE/FALSE

3. Run that function many times

4. Summarize the results

Step 1: do it once


What’s the probability two dice sum to 7?

set.seed(2026)
dice <- sample(1:6, size = 2, replace = TRUE)
dice
sum(dice) == 7
[1] 5 1
[1] FALSE

Step 2: wrap it in a function


roll_is_7 <- function() {
    dice <- sample(1:6, size = 2, replace = TRUE)
    sum(dice) == 7
}

roll_is_7()
[1] TRUE

Step 3: replicate()


replicate(n, expr) runs an expression n times and collects the results:

set.seed(2026)
rolls <- replicate(10000, roll_is_7())

length(rolls)
head(rolls)
[1] 10000
[1] FALSE  TRUE FALSE FALSE FALSE FALSE


This is just a for loop.

Step 4: summarize


mean(rolls)
[1] 0.1654


The true answer is 6/36 = 0.1667.

How many replications?


Run the whole simulation at different sizes and watch it settle:

set.seed(2026)
sapply(
    c(10, 100, 1000, 10000, 100000),
    function(n) mean(replicate(n, roll_is_7()))
)
[1] 0.20000 0.15000 0.16300 0.16760 0.16661


More replications == more precision in your estimate.

Resampling your own data

Randomness from data, not a distribution


set.seed(2026)
sample(counties$income, size = 5)
[1] 62496 57056 50635 57504 69826


Why?

1. You have the population. Draw samples (without replacement) to explore sampling variability. You choose (n), and seeing what changes as (n) changes is the point.


2. You have one sample of a population. Resample its observations with replacement, treating the resample as a stand-in for the population. This is the bootstrap. Variation among the bootstrap iterations helps us quantify uncertainty in the estimate derived from the observed sample.

The bootstrap


You have 50 observations. You want to know how much your estimate would have varied if you’d drawn a different 50.

You can’t — there’s only one dataset, but we can treat the sample as if it were the population and redraw from it at the same size:

set.seed(2026)
my_sample <- sample(counties$income, size = 50)

mean(sample(my_sample, size = 50, replace = TRUE))
[1] 64069.1

Do that 10,000 times.

If you don’t bootstrap with the same size, you have to apply a correction factor.

Useful tools


  • boot::boot() — the classic implementation, several kinds of interval; will apply correction factors
  • rsample::bootstraps() — tidyverse-flavored


Both do exactly what we just did – and other stuff.

Big data paradox


Meng (2018) examines a 2016 election survey with 2.3 million respondents — about 1% of the American electorate.


Its tiny non-randomness meant it carried the same error as a genuine random sample of about 400 people.

Simulation exercises


Render as you go. No AI tools.
princeton-ddss.github.io/r-bootcamp/notes/day-3-sim-ex.pdf


Function What it does
sample(x, size, replace, prob) Draw from a vector
set.seed(n) Make the draws reproducible
runif(n, min, max) Uniform draws
rbinom(n, size, prob) Coin flips / counts of successes
rnorm(n, mean, sd) Normal draws
dnorm() / pnorm() / qnorm() Height / area below / inverse
replicate(n, expr) Run an expression n times, keep results
mean(logical_vector) Proportion TRUE
quantile(x, probs) Cut points of a set of numbers