<- function(seed, n, min = 0, max = 1) {
runif_zx81 <- 2^16 + 1
m if(seed > m) return('Error: seed cannot exceed 65537.')
if(seed <= 0) return('Error: must be positive.')
<- 75
a <- vector(mode = 'numeric', length = n+1)
x 1] <- seed
x[for(i in 1:n) {
+ 1] <- (a * x[i]) %% m
x[i
}+ (x[-1] / m) * (max - min)
min }
Exercise - Monte Carlo On a Desert Island
When it comes to basic numerical methods, “rolling your own” is generally a bad idea. If you need to generate (pseudo)random numbers in everyday life, you should rely on the excellent work of generations of numerical analysts and computer scientists and just use the base R functions runif()
, rnorm()
etc. But what if you were stranded on a desert island and needed to re-create civilization from scratch? Surely you would need a way to simulate uniform and normal random variables! In this exercise you will generate pseudorandom numbers as though you were Robinson Crusoe. Again: unless civilization collapses, you you should not do this in real life. Nevertheless, this exercise will give you a better appreciation for how much work goes on behind the scenes when you run a simulation study!
- The ZX81 was a personal computer sold in the UK in the early 1980s. To simulate standard uniform draws, the ZX81 relied on the Lehmer random number generator discussed in our lecture on Monte Carlo simulation. In particular, it set \(m = 2^{16} + 1\) and \(a = 75\).
- Write a function called
runif_zx81()
that generates uniform draws in the same way as the Z81. Your function should take four arguments. The first argument,seed
, should be the seed of the Lehmer random number generator. The remaining arguments–n
,min
andmax
–should be identical to the corresponding arguments of the base R functionrunif()
. - Choose a seed and make 1000 random draws from
runif_zx81()
. Do these appear to be iid Uniform(0,1) draws? Make a histogram, QQ plot, and time series plot to check.
- Write a function called
- Let \(U_1\) and \(U_2\) be a pair of independent standard uniform random variables. Now define two new random variables: \(R \equiv \sqrt{-2 \log(U_1)}\) and \(\Theta \equiv 2\pi U_2\). Then it can be shown that \(Z_1 \equiv R \cos(\Theta)\) and \(Z_2 \equiv R \sin(\Theta)\) are a pair of independent standard normal random variables. This result is usually called the Box-Muller Transform.
- Write a function called
rnorm_zx81()
that usesrunif_zx81()
along with the the Box-Muller transform to simulate draws from a normal random variable. The first argument,seed
, should be the seed that you pass torunif_zx81()
. The remaining arguments–n
,mean
, andsd
–should match those of the base R functionrnorm()
. - Choose a seed and make 1000 random draws from
rnorm_zx81()
. Do these appear to be iid standard normal draws? Make a histogram, QQ plot, and time series plot to check.
- Write a function called
Solutions
1(a)
1(b)
library(tidyverse)
<- runif_zx81(1983, 1000)
u tibble(u) |>
ggplot(aes(u)) +
geom_histogram(bins = 20)
tibble(u) |>
ggplot(aes(1:length(u), u)) +
geom_line()