```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
knitr::opts_chunk$set(tidy = TRUE)
```

### How to estimate \(\pi\) using Monte Carlo?
 
Introduce a square $\mathcal{S\subseteq}\mathbb{R}^{2}$, the sides being of length 2, with
inscribed disk $\mathcal{D}$ of radius 1. We are
interested in computing the area $I$ of $\mathcal{D}$ using Monte Carlo. We
have
\[\begin{aligned}
I  &  =\pi=\int\int_{\mathcal{D}}dx_{1}dx_{2}\\
&  =\int\int_{\mathcal{S}}\mathbb{I}_{\mathcal{D}}\left(  x_{1},x_{2}\right)
dx_{1}dx_{2}\text{ as }\mathcal{D\subset S}\\
&  =4\int\int_{\mathbb{R}^{2}}\mathbb{I}_{\mathcal{D}}\left(  x_{1}%
                                                             ,x_{2}\right)  \mu\left(  x_{1},x_{2}\right)  dx_{1}dx_{2}%
\end{aligned}
\]
where $\mathcal{S}:=\left[  -1,1\right]  \times\left[  -1,1\right]$ and
\[
 \mu\left(  x_{1},x_{2}\right)  =\frac{1}{4}\mathbb{I}_{\mathcal{S}}\left(
   x_{1},x_{2}\right)
 \]
is the uniform density on the square $\mathcal{S}$. In this case, we have
\[
 \widehat{I}_{n}=4\frac{n_{\mathcal{D}}}{n}%
 \]
where $n_{\mathcal{D}}$ is the number of samples which fell within the disk.

```{r first_sample}
# Asign basic parameters of the experiment.
numPoints <- 1e3 # number of samples
r <- 1 # radius of disk
C <- c(1,1) # center of disk and square
l <- 2 # size of square

# Draw 2-by-numPoints matrix of independent U[0, l] random variables.
U <- matrix(data=l * runif(n=2 * numPoints), nrow=numPoints, ncol=2)

# Compute index mask of all points which are inside the circle.
inside <- (U[,1]-C[1])^2+(U[,2]-C[2])^2 < r^2

# Count the number points in U that fell inside the disk.
M <- sum(inside)
cat("Monte Carlo estimate of pi: ", (M/numPoints) * (l^2) / (r^2))
```

```{r simple_plot, cache = TRUE}
# A simple graph of the points: green for the points inside the circle and black for the others.
plot(U[inside,1], U[inside,2], type = "p", xlab = "X", ylab = "Y", xlim = c(0, l), ylim = c(0,l), col = "green")
points(U[!inside,1], U[!inside,2], col = "black")
```

```{r graph_first_sample, tidy = TRUE, cache = TRUE}
# Graph using ggplot2.
library(ggplot2)
theme_set(theme_bw()) # customize the graph style

# Plot the points.
graph_sample <- qplot(x=U[,1], y=U[,2], colour = inside, geom = "point") # plot points
graph_sample <- graph_sample + xlab("X") + ylab("Y") # change axis labels

# Draw the circle.
graph_sample <- graph_sample + annotate("path",
              x=C[1]+r*cos(seq(0,2*pi,length.out=100)),
              y=C[2]+r*sin(seq(0,2*pi,length.out=100)), colour = "blue", size = 2)

# Draw the rectangle.
graph_sample <- graph_sample + geom_rect(aes(xmin=0,ymin=0,xmax=l,ymax=l), colour = "red", size = 2, fill = NA)

# Show the graph.
print(graph_sample)
```

```{r error, cache = TRUE}
# Repeat the above experiment 100 times.
numRepeats <- 100

# Store the results per experiment and per number of points used.
onlineErrors <- matrix(nrow=numPoints, ncol=numRepeats)
# For each experiment, 
for(i in 1:numRepeats){
  # draw a uniform sample,
  U <- matrix(data=l * runif(n=2 * numPoints), nrow=numPoints, ncol=2)
  # count points that fell inside the disk for each sample size up to numPoints, and
  M <- cumsum((U[,1]-C[1])^2+(U[,2]-C[2])^2 < r^2)
  # compute the relative error = (estimate - true value)/ true value.
  onlineEstimate <- M/(1:numPoints) * (l^2)/(r^2)
  onlineErrors[,i] <- (onlineEstimate - pi)/pi
}
# To plot the results nicely using ggplot2, first reshape the results into a data frame.
library(reshape2)
onlineErrors <- melt(onlineErrors)
names(onlineErrors) <- c("samples", "run", "error")
```

```{r graph, warning = FALSE, cache = TRUE, dependson = "error"}
# Plot error vs number of samples for each experiment.
graph_error <- ggplot(onlineErrors, aes(x = samples, y = error, colour = factor(run), group = factor(run)))
graph_error <- graph_error + geom_line(alpha = 1.0)
graph_error <- graph_error + theme(legend.position = "none") + ylim(-.3,.3)
print(graph_error)
```