This simple Monte Carlo simulation aims to estimate the probability of obtaining a sum of 10 when rolling three fair dice, returning the probability estimate.
In this analysis, we explore the influence of two factors: the number of rolls and the number of estimates per simulation on the resulting probabilities.
#function to simulate the roll of three fair dice, returning the sum
rollDie <- function() {
die <- sample(1:6, 3)
return(die)
}
A function to estimate the probability of a sum of 10 when rolling three fair dice,returning that probability estimate. Include a parameter that indicates the number of rolls used in producing the estimate, defaulting to 1000.
oneEstimate <- function(num_rolls = 1000) {
# this function returns the estimate of a sum of 10
# when rolling three fair dice
count <- 0
for (i in 1:num_rolls) {
if (sum(rollDie())== 10) {
count <- count + 1
}
}
return(count / num_rolls)
}
A function that will create and return a vector of probability estimates, where each estimate is #the result of calling your second function above. The current function should have a parameter that indicates the number of estimates in the vector (defaulting to 100) and the number of rolls used per estimate (defaulting to 1000). Return the vector of estimates at the end
getEstimates <- function(num_estimates = 100, num_rolls = 1000)
{
# this function creates a vector of size num_estimates and fills that vector with estimates of rolling a sum of 10 using
# num_rolls rolls per computed estimate
estimates <- rep(0, num_estimates)
for (i in 1:num_estimates) {
estimates[i] <- oneEstimate(num_rolls)
}
return(estimates)
}
A separate function that, given a vector of probability estimates, will produce a histogram of those estimates, and will superimpose vertical lines for the mean and for the mean +/- two standard deviations.
oneHistogram <- function(num_rolls = 1000)
{
all_estimates <- getEstimates(100, num_rolls)
hist(all_estimates, xlim = c(0.05, 0.3))
m <- mean(all_estimates)
abline(v = m, col = "blue", lwd = 2)
s <- sd(all_estimates)
abline(v = c(m - 2*s, m + 2*s), col = "red", lwd = 2, lty = "dashed")
}
oneHistogram(100)
Fixing the number of estimates at 1000 and changing number of rolls to 10, 1000, and 1000
par(mfrow = c(2,3)) # display 2 rows, 3 columns
oneHistogram <- function(num_estimates = 1000, num_rolls = 1000)
{
all_estimates <- getEstimates(num_estimates, num_rolls)
hist(all_estimates, xlim = c(0, 0.5), main = "")
title(paste(num_rolls, "rolls", num_estimates, "estimates"))
m <- mean(all_estimates)
abline(v = m, col = "blue", lwd = 2)
s <- sd(all_estimates)
abline(v = c(m - 2*s, m + 2*s), col = "red", lwd = 2, lty = "dashed")
print(quantile(all_estimates, c(0.5, 0.75, 0.95, 0.99)))
}
oneHistogram(num_rolls = 10)
## 50% 75% 95% 99%
## 0.1 0.2 0.3 0.4
oneHistogram(num_rolls = 100)
## 50% 75% 95% 99%
## 0.15 0.17 0.21 0.24
oneHistogram(num_rolls = 1000)
## 50% 75% 95% 99%
## 0.149 0.157 0.168 0.175
#Fixing the number of rolls at 1000 and changing number of estimates to 10, 1000, and 1000
oneHistogram(num_estimates = 10)
## 50% 75% 95% 99%
## 0.14500 0.15425 0.15765 0.15873
oneHistogram(num_estimates = 100)
## 50% 75% 95% 99%
## 0.152 0.158 0.167 0.172
oneHistogram(num_estimates = 1000)
## 50% 75% 95% 99%
## 0.150 0.158 0.169 0.176
From the histograms, it becomes evident that the number of rolls has a significant impact on the overall standard deviation of the probability. The fewer the number of rolls, the higher the variance in probability, showing that a higher number of rolls leads to higher accuracy in the data.
On the second row of the histograms, we observe that the number estimates per simulation has a relatively minor impact on the accuracy of the probabilities since the standard deviations do not greatly differ. Nevertheless, increasing the number of estimates tends to reduce the standard deviation in probability values, indicating a trend toward improved precision with a larger number of estimates.