Generating random numbers

Uniform distribution

In a uniform distribution, every number is equally likely. This is the “pick a number between 1 and 10” scenario, or rolling a single die. There are a couple ways to work with a uniform distribution in R: (1) sample() and (2) runif().

possible_answers <- c(1, 2, 3, 4, 5, 6)  # We could also write this as 1:6 instead

# Choose 1 random number
sample(possible_answers, size = 1)

## [1] 4

# Choose 3 random numbers
sample(possible_answers, size = 3)

## [1] 2 6 5

# Choose 10 random numbers, with replacement
sample(possible_answers, size = 10, replace = TRUE)

##  [1] 6 4 6 6 6 4 4 5 4 3

# Choose 8 numbers between 1 and 6, but don't replace them.
# This won't work!
sample(possible_answers, size = 8)

## Error in sample.int(length(x), size, replace, prob): cannot take a sample larger than the population when 'replace = FALSE'

set.seed(1234)
die <- tibble(value = sample(possible_answers, 
                             size = 1000, 
                             replace = TRUE))
die %>% 
  count(value)

## # A tibble: 6 x 2
##   value     n
##   <dbl> <int>
## 1     1   161
## 2     2   153
## 3     3   188
## 4     4   149
## 5     5   157
## 6     6   192

ggplot(die, aes(x = value)) +
  geom_bar() +
  labs(title = "1,000 rolls of a single die")

set.seed(1234)
die <- tibble(value = sample(possible_answers, 
                             size = 100000, 
                             replace = TRUE))

ggplot(die, aes(x = value)) +
  geom_bar() +
  labs(title = "100,000 rolls of a single die")

runif(5)

## [1] 0.09862 0.96294 0.88655 0.05623 0.44452

runif(5, min = 35, max = 56)

## [1] 46.83 42.89 37.75 53.22 46.13

# Generate 5 people between the ages of 18 and 35
round(runif(5, min = 18, max = 35), 0)

## [1] 21 28 33 34 31

set.seed(1234)
lots_of_numbers <- tibble(x = runif(5000, min = 18, max = 35))

ggplot(lots_of_numbers, aes(x = x)) +
  geom_histogram(binwidth = 1, color = "white", boundary = 18)

Normal distribution

The whole “choose a number between 1 and 10” idea of a uniform distribution is neat and conceptually makes sense, but most numbers that exist in the world tend to have higher probabilities around certain values—almost like gravity around a specific point. For instance, income in the United States is not uniformly distributed—a handful of people are really really rich, lots are very poor, and most are kind of clustered around an average.

The idea of having possible values clustered around an average is how the rest of these distributions work (uniform distributions don’t have any sort of central gravity point; all these others do). Each distribution is defined by different things called parameters, or values that determine the shape of the probabilities and locations of the clusters.

A super common type of distribution is the normal distribution. This is the famous “bell curve” you learn about in earlier statistics classes. A normal distribution has two parameters:

1. A mean (the center of the cluster)
2. A standard deviation (how much spread there is around the mean).

In R, you can generate random numbers from a normal distribution with the rnorm() function. It takes three arguments: the number of numbers you want to generate, the mean, and the standard deviation. It defaults to a mean of 0 and a standard deviation of 1, which means most numbers will cluster around 0, with a lot between −1 and 1, and some going up to −2 and 2 (technically 67% of numbers will be between −1 and 1, while 95% of numbers will be between −2–2ish)

rnorm(5)

## [1] -1.3662  0.5392 -1.3219 -0.2813 -2.1049

# Cluster around 10, with an SD of 4
rnorm(5, mean = 10, sd = 4)

## [1]  3.530  7.105 11.227 10.902 13.743

When working with uniform distributions, it’s easy to know how high or low your random values might go, since you specify a minimum and maximum number. With a normal distribution, you don’t specify starting and ending points—you specify a middle and a spread, so it’s harder to guess the whole range. Plotting random values is thus essential. Here’s 1,000 random numbers clustered around 10 with a standard deviation of 4:

set.seed(1234)

plot_data <- tibble(x = rnorm(1000, mean = 10, sd = 4))
head(plot_data)

## # A tibble: 6 x 1
##        x
##    <dbl>
## 1  5.17 
## 2 11.1  
## 3 14.3  
## 4  0.617
## 5 11.7  
## 6 12.0

ggplot(plot_data, aes(x = x)) +
  geom_histogram(binwidth = 1, boundary = 0, color = "white")

Neat. Most numbers are around 10; lots are between 5 and 15; some go as high as 25 and as low as −5.

Watch what happens if you change the standard deviation to 10 to make the spread wider:

set.seed(1234)

plot_data <- tibble(x = rnorm(1000, mean = 10, sd = 10))
head(plot_data)

## # A tibble: 6 x 1
##        x
##    <dbl>
## 1  -2.07
## 2  12.8 
## 3  20.8 
## 4 -13.5 
## 5  14.3 
## 6  15.1

ggplot(plot_data, aes(x = x)) +
  geom_histogram(binwidth = 1, boundary = 0, color = "white")

It’s still centered around 10, but now you get values as high as 40 and as low as −20. The data is more spread out now.

When simulating data, you’ll most often use a normal distribution just because it’s easy and lots of things follow that pattern in the real world. Incomes, ages, education, etc. all have a kind of gravity to them, and a normal distribution is a good way of showing that gravity. For instance, here are 1,000 simulated people with reasonable random incomes, ages, and years of education:

set.seed(1234)

fake_people <- tibble(income = rnorm(1000, mean = 40000, sd = 15000),
                      age = rnorm(1000, mean = 25, sd = 8),
                      education = rnorm(1000, mean = 16, sd = 4))
head(fake_people)

## # A tibble: 6 x 3
##   income   age education
##    <dbl> <dbl>     <dbl>
## 1 21894. 15.4      12.1 
## 2 44161. 27.4      15.6 
## 3 56267. 12.7      15.6 
## 4  4815. 30.1      20.8 
## 5 46437. 30.6       9.38
## 6 47591.  9.75     11.8

fake_income <- ggplot(fake_people, aes(x = income)) +
  geom_histogram(binwidth = 5000, color = "white", boundary = 0) +
  labs(title = "Simulated income")

fake_age <- ggplot(fake_people, aes(x = age)) +
  geom_histogram(binwidth = 2, color = "white", boundary = 0) +
  labs(title = "Simulated age")

fake_education <- ggplot(fake_people, aes(x = education)) +
  geom_histogram(binwidth = 2, color = "white", boundary = 0) +
  labs(title = "Simulated education")

fake_income + fake_age + fake_education

These three columns all have different centers and spreads. Income is centered around $45,000, going up to almost $100,000 and as low as −$10,000; age is centered around 25, going as low as 0 and as high as 50; education is centered around 16, going as low as 3 and as high as 28. Cool.

Again, when generating these numbers, it’s really hard to know how high or low these ranges will be, so it’s a good idea to plot them constantly. I settled on sd = 4 for education only because I tried things like 1 and 10 and got wild looking values (everyone basically at 16 with little variation, or everyone ranging from −20 to 50, which makes no sense when thinking about years of education). Really it’s just a process of trial and error until the data looks good and reasonable.

Truncated normal distribution

Sometimes you’ll end up with negative numbers that make no sense. Look at income in the plot above, for instance. Some people are earning −$10,000 year. The rest of the distribution looks okay, but those negative values are annoying.

To fix this, you can use something called a truncated normal distribution, which lets you specify a mean and standard deviation, just like a regular normal distribution, but also lets you specify a minimum and/or maximum so you don’t get values that go too high or too low.

R doesn’t have a truncated normal function built-in, but you can install the truncnorm package and use the rtruncnorm() function. A truncated normal distribution has four parameters:

1. A mean (mean)
2. A standard deviation (sd)
3. A minimum (optional) (a)
4. A maximum (optional) (b)

For instance, let’s pretend you have a youth program designed to target people who are between 12 and 21 years old, with most around 14. You can generate numbers with a mean of 14 and a standard deviation of 5, but you’ll create people who are too old, too young, or even negatively aged!

set.seed(1234)

plot_data <- tibble(fake_age = rnorm(1000, mean = 14, sd = 5))
head(plot_data)

## # A tibble: 6 x 1
##   fake_age
##      <dbl>
## 1     7.96
## 2    15.4 
## 3    19.4 
## 4     2.27
## 5    16.1 
## 6    16.5

ggplot(plot_data, aes(x = fake_age)) +
  geom_histogram(binwidth = 2, color = "white", boundary = 0)

To fix this, truncate the range at 12 and 21:

library(truncnorm)  # For rtruncnorm()

set.seed(1234)

plot_data <- tibble(fake_age = rtruncnorm(1000, mean = 14, sd = 5, a = 12, b = 21))
head(plot_data)

## # A tibble: 6 x 1
##   fake_age
##      <dbl>
## 1     15.4
## 2     19.4
## 3     16.1
## 4     16.5
## 5     14.3
## 6     18.8

ggplot(plot_data, aes(x = fake_age)) +
  geom_histogram(binwidth = 1, color = "white", boundary = 0)

And voila! A bunch of people between 12 and 21, with most around 14, with no invalid values.

Beta distribution

Normal distributions are neat, but they’re symmetrical around the mean (unless you truncate them). What if your program involves a test with a maximum of 100 points where most people score around 85, but a sizable portion score below that. In other words, it’s not centered at 85, but is skewed left.

To simulate this kind of distribution, we can use a Beta distribution. Beta distributions are neat because they naturally only range between 0 and 1—they’re perfect for things like percentages or proportions or or 100-based exams.

Unlike a normal distribution, where you use the mean and standard deviation as parameters, Beta distributions take two non-intuitive parameters:

1. shape1
2. shape2

What the heck are these shapes though?! This answer at Cross Validated does an excellent job of explaining the intuition behind Beta distributions and it’d be worth it to read it.

Basically, Beta distributions are good at modeling probabilities of things, and shape1 and shape2 represent specific parts of a probability formula.

Let’s say that there’s an exam with 10 points where most people score a 6/10. Another way to think about this is that an exam is a collection of correct answers and incorrect answers, and that the percent correct follows this equation:

\[ \frac{\text{Number correct}}{\text{Number correct} + \text{Number incorrect}} \]

If you scored a 6, you could write that as:

\[ \frac{6}{6 + 4} \]

To make it more general, we can use Greek variable names: \(\alpha\) for the number correct and \(\beta\) for the number incorrect, leaving us with this:

\[ \frac{\alpha}{\alpha + \beta} \]

Neat.

In a Beta distribution, the \(\alpha\) and \(\beta\) in that equation correspond to shape1 and shape2. If we want to generate random scores for this test where most people get 6/10, we can use rbeta():

set.seed(1234)

plot_data <- tibble(exam_score = rbeta(1000, shape1 = 6, shape2 = 4)) %>% 
  # rbeta() generates numbers between 0 and 1, so multiply everything by 10 to
  # scale up the exam scores
  mutate(exam_score = exam_score * 10)

ggplot(plot_data, aes(x = exam_score)) +
  geom_histogram(binwidth = 1, color = "white") +
  scale_x_continuous(breaks = 0:10)

Most people score around 6, with a bunch at 5 and 7, and fewer in the tails. Importantly, it’s not centered at 6—the distribution is asymmetric.

The magic of—and most confusing part about—Beta distributions is that you can get all sorts of curves by just changing the shape parameters. To make this easier to see, we can make a bunch of different Beta distributions. Instead of plotting them with histograms, we’ll use density plots (and instead of generating random numbers, we’ll plot the actual full range of the distribution (that’s what dbeta and geom_function() do in all these examples)).

Here’s what we saw before, with \(\alpha\) (shape1) = 6 and \(\beta\) (shape2) = 4:

ggplot() +
  geom_function(fun = dbeta, args = list(shape1 = 6, shape2 = 4))

Again, there’s a peak at 0.6 (or 6), which is what we expected.

We can make the distribution narrower if we scale the shapes up. Here pretty much everyone scores around 50% and 75%.

ggplot() +
  geom_function(fun = dbeta, args = list(shape1 = 60, shape2 = 40))

So far all these curves look like normal distributions, just slightly skewed. But when if most people score 90–100%? Or most fail? A Beta distribution can handle that too:

ggplot() +
  geom_function(fun = dbeta, args = list(shape1 = 9, shape2 = 1), color = "blue") +
  geom_function(fun = dbeta, args = list(shape1 = 1, shape2 = 9), color = "red")

With shape1 = 9 and shape2 = 1 (or \(\frac{9}{9 + 1}\)) we get most around 90%, while shape1 = 1 and shape2 = 9 (or \(\frac{1}{1 + 9}\)) gets us most around 10%.

Check out all these other shapes too:

ggplot() +
  geom_function(fun = dbeta, args = list(shape1 = 5, shape2 = 5), color = "blue") +
  geom_function(fun = dbeta, args = list(shape1 = 2, shape2 = 5), color = "red") +
  geom_function(fun = dbeta, args = list(shape1 = 80, shape2 = 23), color = "orange") +
  geom_function(fun = dbeta, args = list(shape1 = 13, shape2 = 17), color = "brown")

In real life, if I don’t want to figure out the math behind the \(\frac{\alpha}{\alpha + \beta}\) shape values, I end up just choosing different numbers until it looks like the shape I want, and then I use rbeta() with those parameter values. Like, how about we generate some numbers based on the red line above, with shape1 = 2 and shape2 = 5, which looks like it should be centered around 0.2ish (\(\frac{2}{2 + 5} = 0.2857\)):

set.seed(1234)

plot_data <- tibble(thing = rbeta(1000, shape1 = 2, shape2 = 5)) %>% 
  mutate(thing = thing * 100)
head(plot_data)

## # A tibble: 6 x 1
##   thing
##   <dbl>
## 1 10.1 
## 2 34.5 
## 3 55.3 
## 4  2.19
## 5 38.0 
## 6 39.9

ggplot(plot_data, aes(x = thing)) +
  geom_histogram(binwidth = 2, color = "white", boundary = 0)

It worked! Most values are around 20ish, but some go up to 60–80.

Binomial distribution

Often you’ll want to generate a column that only has two values: yes/no, treated/untreated, before/after, big/small, red/blue, etc. You’ll also likely want to control the proportions (25% treated, 62% blue, etc.). You can do this in two different ways: (1) sample() and (2) rbinom().

set.seed(1234)

# Choose 5 random T/F values
possible_things <- c(TRUE, FALSE)
sample(possible_things, 5, replace = TRUE)

## [1] FALSE FALSE FALSE FALSE  TRUE

set.seed(1234)
candidates <- c("Person 1", "Person 2")
sample(candidates, size = 1, prob = c(0.8, 0.2))

## [1] "Person 1"

set.seed(1234)
fake_elections <- tibble(winner = sample(candidates, 
                                         size = 1000, 
                                         prob = c(0.8, 0.2),
                                         replace = TRUE))
fake_elections %>% 
  count(winner)

## # A tibble: 2 x 2
##   winner       n
##   <chr>    <int>
## 1 Person 1   792
## 2 Person 2   208

ggplot(fake_elections, aes(x = winner)) +
  geom_bar()

set.seed(1234)

rbinom(5, size = 20, prob = 0.6)

## [1] 15 11 11 11 10

set.seed(1234)

rbinom(5, size = 1, prob = 0.6)

## [1] 1 0 0 0 0

set.seed(12345)

plot_data <- tibble(thing = rbinom(2000, 1, prob = 0.6)) %>% 
  # Make this a factor since it's basically a yes/no categorical variable
  mutate(thing = factor(thing))

plot_data %>% 
  count(thing) %>% 
  mutate(proportion = n / sum(n))

## # A tibble: 2 x 3
##   thing     n proportion
##   <fct> <int>      <dbl>
## 1 0       840      0.42 
## 2 1      1160      0.580

ggplot(plot_data, aes(x = thing)) +
  geom_bar()

Poisson distribution

One last common distribution that you might find helpful when simulating data is the Poisson distribution (in French, “poisson” = fish, but here it’s not actually named after the animal, but after French mathematician Siméon Denis Poisson).

A Poisson distribution is special because it generates whole numbers (i.e. nothing like 1.432) that follow a skewed pattern (i.e. more smaller values than larger values). There’s all sorts of fancy math behind it that you don’t need to worry about so much—all you need to know is that it’s good at modeling things called Poisson processes.

For instance, let’s say you’re sitting at the front door of a coffee shop (in pre-COVID days) and you count how many people are in each arriving group. You’ll see something like this:

• 1 person
• 1 person
• 2 people
• 1 person
• 3 people
• 2 people
• 1 person

Lots of groups of one, some groups of two, fewer groups of three, and so on. That’s a Poisson process: a bunch of independent random events that combine into grouped events.

That sounds weird and esoteric (and it is!), but it reflects lots of real world phenomena, and things you’ll potentially want to measure in a program. For instance, the number of kids a family has follows a type of Poisson process. Lots of families have 1, some have 2, fewer have 3, even fewer have 4, and so on. The number of cars in traffic, the number of phone calls received by an office, arrival times in a line, and even the outbreak of wars are all examples of Poisson processes.

You can generate numbers from a Poisson distribution with the rpois() function in R. This distribution only takes a single parameter:

• lambda (\(\lambda\))

The \(\lambda\) value controls the rate or speed that a Poisson process increases (i.e. jumps from 1 to 2, from 2 to 3, from 3 to 4, etc.). I have absolutely zero mathematical intuition for how it works. The two shape parameters for a Beta distribution at least fit in a fraction and you can wrap your head around that, but the lambda in a Poisson distribution is just a mystery to me. So whenever I use a Poisson distribution for something, I just play with the lambda until the data looks reasonable.

Let’s assume that the number of kids a family has follows a Poisson process. Here’s how we can use rpois() to generate that data:

set.seed(123)

# 10 different families
rpois(10, lambda = 1)

##  [1] 0 2 1 2 3 0 1 2 1 1

Cool. Most families have 0–1 kids; some have 2; one has 3.

It’s easier to see these patterns with a plot:

set.seed(1234)

plot_data <- tibble(num_kids = rpois(500, lambda = 1))
head(plot_data)

## # A tibble: 6 x 1
##   num_kids
##      <int>
## 1        0
## 2        1
## 3        1
## 4        1
## 5        2
## 6        1

plot_data %>% 
  group_by(num_kids) %>% 
  summarize(count = n()) %>% 
  mutate(proportion = count / sum(count))

## # A tibble: 6 x 3
##   num_kids count proportion
##      <int> <int>      <dbl>
## 1        0   180      0.36 
## 2        1   187      0.374
## 3        2    87      0.174
## 4        3    32      0.064
## 5        4    11      0.022
## 6        5     3      0.006

ggplot(plot_data, aes(x = num_kids)) +
  geom_bar()

Here 75ish% of families have 0–1 kids (36% + 37.4%), 17% have 2 kids, 6% have 3, 2% have 4, and only 0.6% have 5.

We can play with the \(\lambda\) to increase the rate of kids per family:

set.seed(1234)

plot_data <- tibble(num_kids = rpois(500, lambda = 2))
head(plot_data)

## # A tibble: 6 x 1
##   num_kids
##      <int>
## 1        0
## 2        2
## 3        2
## 4        2
## 5        4
## 6        2

plot_data %>% 
  group_by(num_kids) %>% 
  summarize(count = n()) %>% 
  mutate(proportion = count / sum(count))

## # A tibble: 8 x 3
##   num_kids count proportion
##      <int> <int>      <dbl>
## 1        0    62      0.124
## 2        1   135      0.27 
## 3        2   145      0.290
## 4        3    88      0.176
## 5        4    38      0.076
## 6        5    19      0.038
## 7        6    10      0.02 
## 8        7     3      0.006

ggplot(plot_data, aes(x = num_kids)) +
  geom_bar()

Now most families have 1–2 kids. Cool.