Random Sampling — Generate random samples with...

SECTION 1

Historical Context & Motivation

The ability to generate random samples computationally is one of the most consequential developments in modern statistics and computer science. Before electronic computers, researchers who needed random numbers relied on physical devices—dice, shuffled cards, or tables of digits laboriously compiled by hand. The demand for large-scale random sampling grew dramatically during World War II, when scientists at Los Alamos needed to simulate neutron transport through fissile material, a problem far too complex for closed-form analysis.

That wartime necessity led directly to the Monte Carlo method, pioneered by Stanislaw Ulam and John von Neumann, which depends entirely on the rapid generation of pseudorandom numbers. As statistical computing matured, languages like S (and later R) embedded random number generation as a first-class primitive, giving analysts direct access to functions such as rnorm(), runif(), and sample(). Understanding these functions is prerequisite to simulation studies, bootstrapping, permutation tests, and stochastic optimization—tools that pervade both applied statistics and algorithm design.

1946

Birth of Monte Carlo Methods

Stanislaw Ulam and John von Neumann formalize the Monte Carlo method at Los Alamos, creating the first systematic need for large volumes of pseudorandom numbers in scientific computation.

1951

Lehmer's Linear Congruential Generator

Derrick Lehmer proposes the linear congruential generator (LCG), providing a fast, deterministic algorithm for pseudorandom uniform variates that becomes the backbone of early simulation software.

1976

S Language at Bell Labs

John Chambers and colleagues create the S language, embedding random variate generation (including rnorm and runif predecessors) as core statistical primitives.

1997

Mersenne Twister Published

Matsumoto and Nishimura publish the Mersenne Twister (MT19937), which R later adopts as its default PRNG due to its extremely long period of 2¹⁹⁹³⁷ − 1 and excellent equidistribution properties.

2000

R 1.0 Released

R 1.0.0 ships with a comprehensive suite of distribution functions following the d/p/q/r naming convention, making rnorm(), runif(), and sample() available to the open-source community at scale.

The central question this lesson addresses is both practical and conceptual: how does R translate a deterministic algorithm into draws that behave as if they came from a specified probability distribution, and how should you wield rnorm(), runif(), and sample() to produce correct, reproducible simulations?

SECTION 2

Core Principles & Definitions

Before invoking any sampling function in R, you need to grasp several foundational concepts that govern how pseudorandom numbers are produced and why they are trustworthy enough for rigorous statistical inference. R's random number infrastructure rests on a pseudorandom number generator (PRNG) that deterministically advances an internal state to produce a stream of values indistinguishable, by standard statistical tests, from true randomness. The following principles form the conceptual scaffold for everything that follows.

Pseudorandom Number Generator (PRNG)

A deterministic algorithm that, given a seed, produces a sequence of numbers that passes statistical randomness tests. R's default PRNG is the Mersenne Twister (MT19937) with a period of 2¹⁹⁹³⁷ − 1.

Seed & Reproducibility

Calling set.seed(n) initializes the PRNG state so that subsequent draws are identical across runs and machines, enabling reproducible research and debugging.

The d/p/q/r Convention

For every distribution in R, four functions exist: d (density), p (cumulative probability), q (quantile), and r (random variate). Thus rnorm generates Normal variates, runif generates Uniform variates, and so on.

Continuous vs. Discrete Sampling

rnorm() and runif() draw from continuous distributions (Normal and Uniform, respectively), while sample() draws discrete elements from a vector, optionally with replacement and with user-defined probability weights.

Inverse Transform & Box-Muller

Continuous variates are derived from uniform variates via transformation methods. The inverse transform method applies the quantile function; the Box-Muller transform converts two Uniform(0,1) draws into two independent Normal(0,1) draws.

✦ KEY TAKEAWAY

Think of R's PRNG like a very long, pre-shuffled deck of cards: calling set.seed() selects which deck you pick up, and each call to rnorm() or runif() draws the next card. The sequence looks random, but it is completely determined by the seed—exactly the property you need for reproducible simulation experiments.

SECTION 3

Visual Explanation — How the Three Functions Map to Distributions

The diagram below illustrates the relationship between R's PRNG engine and the three primary sampling functions. The Mersenne Twister produces a stream of uniform variates on (0, 1). From there, runif() rescales to an arbitrary interval [a, b], rnorm() applies a transformation (typically Box-Muller or Kinderman-Ramage) to produce Normal variates, and sample() uses the uniform stream to index discrete elements from a user-supplied vector.

The pipeline starts with set.seed() initializing the Mersenne Twister, which produces a stream of U(0,1) variates. Each function branches from this stream: runif() rescales to [a, b], rnorm() applies a Normal transformation, and sample() selects discrete elements.

Notice how the entire pipeline is deterministic once the seed is fixed. If you call set.seed(42) followed by rnorm(5) on any machine running the same version of R with the default PRNG, you will obtain the same five values. This reproducibility is critical for peer review, unit testing of stochastic algorithms, and debugging simulation code. The diagram also clarifies that sample() is fundamentally different from the other two: it operates on a finite, user-supplied vector rather than a continuous probability distribution.

SECTION 4

Mathematical Framework — Transformations Under the Hood

All of R's random variate generators ultimately rely on draws from a Uniform(0, 1) distribution produced by the PRNG. Understanding the mathematical transformations that convert these raw uniform values into other distributions gives you deeper control over your simulations and helps you reason about edge cases.

runif — Linear Rescaling

UNIFORM RESCALING

X = a + U × (b − a), where U ~ Uniform(0, 1)

Here a is the lower bound (default 0), b is the upper bound (default 1), and X ~ Uniform(a, b). Calling runif(n, min = a, max = b) applies this affine transformation to n independent U(0,1) draws.

rnorm — Box-Muller Transform

BOX-MULLER TRANSFORM

Z₁ = √(−2 ln U₁) × cos(2πU₂), Z₂ = √(−2 ln U₁) × sin(2πU₂)

Given two independent draws U₁, U₂ ~ Uniform(0, 1), the transform produces two independent standard Normal variates Z₁, Z₂ ~ N(0, 1). R then rescales: X = μ + σZ to produce N(μ, σ²) draws via rnorm(n, mean = μ, sd = σ).

sample — Discrete Selection Mechanism

WEIGHTED DISCRETE SAMPLING

P(select element xᵢ) = wᵢ / Σⱼ wⱼ

When a probability vector prob is supplied, R normalizes the weights so they sum to 1 and uses them to bias selection. Without prob, all elements are equally likely. The replace argument controls whether sampling is with or without replacement.

💡 Inverse Transform Method (General)

For any continuous distribution with CDF F, setting X = F⁻¹(U) where U ~ Uniform(0,1) produces a variate X with CDF F. This is the theoretical basis for R's entire r* family. While R may use more efficient algorithms (e.g., Ahrens-Dieter for Gamma), the inverse transform principle underpins them all.

SECTION 5

Detailed Breakdown — Function Signatures and Parameters

Each of the three primary sampling functions in R has a distinct signature with different defaults and semantics. The table below provides a comprehensive reference, followed by a diagram that visualizes the output characteristics of each function side by side.

Primary random sampling functions in R

Function	Signature	Distribution / Domain	Key Parameters
`runif()`	`runif(n, min=0, max=1)`	Continuous Uniform on [min, max]	`n`: count; `min`, `max`: interval bounds
`rnorm()`	`rnorm(n, mean=0, sd=1)`	Normal (Gaussian) on (−∞, ∞)	`n`: count; `mean`: center μ; `sd`: standard deviation σ
`sample()`	`sample(x, size, replace=FALSE, prob=NULL)`	Discrete elements from vector `x`	`x`: source vector; `size`: draw count; `replace`: with replacement?; `prob`: weight vector

Top row: histograms of 10,000 draws from each function illustrate their characteristic shapes—flat for runif(), bell-shaped for rnorm(), and discrete bars for sample(). Bottom row: common use-case guidance for each function.

A few subtleties worth noting for computer science students: when sample() is called with a single integer argument like sample(10), R interprets this as sample(1:10, 10)—a random permutation of 1 through 10. This is a common source of bugs: if x is a length-1 vector containing a large integer, you will not get a single-element sample but rather a permutation of 1:x. Guard against this edge case by always specifying size explicitly and, when appropriate, wrapping your vector in sample.int() or indexing defensively.

SECTION 6

Worked Example — Monte Carlo Estimation of π

A classic demonstration that integrates all three functions is the Monte Carlo estimation of π. We generate random points inside a unit square, count how many fall inside the inscribed quarter-circle, and use the ratio to approximate π. This example combines runif() for coordinate generation, rnorm() for adding Gaussian noise in a follow-up, and sample() for subsampling results.

Estimating π with runif(), then analyzing with sample()

Step 1 — Set the seed for reproducibilityWe begin by fixing the random seed so the results are identical on every run. This is essential for debugging and peer review of simulation studies. set.seed(2024)

Step 2 — Generate random (x, y) coordinates with runif()We draw n = 100,000 uniform random points in the unit square [0, 1] × [0, 1]. Each coordinate pair is an independent draw from Uniform(0, 1).

n <- 100000
x <- runif(n, min = 0, max = 1)
y <- runif(n, min = 0, max = 1)

Step 3 — Test which points fall inside the quarter-circleA point (x, y) lies inside the quarter-circle of radius 1 if x² + y² ≤ 1. We create a logical vector recording this test for all n points. inside <- (x^2 + y^2) <= 1

Step 4 — Estimate π from the ratioThe area of the quarter-circle is π/4, and the area of the unit square is 1. Therefore, the fraction of points inside approximates π/4, and π ≈ 4 × (hits / n).

pi_estimate <- 4 * sum(inside) / n
cat("Estimated pi:", pi_estimate)

Output: Estimated pi: 3.14112 (close to π ≈ 3.14159)

Step 5 — Subsample and add noise with sample() and rnorm()To simulate measurement uncertainty, we select 500 of the inside-points at random using sample(), then perturb their x-coordinates with Gaussian noise via rnorm().

idx_inside <- which(inside)
chosen <- sample(idx_inside, size = 500, replace = FALSE)
x_noisy <- x[chosen] + rnorm(500, mean = 0, sd = 0.01)

The vector x_noisy now contains 500 perturbed x-coordinates drawn from the interior of the quarter-circle—a realistic model of noisy sensor data.

SECTION 7

Strengths, Limitations & Pitfalls

While R's random sampling functions are powerful and convenient, each has specific strengths and limitations that become critical in production-grade simulation pipelines. The table below contrasts the three functions along dimensions relevant to real-world statistical computing and algorithm design.

Comparative analysis of R's three primary random sampling functions

Dimension	runif()	rnorm()	sample()
Domain	Continuous [min, max]	Continuous (−∞, ∞)	Discrete (finite vector)
Speed (n = 10⁷)	Very fast (~0.3 s)	Fast (~0.5 s)	Moderate (~1.0 s with replace)
Common pitfall	Forgetting that endpoints may or may not be included (implementation detail)	Confusing sd (σ) with variance (σ²); passing variance instead of sd	Single-integer trap: sample(10) ≠ sample(c(10), 1)
Extensibility	Basis for inverse transform; use qnorm(runif(n)) for custom transforms	Generates multivariate normal via MASS::mvrnorm(); key to Bayesian MCMC	Powers bootstrap, cross-validation splits, and permutation tests
Limitation	Cannot model natural phenomena that cluster around a mean	Symmetric tails; not suitable for skewed or heavy-tailed data	Without replacement limited to length(x) draws; prob weights must be non-negative

⚡ CHOOSING THE RIGHT TOOL

Think of these three functions as analogous to three data structures: runif() is like a raw array (primitive, general-purpose), rnorm() is like a specialized hash map (optimized for a specific access pattern—the bell curve), and sample() is like a set operation (selecting elements from a collection). Use the most semantically appropriate one for your problem; never use runif() with manual rounding when sample() exists for discrete draws.

SECTION 8

Connection to Advanced Topics — Bootstrapping & MCMC

The three sampling primitives you have learned are not merely pedagogical tools—they are the operational foundation for some of the most powerful techniques in modern statistical computing. Two areas where these functions play a central role are bootstrapping and Markov Chain Monte Carlo (MCMC) methods.

How basic R sampling functions connect to advanced statistical computing

Concept	This Lesson's Functions	Advanced Extension
Bootstrap confidence intervals	`sample(data, n, replace = TRUE)` draws B resamples of the dataset	The `boot` package automates this; BCa intervals, parametric bootstrap, and wild bootstrap extend the basic idea
Permutation tests	`sample(labels)` shuffles group assignments under H₀	Exact permutation enumeration; approximate p-values via Monte Carlo; `coin` package
Metropolis-Hastings (MCMC)	`rnorm()` proposes candidate moves; `runif(1)` decides acceptance	Hamiltonian Monte Carlo (HMC), No-U-Turn Sampler (NUTS) in Stan/RStan; Gibbs sampling
Cross-validation	`sample(1:n, k)` partitions data into k folds	Stratified k-fold, repeated CV, nested CV for hyperparameter tuning via `caret` or `tidymodels`
Stochastic gradient descent	`sample(1:n, batch_size)` selects mini-batches	Adam, RMSProp, and learning rate schedules; integration with `torch` for R

As you advance into Bayesian inference, machine learning pipelines, or large-scale simulation, you will find that nearly every stochastic algorithm reduces to repeated calls to the primitives covered in this lesson. Mastering runif(), rnorm(), and sample() now gives you a composable toolkit for building arbitrarily complex stochastic systems later.

SECTION 9

Practice Problems

PROBLEM 1 — CONCEPTUAL

Explain why calling set.seed(123) before rnorm(5) always produces the same five values. What would happen if you inserted a call to runif(1) between set.seed(123) and rnorm(5)? Would the five Normal values change? Why?

PROBLEM 2 — BASIC CALCULATION

Write R code to generate 1,000 draws from a Uniform distribution on [5, 15] and verify empirically that the mean is approximately 10 and the standard deviation is approximately 2.887 (which equals (15 − 5) / √12). Use set.seed(1) for reproducibility.

PROBLEM 3 — INTERMEDIATE

Suppose you have a dataset of 200 observations stored in a vector data. Write R code to perform a nonparametric bootstrap with B = 5,000 replicates to estimate the 95% confidence interval for the median. Use sample() with replacement inside a replicate() call.

PROBLEM 4 — APPLIED

You are simulating a sensor network where each of 50 sensors reports a temperature reading corrupted by Gaussian noise. The true temperature is 22.5°C, and each sensor has independent noise with mean 0 and standard deviation 0.8°C. Write R code to: (a) simulate one round of readings for all 50 sensors, (b) compute the sample mean and its standard error, and (c) use sample() to randomly select 10 sensors whose readings you will transmit (simulating bandwidth constraints).

PROBLEM 5 — CRITICAL THINKING

A colleague claims: 'I can replace rnorm(n, 5, 2) with 5 + 2 * runif(n, -3, 3) because both produce values centered at 5 with a similar spread.' Critically evaluate this claim. What are the distributional differences, and in what analysis scenario would substituting one for the other lead to incorrect inferences?

SUMMARY

Lesson Summary

R's random sampling infrastructure is built on the Mersenne Twister PRNG, which produces a deterministic stream of Uniform(0, 1) variates initialized by set.seed(). From this stream, runif(n, min, max) rescales to any continuous uniform interval, rnorm(n, mean, sd) applies the Box-Muller transform to generate Gaussian variates, and sample(x, size, replace, prob) selects discrete elements with optional weights and replacement.

These three functions form the composable primitives underlying Monte Carlo simulation, bootstrap resampling, MCMC sampling, and cross-validation splits. Always set a seed for reproducibility, watch for the single-integer trap in sample(), and remember that rnorm() takes standard deviation (σ), not variance (σ²).