# 1. Introduction to Monte Carlo Methods Question 1
n = 1e+06
x = runif(n, -5, 3)
y = runif(n, -5, 3)
mean((2 * x^2 + x * y) < 4)

## [1] 0.3856

# Answer is 0.385523

# Question 2
u = runif(1e+06, -2, 2)
4 * mean(exp(u + (u^2)))

## [1] 93.48

# Answer is 93.41703

# 2. Simulating Samples from Arbitrary Probability Distributions Question 1
x = c(-3, 0, 2)
p = c(0.45, 0.3, 0.25)
cdf = cumsum(p)
plot(x, p, "h")

plot of chunk unnamed-chunk-1


u = runif(1e+06)
sample = x[1] * (u <= cdf[1]) + x[2] * (u > cdf[1] & u <= cdf[2]) + x[3] * (u > 
    cdf[2] & u <= cdf[3])
hist(sample, 3, freq = F, col = "blue")

plot of chunk unnamed-chunk-1


# I have to do the following on paper
mean(sample < 2)

## [1] 0.7508

# Answer is 0.750604
mean(sample <= 2)

## [1] 1

# Answer is 1
mean(sample)

## [1] -0.8555

# Answer is -0.852132
mean(((sample)^2) - (2 * sample))

## [1] 6.77

# Answer is 6.75462

# Question 2
u = runif(1e+06)
sample = 5 * u^(1/3)
hist(sample, freq = F, 100, col = "blue")

f = function(x) {
    3/125 * x^2
}
curve(f, 0, 5, add = T, lwd = 4, col = "red")

plot of chunk unnamed-chunk-1


mean(sample)

## [1] 3.752

# Answer is 3.748957
mean(exp((-sample)^2 - sample))

## [1] 31649180

# Answer is 31565967
mean(sample > 2)

## [1] 0.9363

# Answer is 0.935949
mean(1 < sample & sample < 4)

## [1] 0.5035

# Answer is 0.504287

# Question 3
u = runif(1e+06)
sample = ((log(1 - u))^2)/9
hist(sample, 100, freq = F, col = "blue")

f = function(x) {
    3 * exp(-3 * sqrt(x))/(2 * sqrt(x))
}
curve(f, 0, 1, add = T, lwd = 4, col = "red")

plot of chunk unnamed-chunk-1


# Question 4
u = runif(1e+06)
sample = -1/2 + 1/2 * (sqrt(1 + 8 * u))
hist(sample, 100, freq = F, col = "blue")

f = function(x) {
    1/2 + x
}
curve(f, 0, 1, add = T, lwd = 4, col = "red")

plot of chunk unnamed-chunk-1


# Question 6
f = function(x) {
    (x/24 + 2/24)
}
plot(f, 2, 6)

plot of chunk unnamed-chunk-1


u = runif(1e+06)
sample = 4 * sqrt(1 + 3 * u) - 2
hist(sample, 100, freq = F, col = "blue")

plot of chunk unnamed-chunk-1

mean(sample)

## [1] 4.223

# Answer is 4.221868

f <- function(x) {
    (((x - 2)/12) - ((x - 6)/24))
}
x = runif(1e+06, 2, 6)
y = runif(1e+06, 0, f(0))
test = (y <= f(x))
sample = x[test]
hist(sample, 100, freq = F, col = "blue")

plot of chunk unnamed-chunk-1

mean(sample)

## [1] 4

# Answer is 4.000738
var(sample)

## [1] 1.333

# Answer is 1.330839