# Question 1b

# Function draws a pair of N(0,1) random numbers using the Polar Marsaglia method
polarMarsaglia <- function() {
  repeat {                                    # Keep sampling until inside unit circle
        u <- -1 + 2*runif(2)                  # Draw U_1,U_2 ~ U[-1,1]
        r2 <- sum(u**2)
        if (r2<=1)                            # If radius is less than one, we are done
            break
  }
  factor <- sqrt(-2*log(r2)/r2)
  u * factor
}

x <- numeric(100)                             # Create vector to store result
for (i in 1:50)
  x[2*i+c(-1,0)] <- polarMarsaglia()          # Fill with pairs of values

####################################################

impsamp.var <- function(sigma2) {
  sigma2 / sqrt(2*sigma2-1) * exp(1/(2*sigma2-1)) - 1
}

optimize(impsamp.var, c(0.5,5))

####################################################

# Question 5a

rejection.example <- function(n, sigma2) {
  x <- numeric(n)                              # Create vector to store result
  sigma2 <- (sqrt(5)+3)/2                      # Optimal variance
  M <- sqrt(sigma2) * exp(1/(2*(sigma2-1)))    # Optimal M
  for (i in 1:n) {
    repeat {                                   # Keep proposing ...
      x.new  <- rnorm(1,1,sqrt(sigma2))        # Generate x.new ~ N(1,sigma2)
      accept <- dnorm(x.new,0,1) / (M*dnorm(x.new,1,sqrt(sigma2)))
                                               # Compute probability of acceptance
      if (runif(1)<accept)                     # ... until a proposed value is accepted
        break
    }
    x[i] <- x.new
  x
}

n <- 1000
x <- rejection.example(n, (sqrt(5)+3)/2)       # Generate sample
hist(x, freq=FALSE)                            # Draw histogram
t <- seq(-3, 3, length.out=1000)
lines(t, dnorm(t,0,1), lwd=2)                  # Add density of target dist'n

####################################################

# Question 5b
importance.sample <- function(n, sigma2) {
   x <- rnorm(n,1,sqrt(sigma2))                # Draw from instrumental dist'n
   w <- dnorm(x,0,1) / dnorm(x,1,sqrt(sigma2)) # Compute weights
   list(x=x, w=w)
}

n <- 1000
s = importance.sample(n, 2.2808)               # Generate weighted sample
sum(s$x*s$w)/sum(s$w)                          # Self-normalised estimate of E(X)
sum(s$x**2*s$w)/sum(s$w)                       # Self-normalised estimate of E(X^2)

####################################################

# Question 7a

sample.truncated.normal.a <- function(n, mu, sigma2, tau) {
  x <- numeric(n)                                 # Create vector to store result
  num.proposed <- 0                               # Set counter of proposed values
  for (i in 1:n) { 
    repeat {                                      # Keep proposing ...
      x.new <- rnorm(1,mu,sqrt(sigma2))           # Generate x.new ~ N(mu,sigma2)
      num.proposed <- num.proposed + 1            # Increment counter
      if (x_new>tau)                              # ... until x_new > tau (accepted)
        break
    }
    x[i] <- x.new
  }
  cat("Proportion of accepted values\n")          # Print proportion of accepted values         
  print(n/num.proposed) 
  x
}

####################################################

# Question 7b

x <- sample.truncated.normal.a(10, 0, 1, 4)

####################################################

# Question 7c

sample.truncated.normal.c <- function(n, tau) {
  M <- 1 / ( 2 * (1-pnorm(tau) )) * exp(-tau^2/2) 
                                                  # Compute M
  x <- numeric(n)                                 # Create vector to store result
  num.proposed <- 0                               # Set counter of proposed values
  for (i in 1:n) { 
    repeat {                                      # Keep proposing ...
      x.new <- tau + abs(rnorm(1))                # Draw new value from instrumental dist'n
      num.proposed <- num.proposed + 1            # Increment counter
      f <- dnorm(x.new) / (1-pnorm(tau))          # Evaluate target density}
      g <- 2 * dnorm(x.new,mean=tau)              # Evaluate instrumental density
      accept <- f / ( M * g )                     # Probability of accepting
      if (runif(1)<accept)                        # ... until one value is accepted
        break
    }
    x[i] <- x.new
  }    
  cat("Proportion of accepted values\n")          # Print proportion of accepted values         
  print(n/num.proposed) 
  x  
}

x <- sample.truncated.normal.c(10, 4)

####################################################

# Question 8a and b

is.beta <- function (n, alpha, beta) {
 x <- runif(n)                                 # Draw from instrumental dist'n
 w <- dbeta(x,alpha,beta)                      # Compute weights
 list(x=x, w=w)
}

n <- 10
N <- 100000
result.selfnorm = numeric(N)
result.simple = numeric(N)

for (i in 1:N) {
    s <- is.beta(n, 2, 3)                      # Draw weighted sample
    result.selfnorm[i] <- sum(s$x*(1-s$x)*s$w)/sum(s$w)
                                               # Compute self-normalised estimate
    result.simple[i] <- sum(s$x*(1-s$x)*s$w)/n # Compute simple estimate
}
bias.selfnorm <- mean(result.selfnorm) - 0.2   # Compute the bias
bias.simple <- mean(result.simple) - 0.2

var.selfnorm <- var(result.selfnorm)           # Compute the variance
var.simple <- var(result.simple) 

mse.selfnorm <- bias.selfnorm^2 + var.selfnorm # Compute m.s.e. 
mse.simple <- bias.simple^2 + var.simple

boxplot(as.data.frame(cbind(result.simple, result.selfnorm)),
        names=c("simple","self-normalised"))