Problem Set - 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 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$. You will write a function called runif_zx81() that generates uniform draws in the same way as the Z81.
1. Type out the equation for the Lehmer random number sequence in your RMarkdown/Quarto document, using LaTeX syntax to typeset mathematical notation. (E.g. $x^2$ produces $x^2$.) Then define the variables m and a in R along with seed <- 42. Now write R code that prints out the next element of the sequence.
2. Now imagine you have done this five times, and you are left with the sequence vector Lehmer_seq <- c(42, 3150, 39639, 23760, 12501, 20057). How can you transform this vector to live in a $[0, 1]$ interval? What about a $[3, 5]$ interval?
3. Use the code from the preceding part as a starting point to write the function runif_zx81(). 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 and max – should be identical to the corresponding arguments of the base R function runif(). Set min = 0, max = 1 as default values. Your function should look something like the pseudo code below.
4. Choose a seed and make n=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. Hint: use the tibble() function to transform your random draws into the required input for ggplot2.

# Pseudocode for 1)c)
runif_zx81 <- function(seed, n, min = 0, max = 1) {
  
  # Set the a and m parameters as specified.
  
  # Add warning messages in case the seed input is negative or larger than m.
  
  # Initialise an empty vector to save the draws, and save the start of the 
  # sequence as the first draw. Hint: how do you access elements of vectors?
  
  # Run a for loop to construct as many elements of the sequence as specified by
  # the number of draws n.
  
  # Adjust the interval of your vector to run from min to max.
  
  # Return the vector of pseudorandom numbers.

}

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 uses runif_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 to runif_zx81(). The remaining arguments –n, mean, and sd– should match those of the base R function rnorm().
1. Think about the output of runif_zx81(). It returns one vector, but the Box-Muller Transform requires two independent standard uniform vectors, $U_1$ and $U_2$. The transform also gives you two vectors, $Z_1$ and $Z_2$, but rnorm_zx81() should only return one vector combining both. Write some bullet points of pseudo code for what rnorm_zx81(seed, n, mean, sd) does. How will it combine runif_zx81() and the Box-Muller Transform to produce the desired output?
2. Now suppose I give you a sequence of uniform draws to start off with, namely unif_seq <- c(0.5600805, 0.5767570, 0.8858708, 0.9313472, 0.7665961, 0.9763004). Write code to implement the Box-Mueller transform to convert this particular vector of six uniform draws into a vector of six normal draws with mean 2 and standard deviation 0.5.
3. Using what you’ve learned from the preceding parts, implement the function rnorm_zx81(). Hint: think carefully about how to ensure that your function still works if the user requests an odd number of random normal draws n.
4. Choose a seed and generate 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.
Lest you take the wrong lesson away from this problem set, it would be remiss of me not to point out that the pseudo-random number generator you have implemented here is seriously flawed. To see the problem, start by making it a bit starker: write a new function that does exactly the same thing as runif_zx81() but sets $a = 66$ and $m = 401$. Then set a seed and generate 1000 supposedly iid Uniform(0,1) draws. Finally, make a make a 2D scatter plot of the pairs $(U_1, U_2), (U_2, U_3), (U_3, U_4), \dots, (U_{999}, U_{1000})$. Do you notice anything strange? Comment briefly, perhaps after peeking at this classic article by (Marsaglia, 1968).