Introduction: The Actuarial "Flight Simulator"
Welcome to one of the most practical chapters in the CS1 curriculum! In Paper B, you aren't just calculating probabilities; you are building models. Think of simulating random variables as a flight simulator for actuaries. Before an insurance company launches a new product, they "test fly" it by simulating thousands of possible claims to see if the company stays solvent. In this chapter, we will learn how to use R to generate these "test flights" using built-in functions and the Inverse Transform Method.
1. The "r" Family of Functions
In R, most probability distributions have a dedicated function for generating random samples. These functions always start with the letter r (for "random").
The general structure is usually: r[distribution_name](n, [parameters]), where \(n\) is the number of observations you want to generate.
Common Discrete Distributions
- Binomial: rbinom(n, size, prob) — generates \(n\) samples from \(X \sim Bin(size, prob)\).
- Poisson: rpois(n, lambda) — generates \(n\) samples from \(X \sim Pois(\lambda)\).
- Geometric: rgeom(n, prob) — generates the number of failures before the first success.
- Negative Binomial: rnbinom(n, size, prob) — generates the number of failures before a certain number of successes occur.
Common Continuous Distributions
- Normal: rnorm(n, mean, sd) — Note that R uses standard deviation \(\sigma\), not variance \(\sigma^2\).
- Exponential: rexp(n, rate) — Uses the rate parameter \(\lambda\).
- Gamma: rgamma(n, shape, rate) — Can also use scale (\(1/rate\)).
- Uniform: runif(n, min, max) — Generates values between a lower and upper bound.
- Lognormal: rlnorm(n, meanlog, sdlog).
- Beta: rbeta(n, shape1, shape2).
Quick Review: If you need to simulate 100 claims from a Poisson distribution with a mean of 5, you would use rpois(100, 5).
2. Reproducibility: The set.seed() Function
Computers don't actually generate truly "random" numbers; they use algorithms to create "pseudo-random" numbers. If you run rnorm(5, 0, 1) twice, you will get different results. This is a problem for examiners who need to mark your work!
The Solution: Use set.seed(). By typing set.seed(123) (or any integer) before your simulation code, you ensure that R produces the exact same "random" numbers every time the code is run.
Important Exam Tip: Always check if the exam paper specifies a seed. If it says "Use a seed of 42," your first line of code must be set.seed(42).
3. The Inverse Transform Method (ITM)
Sometimes, you might be asked to generate a random variable manually without using the specific "r" function. The Inverse Transform Method is the standard way to do this for both discrete and continuous variables.
The Logic Behind ITM
Every Cumulative Distribution Function (CDF), \(F(x)\), outputs a value between 0 and 1. If we take a random number \(U\) from a Uniform(0,1) distribution, we can "work backward" to find the value of \(x\) that corresponds to that probability.
Step-by-Step Process:
- Generate a random number \(u\) from \(U \sim Uniform(0,1)\) using runif(1).
- Set \(F(x) = u\).
- Solve for \(x\) by calculating the inverse: \(x = F^{-1}(u)\).
Example for Exponential Distribution:
The CDF is \(F(x) = 1 - exp(-\lambda x)\).
1. Set \(u = 1 - exp(-\lambda x)\)
2. \(1 - u = exp(-\lambda x)\)
3. \(\ln(1 - u) = -\lambda x\)
4. \(x = -\frac{1}{\lambda} \ln(1 - u)\)
In R, if you have a vector of uniform random variables u, you would code this as: x <- - (1/lambda) * log(1 - u).
Did you know? Since \(U\) and \(1-U\) both follow the \(Uniform(0,1)\) distribution, actuaries often simplify the formula to \(x = -\frac{1}{\lambda} \ln(u)\).
4. Sampling from Existing Data
The syllabus mentions generating samples not just from theoretical distributions, but from existing data sets. We use the sample() function for this.
Key Arguments for sample():
- x: The data or vector to sample from.
- size: How many items to pick.
- replace: Should we put the item back after picking it? (TRUE or FALSE).
Real-World Application:
The Bootstrap: To estimate the properties of an estimator (like its bias or variance), we can use the "Bootstrap method." This involves sampling from our data with replacement (replace = TRUE). This mimics taking new samples from the population.
Permutation Tests: To perform non-parametric hypothesis tests, we often sample without replacement (replace = FALSE) to "shuffle" the data and see if patterns happen by chance.
5. Comparing Simulations to Theory
Once you have simulated data, the syllabus requires you to compare it to known distributions, especially the Normal Distribution as part of the Central Limit Theorem (CLT).
If you simulate 1,000 sample means from any distribution, the CLT tells us those means should look like a Normal curve. In R, you can check this by:
- Plotting a histogram of your simulated data: hist(sim_data).
- Overlaying a density curve: lines(density(sim_data)).
- Using a Q-Q plot to see if the points fall on a straight line (covered further in the "Exploratory plots" chapter).
Key Takeaway: Simulation allows us to verify the CLT. Even if the underlying data is heavily skewed (like an Exponential distribution), the average of those simulations will start to look "Normal" as your sample size increases.
Summary & Quick Review
Common Pitfalls to Avoid:
- Standard Deviation vs Variance: Remember rnorm() needs \(\sigma\), but many exam questions provide \(\sigma^2\). Use sqrt() if necessary!
- Forgetting set.seed(): If you don't set a seed, your numerical answers will differ from the mark scheme.
- ITM Logic: Remember that ITM always starts with runif(). You are transforming a "probability" back into a "value."
Note: For help on calculating specific probabilities or summary statistics of these simulated samples, please refer to the chapter on "Probabilities, quantiles and summary statistics in R".