Website (suggested) version :

https://rpubs.com/Isaiah-Mireles/1447225

1 Q1

1.1 Choleski factorization method

Sigma <- matrix(c(
  2.8, 0,   0.2,  2,
  0,   1.7, 2,    0,
  0.2, 2,   3.6, -1.2,
  2,   0,  -1.2,  3
), nrow = 4, byrow = TRUE)

Sigma

##      [,1] [,2] [,3] [,4]
## [1,]  2.8  0.0  0.2  2.0
## [2,]  0.0  1.7  2.0  0.0
## [3,]  0.2  2.0  3.6 -1.2
## [4,]  2.0  0.0 -1.2  3.0

A <- Sigma |> chol()
t(A)%*%A # worked

##      [,1] [,2] [,3] [,4]
## [1,]  2.8  0.0  0.2  2.0
## [2,]  0.0  1.7  2.0  0.0
## [3,]  0.2  2.0  3.6 -1.2
## [4,]  2.0  0.0 -1.2  3.0

Sigma

##      [,1] [,2] [,3] [,4]
## [1,]  2.8  0.0  0.2  2.0
## [2,]  0.0  1.7  2.0  0.0
## [3,]  0.2  2.0  3.6 -1.2
## [4,]  2.0  0.0 -1.2  3.0

1.1.1 Simulate Rand. Normal w/ unif(0,1)

set.seed(123)
n <- 1000
U <- runif(n)
V <- runif(n)
Z1 <- sqrt(-2*log(U))*cos(2*pi*V)
Z2 <- sqrt(-2*log(U))*sin(2*pi*V)
plot(Z1, Z2)

clearly it worked – characteristic football shape

U2 <- runif(n)
V2 <- runif(n)

Z3 <- sqrt(-2 * log(U2)) * cos(2 * pi * V2)
Z4 <- sqrt(-2 * log(U2)) * sin(2 * pi * V2)

Z <- cbind(Z1, Z2, Z3, Z4)

So now we just need to move it and stretch it – translation & dilation

dim(A)

## [1] 4 4

dim(Z)

## [1] 1000    4

# dilate 
X <- Z %*% A 

var(X[,1])

## [1] 2.895608

# notice the variance is close, so good

# translate
mu <- c(2, 1.5, 3, 1) 
X <- X + mu
X[,1] |> hist() # centered @ 2

1.2 pairs()

pairs(X)

notice they are all shaped like footballs, centered nicely and we can see varying as much as we should expect – stretching our eigen vectors to fit the variance of each random variable, orienting at an angle as expected with correlated variables

2 Q2

2.1 (b) hit-or-miss approach

2.1.1 Composition

for this course, we will assume we are stripped of our typical stats packages and must build from the ground up–

and so, remember :

set.seed(123)

x <- rgamma(1000, shape = 3, rate = 2)

hist(x, probability = TRUE)
curve(dgamma(x, shape = 3, rate = 2),
      add = TRUE, col = "red", lwd = 2)

cheating^

so, recall : convolution :

and,

# exponentials / inv. cdf method : 
set.seed(123)
ExpGen <- function(n, alpha, beta) {
  E <- matrix(0, nrow = n, ncol = alpha)

  for (i in 1:alpha) {
    U <- runif(n)
    E[, i] <- -log(1 - U) / beta
  }

  return(E)
}
ExpGen(1000, 3, 2)[,2] |> hist()

yeh – exponential

# convolution
E <- ExpGen(1000, 3, 2)
G <- rowSums(E) 
hist(G, probability=T, main = "Gamma")
curve(dgamma(x, shape = 3, rate = 2),
      add = TRUE, col = "red", lwd = 2)

yep, its a fit – so now we just need to make beta (B) dist from gamma (G)

set.seed(69)
G2 <- ExpGen(1000, 2, 2) |> rowSums()
B <- G / (G + G2)
hist(B, probability = TRUE,
     main = "Beta(3,2) via Gamma",
     xlab = "B")

curve(dbeta(x, shape1 = 3, shape2 = 2),
      from = 0, to = 1,
      add = TRUE,
      col = "red",
      lwd = 2)

now that i simulated the dist. I must just count whats below.

x <- seq(0.1, 0.9, by = 0.1)
hm.estimates <- numeric(length(x)) #bucket


for(i in 1: length(x)){
  hm.estimates[i] <- mean(B <= x[i])
}

hm.estimates

## [1] 0.001 0.022 0.070 0.191 0.301 0.477 0.648 0.817 0.954

2.2 (a) Monte-Carlo Method

here we derive the equation to simulate below :

notice at first i forgot to mult. by x but corrected that.

set.seed(102)
# numb. samp.
m <- 100000
# input values : 
x <- seq(0.1, 0.9, by = 0.1)

u <- runif(m) #uniform
mc <- numeric(length(x)) #storage

for (i in 1:length(x)) {

  g <- 12 * x[i]^3 * u^2 * (1 - x[i] * u)

  mc[i] <- mean(g)

}

mc

## [1] 0.003708779 0.027261161 0.083875794 0.179544104 0.313030295 0.475871351
## [7] 0.652377031 0.819629874 0.947485195

2.3 (c) compare

set.seed(1243)
# exact values
exact <- pbeta(x, shape1 = 3, shape2 = 2)

results <- data.frame(
  x = x,
  Hit.or.Miss = hm.estimates,
  Monte.Carlo = mc,
  Exact = exact,
  HM.Error = hm.estimates - exact,
  MC.Error = mc - exact
)

results

##     x Hit.or.Miss Monte.Carlo  Exact HM.Error      MC.Error
## 1 0.1       0.001 0.003708779 0.0037  -0.0027  8.779383e-06
## 2 0.2       0.022 0.027261161 0.0272  -0.0052  6.116103e-05
## 3 0.3       0.070 0.083875794 0.0837  -0.0137  1.757936e-04
## 4 0.4       0.191 0.179544104 0.1792   0.0118  3.441036e-04
## 5 0.5       0.301 0.313030295 0.3125  -0.0115  5.302954e-04
## 6 0.6       0.477 0.475871351 0.4752   0.0018  6.713513e-04
## 7 0.7       0.648 0.652377031 0.6517  -0.0037  6.770315e-04
## 8 0.8       0.817 0.819629874 0.8192  -0.0022  4.298740e-04
## 9 0.9       0.954 0.947485195 0.9477   0.0063 -2.148053e-04

plot(x, exact,
     type = "l", lwd = 2, col = "black",
     ylim = c(0, 1),
     ylab = "CDF", xlab = "x")

lines(x, hm.estimates,
      type = "b", pch = 16, col = "red")

lines(x, mc,
      type = "b", pch = 17, col = "blue")

legend("bottomright",
       legend = c("Exact", "Hit-or-Miss", "Monte Carlo"),
       col = c("black", "red", "blue"),
       lty = c(1, 1, 1),
       pch = c(NA, 16, 17),
       cex = 0.4)

what we see is that none of the methods did terribly bad but its clear that our monte carlo method is doing better. Our hit or miss method appears to be missing at the top and bottom values