A3_Quadrature_MC.Rmd

---
title: "STAT 812 - Assignment 3"
author: "Kangjie Zhang"
date: "November 29, 2018"
output:
  md_document:
    variant: markdown_github
    df_print: paged
    toc: yes 
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## Method (1) 
```{r}
## midpoint rule
## a function computing the sum of numbers represented with logarithm
## lx     --- a vector of numbers, which are the log of another vector x.
## the log of sum of x is returned
log_sum_exp <- function(lx)
{   mlx <- max(lx)
mlx + log(sum(exp(lx-mlx)))
}

## the generic function for approximating 1-D integral with midpoint rule
## the logarithms of the function values are passed in
## the log of the integral result is returned
## log_f  --- a function computing the logarithm of the integrant function
## range  --- the range of integral varaible, a vector of two elements
## n      --- the number of points at which the integrant is evaluated
## ...    --- other parameters needed by log_f
log_int_mid <- function(log_f, range, n,...)
{   if(range[1] >= range[2]) 
  stop("Wrong ranges")
  h <- (range[2]-range[1]) / n
  v_log_f <- sapply(range[1] + (1:n - 0.5) * h, log_f,...)
  log_sum_exp(v_log_f) + log(h)       
}

#Find the true value of E[a(x)] by midpoint rule
a_tnorm_mid <- function(fa,l,u,n)
{
  log_f <- function(x)  log(fa(x))+dnorm(x,log = TRUE)
  log_integral<-log_int_mid(log_f=log_f, range=c(l,u), n=n)
  exp(log_integral)/(pnorm(u)-pnorm(l))
}
```

## Method (2) 
```{r}
library(truncnorm)
#1 Draw samples from a truncated normal distribution
a_tnorm_naive1 <- function(fa,l,u,n)
{
  x<-rtruncnorm(n, a=l, b=u, mean = 0, sd = 1)
  mean(fa(x))  
}

#2 Draw samples from N(0,1), and keep samples within interval(l,u)
a_tnorm_naive <- function(fa,l,u,n)
{
  x <- rep (0, n)
  for (i in 1:n)
  {
    rej <- TRUE
    while (rej)
    {
      x[i] <- rnorm (1)
      if (x[i] >= l & x[i] <= u) rej <- FALSE
    }
  }
  mean(fa(x))  
}
```

## Method (3) 
### $l(x)$ and $\tilde{f}(x)$
$$(l,\ \log(l)), \ (u, \ \log(u)), \ \ \ \ \ \ \ let\ \ m=\frac{l+u}{2} $$

$$l(x)=-\frac{1}{2}(l+u)x+\frac{ul}{2}-\frac{log(2\pi)}{2}=-mx+\frac{ul}{2}-\frac{log(2\pi)}{2}$$ 

$$\tilde{f}(x)=\exp(-mx+\frac{ul}{2}-\frac{log(2\pi)}{2}), \ \ \ x \in(l,u)$$

### PDF of g(x)
$$g(x)=a\tilde{f}(x), \ a=\frac{1}{\int_{l}^{u}\tilde{f}(x)}=\frac{m}{c(e^{-ml}-e^{-mu})}, \ \ c=\exp(\frac{ul}{2}-\frac{log(2\pi)}{2})$$

### CDF of g(x)
$$G(x)=\frac{e^{-ml}-e^{-mx}}{e^{-ml}-e^{-mu}},\ \ \ x \in(l,u)$$

### The inverse of G(x) 
$$m\neq0, \ \ \ G^{-1}(v) = -\frac{1}{m}\log(e^{-ml}(1-v)+e^{-mu}v), \ \ \ \ \ \ x \in(l,u)$$
$g(x)$ follows a uniform distribution unif(l,u) when $m=0$
```{r}
## Importance sampling
a_tnorm_is <- function(fa, l, u, n)
{
  x <- rep(0,n) 
  m <- (l+u)/2
  v <- runif(n)
  if (m!=0){
    x <- -1/m*log(exp(-m*l)*(1-v)+v*exp(-m*u)) #m cannot be zero
  } else {
    x <- runif(n,l,u)#when m=0, l=-u, gx~unif(l,u)
  }
  fx <- dnorm(x)/(pnorm(u)-pnorm(l))
  gx <- exp(-m*x+0.5*u*l-0.5*log(2*pi))
  W <- fx/gx
  ahat <- sum (fa(x) * W) / sum (W)
  attr(ahat, "effective sample size") <- 1/sum((W/sum(W))^2)
  ahat
}
```

## Method (4) 
### Tangent line of $log(\phi(x))$ at $\frac{l+u}{2}$
let $m=\frac{l+u}{2}$,
$$\log \phi(x)=-\frac{1}{2}\log(2\pi)-\frac{x^2}{2}$$
$$[\log \phi(x)]'=-x$$
$$\log l(x)=-m(x-m)-\frac{1}{2}\log(2\pi)-\frac{m^2}{2}=-mx+\frac{m^2}{2}-\frac{1}{2}\log(2\pi)$$

### PDF of g(x)
$$l(x)=\exp(-mx+\frac{m^2}{2}-\frac{1}{2}\log(2\pi))$$
$$g(x)=$$
$$g(x)=al(x), \ a=\frac{1}{\int_{l}^{u}l(x)}=\frac{m}{c(e^{-ml}-e^{-mu})}, \ \ c=\exp(\frac{m^2}{2}-\frac{log(2\pi)}{2})$$

### CDF of g(x)
$$G(x)=\frac{e^{-ml}-e^{-mx}}{e^{-ml}-e^{-mu}},\ \ \ x \in(l,u)$$

### The inverse of G(x) 
$$m\neq0, \ \ \ G^{-1}(v) = -\frac{1}{m}\log(e^{-ml}(1-v)+e^{-mu}v), \ \ \ \ \ \ x \in(l,u)$$
$g(x)$ follows a uniform distribution unif(l,u) when $m=0$

```{r}
## log(g(x))
log_g <- function(x, l, u)
{
  m <- (l+u)/2
  -m*x+0.5*m^2-0.5*log(2*pi)
}

## Rejection sampling
a_tnorm_rej <- function(fa, l, u, n)
{
  sample<-rep(0,n)
  for (i in 1:n) {
    rejected <- TRUE
    while (rejected) {
      v <- runif(1)
      m <- (l+u)/2
      if (m!=0){
        sample[i] <- -1/m*log(exp(-m*l)*(1-v)+v*exp(-m*u)) #m cannot be zero
      } else {
        sample[i] <- runif(1,l,u)#when m=0, l=-u, gx is a uniform dist
      }
      U <- runif(1)
      rejected <- (log(U) > dnorm(sample[i],log = TRUE) - log_g(sample[i],l,u) )
    }
  }
  mean(fa(sample))
}
```



# Question (a)
```{r}
## a(x)=x^2
fa <- function(x) x^2
## true value of E[a(X)]
A1 <- a_tnorm_mid (fa,1,2,1000000)
A2 <- a_tnorm_mid (fa,-1,1,1000000)
A3 <- a_tnorm_mid (fa,-1,0,1000000)
```


# Question (b)
### Test 1 (Use A1, l=1, u=2)
```{r}
system.time(
a_tnorm_naive(fa,1,2,1000) ## naive
)

system.time(
a_tnorm_is(fa, 1, 2, 1000) ## importance sampling
)

system.time(
a_tnorm_rej(fa, l=1, u=2, 1000) ## rejection sampling
)

## replicate 5000 times to find MSE
times.naive1 <- system.time(
naive1 <- replicate (5000, a_tnorm_naive(fa,1,2,1000)) 
)

times.is1 <- system.time(
is1 <- replicate (5000, a_tnorm_is(fa, 1, 2, 1000)) 
)

times.rej1 <- system.time(
rej1 <- replicate (5000, a_tnorm_rej(fa, 1, 2, 1000)) 
)

MSE.naive1 <- mean ((naive1-A1)^2)
MSE.is1 <- mean ((is1-A1)^2)
MSE.rej1 <- mean ((rej1-A1)^2)

MSE.naive1; times.naive1
MSE.is1; times.is1
MSE.rej1; times.rej1

ylim <- range (naive1, is1, rej1)
plot (naive1, ylim = ylim, main = "Naive MC") 
abline (h = A1)
plot (is1, ylim = ylim, main = "Importance Sampling") 
abline (h = A1)
plot (rej1, ylim = ylim, main = "Rejection Sampling") 
abline (h = A1)
```

### Test 2 (Use A2, l=-1, u=1)
```{r}
## replicate 5000 times to find MSE
times.naive2 <- system.time(
naive2 <- replicate (5000, a_tnorm_naive(fa,-1,1,1000)) 
)

times.is2 <- system.time(
is2 <- replicate (5000, a_tnorm_is(fa, -1, 1, 1000)) 
)

times.rej2 <- system.time(
rej2 <- replicate (5000, a_tnorm_rej(fa, -1, 1, 1000)) 
)

MSE.naive2 <- mean ((naive2-A2)^2)
MSE.is2 <- mean ((is2-A2)^2)
MSE.rej2 <- mean ((rej2-A2)^2)

MSE.naive2; times.naive2
MSE.is2; times.is2
MSE.rej2; times.rej2

ylim <- range (naive2, is2, rej2)
plot (naive2, ylim = ylim, main = "Naive MC") 
abline (h = A2)
plot (is2, ylim = ylim, main = "Importance Sampling") 
abline (h = A2)
plot (rej2, ylim = ylim, main = "Rejection Sampling") 
abline (h = A2)
```

### Test 3 (Use A3, l=-1, u=0)
```{r}
## replicate 5000 times to find MSE
times.naive3 <- system.time(
naive3 <- replicate (5000, a_tnorm_naive(fa,-1,0,1000)) 
)

times.is3 <- system.time(
is3 <- replicate (5000, a_tnorm_is(fa, -1, 0, 1000)) 
)

times.rej3 <- system.time(
rej3 <- replicate (5000, a_tnorm_rej(fa, -1, 0, 1000)) 
)

MSE.naive3 <- mean ((naive3-A3)^2)
MSE.is3 <- mean ((is3-A3)^2)
MSE.rej3 <- mean ((rej3-A3)^2)

MSE.naive3; times.naive3
MSE.is3; times.is3
MSE.rej3; times.rej3

ylim <- range (naive3, is3, rej3)
plot (naive3, ylim = ylim, main = "Naive MC") 
abline (h = A3)
plot (is3, ylim = ylim, main = "Importance Sampling") 
abline (h = A3)
plot (rej3, ylim = ylim, main = "Rejection Sampling") 
abline (h = A3)
```

### Test other intervals (l=-5, u=-3)
```{r}
a_tnorm_mid (fa,-5,-3,1000000)
a_tnorm_naive(fa,-5,-3,1000)
a_tnorm_is(fa,-5,-3,1000)
a_tnorm_rej(fa,-5,-3,1000)
```

## Comments
Above all, methods 2-4 all work well for estimating the E[a(x)], all of them got very close estimations to the true values calculated by the midpoint rule method. For three different sets of intervals from test1 to test3, Method (4) (Rejection sampling) is the most time consuming one; Method (3) (Importance sampling) is the most efficiency method among these three, which has the fastest speed. Also difference values of l and u would lead to difference speed for 3 methods.