We’re in the 21st century!
Simulations can often be
easier than hand calculations
made more realistic than hand calculations
R provides unique access to great (statistical) simulation tools (compared to other languages)
PDF: \(\displaystyle f(x \mid \mu, \sigma^2)=\frac{1}{\sigma \sqrt{2 \pi}} e^{-\frac{1}{2}\left(\frac{x-\mu}{\sigma}\right)^{2}}\); \(x \in (- \infty, \infty)\)
We write \(X \sim N(\mu, \sigma^2)\)
Standard normal distribution: \(N(0, 1)\)
library(tidyverse)
theme_set(theme_light())
tibble(x = c(-5, 5)) |>
ggplot(aes(x)) +
stat_function(fun = dnorm, color = "blue",
args = list(mean = 0, sd = 1)) +
stat_function(fun = dnorm, color = "red",
args = list(mean = -2, sd = sqrt(0.5)))
PMF: \(\displaystyle f(x \mid n, p)= \binom{n}{x} p^{x}(1-p)^{n-x}\); \(x = 0, 1, \dots, n\)
Model for the probability of \(x\) successes in \(n\) independent trials, each with success probability \(p\)
We write \(X \sim \text{Binomial}(n,p)\)
tibble(x = 0:20) |>
mutate(binom1 = dbinom(x, size = 20, prob = 0.5),
binom2 = dbinom(x, size = 20, prob = 0.1)) |>
ggplot(aes(x)) +
geom_point(aes(y = binom1), color = "blue", size = 4) +
geom_point(aes(y = binom2), color = "red", size = 4)
PMF: \(\displaystyle f(x \mid \lambda)= \frac{ e^{-\lambda} \lambda^x}{x!}\); \(x = 0, 1, 2, \dots\) and \(\lambda > 0\)
Model for the counts of an event in a fixed period of time, with a rate of occurrence parameter \(\lambda\)
We write \(X \sim \text{Poisson}(\lambda)\)
tibble(x = 0:10) |>
mutate(y = dpois(x, 1)) |>
ggplot(aes(x, y)) +
geom_point(size = 4)
Recall that we can use the ECDF to estimate the true cumulative distribution function
z_ecdf <- ecdf(z)
z_ecdf(0) # should get close to 1/2[1] 0.503class(z_ecdf)[1] "ecdf" "stepfun" "function"normal_tbl <- tibble(z = sort(z)) |>
mutate(empirical = z_ecdf(z),
true = pnorm(z))
normal_tbl# A tibble: 1,000 × 3
z empirical true
<dbl> <dbl> <dbl>
1 -3.36 0.001 0.000383
2 -3.17 0.002 0.000757
3 -2.93 0.003 0.00169
4 -2.62 0.004 0.00436
5 -2.56 0.005 0.00520
6 -2.43 0.006 0.00758
7 -2.37 0.007 0.00881
8 -2.37 0.008 0.00901
9 -2.31 0.009 0.0104
10 -2.23 0.01 0.0128
# ℹ 990 more rowsnormal_tbl |>
pivot_longer(!z,
names_to = "method",
values_to = "val") |>
ggplot(aes(x = z, y = val, color = method)) +
geom_line()
n_points <- 10000
x <- runif(n_points, -1, 1)
y <- runif(n_points, -1, 1)
inside <- ifelse(x^2 + y^2 <= 1, 1, 0)
4 * mean(inside)[1] 3.1532tibble(x, y, inside) |>
ggplot(aes(x, y, color = factor(inside))) +
geom_point(show.legend = FALSE) +
coord_equal()
Random numbers generated in R (or any language) are not truly random; they are what we call pseudorandom
These are numbers generated by computer algorithms that very closely mimick truly random numbers
The default algorithm in R is called the Mersenne Twister
# suppose each group has 50 subjects
set.seed(100)
n_subjects <- 50
mean_drug <- 2
mean_nodrug <- runif(n_subjects, 0, 1)
x_drug <- 100 * rexp(n_subjects, 1 / mean_drug)
x_nodrug <- 100 * rexp(n_subjects, 1 / mean_nodrug)
tibble(x_drug, x_nodrug) |>
pivot_longer(everything(),
names_to = "group",
names_prefix = "x_",
values_to = "reduction") |>
ggplot(aes(x = reduction, y = after_stat(density),
color = group)) +
geom_histogram(aes(fill = group),
alpha = 0.5, color = "black",
position = "identity") +
geom_density(aes(color = group))
Not only computations take time, memory allocations do too
Changing the size of a vector takes just about as long as creating a new vector does
Each time the size changes, R needs to reconsider its allocation of memory to the object
Never reallocate a vector after each iteration
