-
Notifications
You must be signed in to change notification settings - Fork 0
/
Pi Day.Rmd
173 lines (134 loc) · 4.93 KB
/
Pi Day.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
---
title: "Celebrating Pi Day <br> with R!"
author: "Alyssa Columbus"
date: "March 14, 2024"
output:
ioslides_presentation:
css: styles.css
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE)
```
## Overview
Monte Carlo simulations are powerful tools that leverage randomness to model various statistical phenomena. In honor of Pi ($\pi$) Day, we will use this type of method to estimate the value of $\pi$.
```{r echo=TRUE}
pi
```
```{r echo=FALSE, warning=FALSE, message=FALSE}
# Import necessary libraries
library(knitr)
library(grid)
library(plotrix)
# Define helper functions
is_in_circle <- function(x,y) {
sqrt((x-0.5)^2 + (y-0.5)^2) <= 0.5
}
estimate_pi <- function() {
4*length(darts[darts])/num_attempts
}
throw_darts <- function(n) {
num_attempts <<- n
dart_positions <<- matrix(runif(2*num_attempts), ncol=2)
darts <<- is_in_circle(dart_positions[,1], dart_positions[,2])
}
plot_darts <- function() {
plot(c(0, 1), c(0,1), type = "n", asp = 1, main="The Dartboard",xlab='Unit Square/Circle',ylab='',yaxt='n',frame.plot = FALSE)
rect(0,0,1,1,border='black')
draw.circle(0.5,0.5,0.5,border='black')
points(dart_positions[!darts,1], dart_positions[!darts,2], col='grey', pch=20)
points(dart_positions[darts,1], dart_positions[darts,2], col='purple', pch=20)
}
calculate_pi <- function(n) {
throw_darts(n)
estimate_pi()
}
```
## Estimating $\pi$ with Monte Carlo Simulation
Imagine you’re engaged in a friendly game of darts, where a circular dartboard is mounted on a square piece of wood of the same diameter. Your friend claims that they can estimate the number $\pi$ by randomly throwing darts at the board.
You’re skeptical, but you decide to see how it plays out!
## Dart Throwing in Action -- 100 Darts!
```{r}
throw_darts(100)
plot_darts()
paste0("Pi Estimate: ", estimate_pi())
```
## Dart Throwing in Action -- 10,000 Darts!
```{r}
throw_darts(10000)
plot_darts()
paste0("Pi Estimate: ", estimate_pi())
```
## Dart Throwing in Action -- 30,000 Darts!
```{r}
throw_darts(30000)
plot_darts()
paste0("Pi Estimate: ", estimate_pi())
```
## Understanding the Circle-Square Area Relationship
As the darts fly, you notice your friend tallying the darts landing inside vs. outside of the circle (but still within the square). They keep mentioning numbers that seem to gradually get closer and closer to $\pi$.
To understand what your friend is doing, let’s examine the relationship between the **unit circle** and the **unit square**. We can start with the familiar formula for the area of a (unit) circle, $\pi r^2$. Similarly, we can find that the area of a unit square is $(2r)^2$.
## The Ratio of the Area of the Circle to the Area of the Square
Since the area of the circle is $\pi r^2$ and the area of the square is $(2r)^2$, the ratio of their areas can be calculated as follows:
$$\pi r^2 / (2r)^2 = \pi r^2 / 4 r^2 = \pi/4.$$
```{r echo = TRUE}
pi/4
```
If we take this ratio and multiply it by 4, we get an estimate of $\pi$!
## Estimating $\pi$
The more darts your friend throws, the closer their estimate gets to the actual value of $\pi$.
```{r warning=FALSE}
plot(c(0,3000),c(-0.5,5), type = "n", main=NULL, xlab='Number of Dart Throws',
ylab='Estimate of Pi', frame.plot = FALSE)
abline(h=pi,col=2,lty=2)
lines(sapply(1:3000,calculate_pi))
abline(h=pi,col=2,lty=1)
```
## Let's Join the Dart Throwing!
```{r echo=TRUE,eval=FALSE}
# Import necessary libraries
library(grid)
library(plotrix)
# Define helper functions
is_in_circle <- function(x,y) {
sqrt((x-0.5)^2 + (y-0.5)^2) <= 0.5
}
estimate_pi <- function() {
4*length(darts[darts])/num_attempts
}
throw_darts <- function(n) {
num_attempts <<- n
dart_positions <<- matrix(runif(2*num_attempts), ncol=2)
darts <<- is_in_circle(dart_positions[,1], dart_positions[,2])
}
```
## Let's Join the Dart Throwing!
```{r echo=TRUE,eval=FALSE}
plot_darts <- function() {
plot(c(0, 1), c(0,1), type = "n", asp = 1,
main="The Dartboard",
xlab='Unit Square/Circle',
ylab='',yaxt='n', frame.plot = FALSE)
rect(0,0,1,1,border='black')
draw.circle(0.5,0.5,0.5,border='black')
points(dart_positions[!darts,1],
dart_positions[!darts,2], col='grey', pch=20)
points(dart_positions[darts,1],
dart_positions[darts,2], col='purple', pch=20)
}
calculate_pi <- function(n) {
throw_darts(n)
estimate_pi()
}
```
## Let's Throw $\pi \times 10,000$ Darts!
```{r echo=TRUE, eval=FALSE}
throw_darts(31416); plot_darts(); estimate_pi()
```
## Let's Throw $\pi \times 10,000$ Darts!
```{r}
throw_darts(31416); plot_darts(); estimate_pi()
```
## Thank You!
**Enjoy your Pi Day!**
Demonstration inspired by [this presentation by Jamin Ragle](https://rpubs.com/jaminragle/97163).
Originally presented at [R-Ladies Irvine](https://rladiesirvine.org).