#
# CHAPTER 25: Simulation
#
# Subtopic 1: prime modulus multiplicative
#
runif(14)
#
# linear congruential generators (Lehmer, 1951)
#
# x_i = (a * x_{i - 1}) mod m
#
# for carefully selected positive integers a and m and
# seed x_0. The "modulus" m is typically big and prime.
#
# The random numbers are u_i = x_i / m
#
#########################################################
#
# Example 1: a = 5, m = 7, x_0 = 3
#
# x_i = (5 * x_{i - 1}) mod 7
#
# 3, 1, 5, 4, 6, 2, 3, ...
# 1/7, 5/7, 4/7, 6/7, 2/7, 3/7, ...
#
#########################################################
#
# Example 2: a = 7, m = 61, x_0 = 3
# Generate 14 pseudorandom numbers
# x_i = (7 * x_{i - 1}) mod 61
#
x = numeric(14)
x[1] = 3
for (i in 2:14) x[i] = (7 * x[i - 1]) %% 61
u = x / 61
x
u
#
# Example 3: make a plot of overlapping pairs from the
# generator with a = 7, m = 61, x_0 = 3
#
x = numeric(60)
x[1] = 3
for (i in 2:60) x[i] = (7 * x[i - 1]) %% 61
u = x / 61
plot(0, 0, type = "n", xlim = c(0, 1), ylim = c(0, 1))
for (i in 2:60) points(u[i - 1], u[i], pch = 16)
points(u[60], u[1], pch = 16)
#
# "Random numbers fall mainly in the plane" --- Marsaglia
#
u = runif(60)
plot(0, 0, type = "n", xlim = c(0, 1), ylim = c(0, 1))
for (i in 2:60) points(u[i - 1], u[i], pch = 16)
points(u[60], u[1], pch = 16)
#
# Subtopic 2: Bernoulli trials
#
x = runif(14) < 0.3
x
x = as.numeric(x)
x
sum(x)
rbinom(1, 14, 0.3)
rbinom(17, 14, 0.3)
#
# Subtopic 3: Monte Carlo simulation
#
#
# Example 1: Three men and two women sit in a row of chairs in
# a random order. Use Monte Carlo simulation to estimate the
# probability that men and women alternate. Women even, men odd.
# Analytic solution: 1 / 10
#
nrep = 100000
count = 0
for (i in 1:nrep) {
x = sample(c("M", "W", "M", "W", "M"))
if (all(x[c(1, 3, 5)] == "M")) count = count + 1
}
print(count / nrep)
nrep = 100000
count = 0
for (i in 1:nrep) {
x = sample(5)
if (x[1] * x[3] * x[5] == 15) count = count + 1
}
print(count / nrep)
#
# Example 2:
# Find the probability of rolling a total of 10 with three rolls of
# a fair die.
#
nrep = 100000
count = 0
for (i in 1:nrep) {
x = sample(1:6, 3, replace = TRUE)
if (sum(x) == 10) count = count + 1
}
print(count / nrep)
#
# Example 3:
# Urn 1 contains three white balls and two black balls.
# Urn 2 contains one white ball and three black balls.
# A ball is selected at random from Urn 1 and transferred
# to Urn 2. Next, a ball is selected at random from Urn 2
# and transferred to Urn 1. Finally, a ball is selected
# at random from Urn 1. Find the probability that the
# ball selected on this third draw is black.
#
nrep = 100000
count = 0
for (i in 1:nrep) {
urn1 = c(0, 0, 0, 1, 1)
urn2 = c(0, 1, 1, 1)
x = sample(5, 1)
urn2 = c(urn2, urn1[x])
y = sample(5, 1)
urn1[x] = urn2[y]
if (sample(urn1, 1) == 1) count = count + 1
}
print(count / nrep)
#
# Example 4:
# What is the 90th percentile of the distance between two points
# chosen at random in the interior of a unit square?
#
nrep = 1000000
x1 = runif(nrep)
x2 = runif(nrep)
y1 = runif(nrep)
y2 = runif(nrep)
d = sqrt((x1 - x2) ^ 2 + (y1 - y2) ^ 2)
sort(d)[0.9 * nrep]
#
# Example 5:
# Xavier and Yolanda agree to meet at the sun dial at 2:00 PM.
# Neither, however, is particularly punctual, and their actual arrival
# times to the sun dial are independent and uniformly distributed between
# 2:00 PM and 3:00 PM.
# Assuming that each will wait~15 minutes for the other, find the
# probability that they will actually meet.
#
nrep = 1000000
x = runif(nrep, 0, 60)
y = runif(nrep, 0, 60)
sum(abs(x - y) < 15) / nrep
#
# Example 6:
# Let the continuous random variables A, B, and C
# be mutually independent U(0, 1) random variables.
# Find the probability that the random quadratic equation
# A x^2 + B x + C = 0 has real roots.
# Analytic: 0.2544
#
nrep = 1000000
A = runif(nrep)
B = runif(nrep)
C = runif(nrep)
mean(B ^ 2 - 4 * A * C > 0)
#
# Example 7:
# Russell on a 7 x 16 grid
# six north moves, 15 east moves
# mailbox at 8 north and 2 east moves
#
nrep = 100000
count = 0
moves = c(rep(1, 6), rep(0, 15))
for (i in 1:nrep) {
x = sample(moves)
if (cumsum(x)[10] == 2) count = count + 1
}
print(count / nrep)
# Analytic: choose(10, 8) * choose(11, 7) / choose(21, 15)
# 0.2737
#
# Example 2 revisited with replicate
#
roll = function() sum(sample(6, 3, replace = TRUE))
y = replicate(100000, roll())
mean(y == 10)
#
# one line version
#
mean(replicate(100000, sum(sample(6, 3, replace = TRUE)) == 10))
#
# Example 7 revisited with replicate
#
moves = c(rep(1, 6), rep(0, 15))
mean(replicate(100000, cumsum(sample(moves))[10] == 2))