\documentclass[12pt]{article} \SweaveOpts{echo=T,pdf=T,eps=F,eval=T,keep.source=T} \usepackage{Sweave} \usepackage{mathpazo} \usepackage{hyperref,url} \usepackage[a4paper,margin=1.5cm]{geometry} \begin{document} \title{Example Sweave document: estimating $\pi$} \author{Stephen J Eglen} \date{\today} \maketitle \section{Introduction} This is an example document created using the Sweave system (\url{http://www.statistik.lmu.de/~leisch/Sweave/}). Sweave is a tool for combining both \LaTeX\ documentation and R code within the same file. For this document, the master file is \url{estimate.Rnw}. This is processed by the Sweave system in R, which runs the R code to generate textual/graphical output, and also creates a \LaTeX\ document. The \LaTeX\ document is then typeset to create the pdf document. On unix/macintosh, the following commands should recreate the pdf file: \begin{verbatim} $R CMD Sweave estimate.Rnw$ pdflatex estimate.tex \end{verbatim} \section{Task: estimate the value of $\pi$} Our task is to estimate the value of $\pi$ by simulating darts being thrown at a dartboard. Imagine that the person throwing the darts is not very good, and randomly throws each dart so that it falls uniformly within a square of side length $2r$, with the dartboard of radius $r$ centred within that square. If the player throws $n$ darts, and $d$ of them hit the dartboard, then for large enough $n$, the ratio $d/n$ should approximate the ratio of the area of the dartboard to the enclosing square, $\pi r^2 / 4 r^2 \equiv \pi/4$. From this, we can estimate $\pi \approx 4d/n$. We start with an example, using R to draw both the dartboard and the surrounding square, together with $n=50$ darts. The radius of the dartboard here is 1 unit, although the value is not important. \setkeys{Gin}{width=0.6\textwidth} \begin{center} <>= r <- 1 n <- 50 par(las=1) plot(NA, xlim=c(-r,r), ylim=c(-r,r), asp=1, bty='n', xaxt='n', yaxt='n', xlab='', ylab='') axis(1, at=c(-r,0,r)); axis(2, at=c(-r,0,r)) symbols(x=0, y=0, circles=r, inch=F, add=T) x <- runif(n, -r, r); y <- runif(n, -r, r) inside <- x^2 + y^2 < r^2 d <- length(which(inside)) points(x, y, pch=ifelse(inside, 19, 4)) rect(-r, -r, r, r, border='blue', lwd=2) @ \end{center} A dart is drawn as a filled circle if it falls within the dartboard, else it is drawn as a cross. In this case the number of darts within the circle is \Sexpr{d}, and so the estimated value is $\pi \approx \Sexpr{4*d/n}$. The estimate of $\pi$ should improve as we increase the number of darts thrown at the dartboard. To verify this, we write a short function that, given the number of darts to throw, $n$, returns an estimate of $\pi$. <<>>= estimate.pi <- function(n=1000) { ## Return an estimate of PI using dartboard method ## with N trials. r <- 1 # radius of dartboard x <- runif(n, min=-r, max=r) y <- runif(n, min=-r, max=r) l <- sqrt(x^2 + y^2) d <- length(which(l>= replicate(9, estimate.pi()) @ Finally, for a given value of $n$, we can show 99 estimates of $\pi$, as clearly the estimate will vary from run to run. In the following plot, we compare the estimates of $\pi$ for three different values of $n$: \begin{center} <>= ns <- 10^c(2,3,4) res <- lapply(ns, function(n) replicate(99, estimate.pi(n))) par(las=1, bty='n') stripchart(res, method="jitter", group.names=ns, xlab="number of darts", ylab=expression(paste('estimate of ', pi)), vert=T, pch=20, cex=0.5) abline(h=pi, col='red') @ \end{center} As the number of darts increases, the estimate of $\pi$ gradually converges onto the actual value of $\pi$ (shown by the solid red line). \end{document}