# Lecture 2.1 ---------------------------------------------------

# A linear congruential generator
zx81.seed <- 1L
ZX81 <- function(n=1L, a=75L, c=74L, m=2L^16L) {
  x <- NULL
  for (i in 1:n) {
    zx81.seed <<- (a*zx81.seed + c) %% m
    x[i] <- zx81.seed
  }
  x/2^16
}
x <- ZX81(2^16)
x[1:10]
acf(x) # looks good but
hist(x[1:(2^16)])


# Generating samples from exponential distribution using inversion method 
rexp2 <- function(n, rate=1) {
  u <- runif(n)
  -log(u)/rate
}
x <- rexp2(1e+3, rate=.1)
hist(x)
microbenchmark::microbenchmark(rexp(1e5, 2.3), rexp2(1e5, 2.3))

# Inversion method for geometric distribution using explicit formula for the quantile function
rgeom2 <- function(n, prob) {
  u <- runif(n)
  ceiling(log(1-u)/log(1-prob))
}
x <- rgeom2(100,.3)
x
table(x)


# Lecture 2.2 ---------------------------------------------------

# Box-muller is faster on ARM64 processors than the default inversion method
# but maybe not on other platforms back in 2015 according to 
# https://stats.stackexchange.com/questions/132556/advantages-of-box-muller-over-inverse-cdf-method-for-simulating-normal-distribut
microbenchmark::microbenchmark(
  {RNGkind(normal.kind = "Box-Muller"); rnorm(1e4)},
  {RNGkind(normal.kind = "Inversion"); rnorm(1e4)} 
)


# Ratio-of-uniforms method, normal distribution
rnorm_ru <- function(n) {
  y <- NULL
  for (i in 1:n) {
    repeat {
      x1 <- runif(1)
      x2 <- runif(1,-2*exp(-.5), 2*exp(-.5))
      if (x2^2 <= -4*x1^2*log(x1))
        break()
    }
    y[i] <- x2/x1
  }
  y
}
hist(x <- rnorm_ru(1e+5), prob=TRUE, breaks=50)
mean(x)
sd(x)
curve(dnorm(x), add=TRUE)

# Ratio-of-uniforms method, gammadistribution assuming shape parameter alpha >= 1
rgamma_ru <- function(n, alpha) {
  # Store accepted samples in y
  y <- NULL 
  
  repeat {
    # Calculate number of remaining samples we need to accept
    m <- n - length(y) 
    
    # If m = 0 we are done
    if (m==0) 
      break()
    
    # Simulate m proposals from inside rectangle containing the set C
    x1 <- runif(m) 
    x2 <- runif(m,0, ((alpha+1)/exp(1))^((alpha+1)/2))
    
    # Compute logical vector indicating which proposals to accept
    accept <- x1 < sqrt((x2/x1)^(alpha-1)*exp(-x2/x1)) 
    
    # Compute the ratios and append these to y
    y <- c(y, (x2/x1)[accept])
  }
  
  # Return y
  y
}

x <- rgamma_ru(1e+4,1.5)
hist(x,prob=TRUE, breaks=20)
curve(dgamma(x, shape=1.5),add=TRUE)


# Lecture 3.2 ---------------------------------------------------

# Weighted resampling toy example
n <- 1000
f <- function(x) dnorm(x, mean=2, sd=.5)
g <- function(x) dnorm(x, mean=2.3, sd=.6)
x <- rnorm(n, mean=2.3, sd=.6)
w <- f(x)/g(x)/(sum(f(x)/g(x)))
m <- 50
y <- sample(x, size=m, prob=w, replace = TRUE)

par(mfrow=c(2,1))
curve(f(x), xname="x", -2, 6, col="green")
curve(g(x), xname="x", col="red", add=TRUE)
points(x, w/max(w)*.2, type="h", col="blue")
points(y, rep(0, m))
F <- function(x) pnorm(x, mean=2, sd=.5) 
curve(F(x), col="green", -2, 6)
neworder <- order(x)
w <- w[neworder]
x <- x[neworder]
Fx <- cumsum(w)
points(x,Fx, type="s", col="blue")

# Monte Carlo integration toy example
par(mfrow=c(1,1))
f <- dnorm
h <- function(x) x > 1 & x < 2
curve(f,-1,3, col="green",ylim=c(0,1.1))
curve(h, add=TRUE, col="green", type="s")
# Exact answer
pnorm(2)-pnorm(1)
# Monte Carlo integration is quite a bit off
n <- 1e3
x <- rnorm(n)
mean(h(x))
t.test(h(x))$conf
# Importance sampling estimate gets closer exact value
x <- rnorm(n, 1.5, .4)
g <- function(x) dnorm(x, 1.5, .4)
curve(g, add=TRUE, col="red")
mean(h(x)*f(x)/g(x))
t.test(h(x)*f(x)/g(x))$conf



# Lecture 4.2 ---------------------------------------------------
# Antithetic sampling

# Ex. 1: X ~ exp(1), Find E(sqrt(X))
# ordinary Monte Carlo
n <- 1000
h <- function(x) sqrt(x)
x <- rexp(2*n)
mean(h(x))
sd(h(x))/sqrt(2*n)

# Estimate using antithetic sampling
u <- runif(n)
h <- function(x) sqrt(x)
Q <- function(u) qexp(u)
mean(c(h(Q(u)), h(Q(1 - u)))) # estimate
sd((h(Q(u)) + h(Q(1 - u)))/2)/sqrt(n) # combine observations in iid pair averages
plot(h(Q(u)), h(Q(1 - u)))
cor(h(Q(u)), h(Q(1 - u)))

# Ex. 2:  X ~ N(0, 1), h(x)=x/(2^x-1), Find E(h(X))
#
# ordinary Monte Carlo estimate of E(h(X))
n <- 1e3
h <- function(x) x/(2^x-1)
x <- rnorm(2*n)
mean(h(x)) # estimate
sd(h(x))/sqrt(2*n) # standard error 

# Antithetic
x <- rnorm(n)
y <- ((h(x)+h(-x))/2)
mean(y)
sd(y)/sqrt(n)
cor(h(x),h(-x))
plot(h(x),h(-x))