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 Transformation Method

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 :

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)
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)

2.2 Monte-Carlo Method

2.3 hit-or-miss approach

Website (suggested) version :

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