```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE) knitr::opts_chunk$set(tidy = TRUE) ``` ### Box-Muller algorithm As described in Lecture 2. ```{r function} box_muller <- function(size) { # Since Box-Muller generates variates in pairs, we round up from the requested size. num_pairs <- ceiling(size / 2) R <- sqrt(-2 * log(runif(num_pairs))) theta <- 2 * pi * runif(num_pairs) X <- R * cos(theta) Y <- R * sin(theta) c(X, Y)[0:size] } ``` Let's test the function. First, run it. ```{r run} samples <- box_muller(10001) ``` Do we have the right number of samples? ```{r length} cat("length(samples): ", length(samples)) ``` Plot a histogram of the samples. ```{r histogram} library(ggplot2) theme_set(theme_bw()) graph_hist <- qplot(x = samples, geom = "blank") + geom_histogram(aes(y = ..density..), binwidth = 0.1) # Superimpose a line showing the correct density. xs <- seq(min(samples), max(samples), 0.05) ys <- dnorm(xs) df <- data.frame(xs, ys) graph_hist <- graph_hist + geom_line(aes(x = xs, y = ys), df, colour = "red") graph_hist ``` Finally, perform a Kolmogorov-Smirnov test. ```{r ks} ks.test(samples, "pnorm") ```