Statistical Computing and Simulation

library(ggplot2)
library(kableExtra)
library(dplyr)
library(magrittr)
library(kableExtra)
library(knitr)
options(scipen = 10)

Question 01

\(I = \int_{-20}^{20} e^{-|x|}\,dx.\)

Find I by calculus techniques.

set.seed(3)
# Find I by calculus techniques.
I <- function(x){
  return(exp(-abs(x)))
}
Ans1 <- integrate(I, -20, 20)
Ans1

## 2 with absolute error < 0.0000003

Find I by simultaion

N <- 5000 # 5000 samples of size 1000.
n <- 1000

# (1) Monte Carlo

theta.hat <- NULL
for( i in 1:N){
x <- runif(n, -20, 20)
theta.hat[i] <- mean(40*I(x))
}
plot(theta.hat, main = "Monte Carlo",ylab = "I_hat")

Mc.mean <- mean(theta.hat)
Mc.var <- var(theta.hat)

# (2) Importance Sampling Monte Carlo
# N(0,1)

theta.hat <- NULL

for( i in 1:N){
x <- rnorm(n,0,1)
theta.hat[i] <- mean(40*I(x)*dunif(x,-20,20)/dnorm(x,0,1))

}
plot(theta.hat, main = "Importance Sampling Monte Carlo\n N(0,1)",ylab = "I_hat")

t1 <- theta.hat[-which(theta.hat>3)]
plot(t1, main = "Importance Sampling Monte Carlo\n N(0,1)",ylab = "I_hat")

N1.mean <- mean(theta.hat)
N1.var <- var(theta.hat)

# (3) Importance Sampling Monte Carlo
# N(0,5)

theta.hat <- NULL
for( i in 1:N){
  x <- rnorm(n,0,20)
  theta.hat[i] <- mean(40*I(x)*dunif(x,-20,20)/dnorm(x,0,20))
}
plot(theta.hat, main = "Importance Sampling Monte Carlo\n N(0,5)",ylab = "I_hat")

N2.mean <- mean(theta.hat)
N2.var <- var(theta.hat)

# Compare three methods
Q1M <- cbind(Mc.mean, N1.mean, N2.mean) %>% round(., digits = 5)
Q1V <- cbind(Mc.var, N1.var, N2.var) %>% round(., digits = 5)
Q1 <- rbind(Q1M,Q1V)
rownames(Q1) <- c("Mean", "Variance")
colnames(Q1) <- c("MC", "N(0,1)", "N(0,5)")
kable(Q1) %>% kable_styling()

## Warning in kable_styling(.): Please specify format in kable. kableExtra can
## customize either HTML or LaTeX outputs. See https://haozhu233.github.io/
## kableExtra/ for details.

	MC	N(0,1)	N(0,5)
Mean	1.99851	1.98926	1.99481
Variance	0.03619	0.14084	0.04522

Question 02

# Find N(0,1) mean and variance by important sampling with t3

set.seed(3)
ME <- function(x){
  1/sqrt((2*pi))*exp(-x^2/2)
}
N = 5000
# Mean
# diffn
diffn_mean <- NULL
diffn_var <- NULL
n <- seq(100,1000,by = 100)
for(j in 1:length(n)){
theta.hat <- NULL
for(i in 1:N){
x <- rt(n[j],3)
theta.hat[i] <- mean(x*ME(x)*dunif(x,-30,30)/dt(x,3)*60)
}
diffn_mean <- c(diffn_mean,mean(theta.hat))
diffn_var <- c(diffn_var,var(theta.hat))
}

diffn <- data.frame(c(1:10),diffn_mean, diffn_var) %>% round(., digits = 5)
t(diffn)

##               [,1]     [,2]    [,3]     [,4]     [,5]     [,6]     [,7]
## c.1.10.    1.00000  2.00000 3.00000  4.00000  5.00000  6.00000  7.00000
## diffn_mean 0.00010 -0.00050 0.00046 -0.00059 -0.00009 -0.00043 -0.00029
## diffn_var  0.00946  0.00463 0.00318  0.00231  0.00189  0.00156  0.00133
##                [,8]    [,9]    [,10]
## c.1.10.     8.00000 9.00000 10.00000
## diffn_mean -0.00059 0.00054 -0.00023
## diffn_var   0.00118 0.00104  0.00093

# 繪圖
# mean
ggplot(diffn) + 
  geom_point(aes(x = diffn$c.1.10., y = diffn_mean)) +
  geom_smooth(aes(x = diffn$c.1.10., y = diffn_mean), se = FALSE) +
  labs(title="N(0,1) mean of important sampling with t(3)\n Mean",
       x="Index",
       y="Mean of estimator of mean") +
  theme_bw() 

## `geom_smooth()` using method = 'loess'

# Variance
ggplot(diffn) + 
  geom_point(aes(x = diffn$c.1.10., y = diffn_var)) +
  geom_smooth(aes(x = diffn$c.1.10., y = diffn_var),se = FALSE) +
  labs(title="N(0,1) mean of important sampling with t(3)\n Variance",
       x="Index",
       y="Variance of estimator of mean") +
  theme_bw() 

## `geom_smooth()` using method = 'loess'

# Variance
diffn_mean <- NULL
diffn_var <- NULL
n <- seq(100,1000,by = 100)
for(j in 1:length(n)){
theta.hat <- NULL
for(i in 1:N){
  x <- rt(n[j],3)
  theta.hat1 <- mean(x^2*ME(x)*dunif(x,-30,30)/dt(x,3)*60)
  theta.hat2 <- mean(x*ME(x)*dunif(x,-30,30)/dt(x,3)*60)
  theta.hat[i] <- theta.hat1 - (theta.hat2)^2
}
diffn_mean <- c(diffn_mean,mean(theta.hat))
diffn_var <- c(diffn_var,var(theta.hat))
}

diffn <- data.frame(c(1:10),diffn_mean, diffn_var) %>% round(., digits = 5)
t(diffn)

##               [,1]    [,2]    [,3]    [,4]    [,5]    [,6]    [,7]    [,8]
## c.1.10.    1.00000 2.00000 3.00000 4.00000 5.00000 6.00000 7.00000 8.00000
## diffn_mean 0.99049 0.99495 0.99774 0.99830 0.99883 0.99896 0.99943 0.99884
## diffn_var  0.01102 0.00540 0.00374 0.00266 0.00222 0.00183 0.00159 0.00142
##               [,9]    [,10]
## c.1.10.    9.00000 10.00000
## diffn_mean 0.99857  0.99917
## diffn_var  0.00126  0.00109

# 繪圖
# mean
ggplot(diffn) + 
  geom_point(aes(x = diffn$c.1.10., y = diffn_mean)) +
  geom_smooth(aes(x = diffn$c.1.10., y = diffn_mean), se = FALSE) +
  labs(title="N(0,1) variance of important sampling with t(3)\n Mean",
       x="Index",
       y="Mean of estimator of variance") +
  theme_bw() 

## `geom_smooth()` using method = 'loess'

# Variance
ggplot(diffn) + 
  geom_point(aes(x = diffn$c.1.10., y = diffn_var)) +
  geom_smooth(aes(x = diffn$c.1.10., y = diffn_var), se = FALSE) +
  labs(title="N(0,1) variance of important sampling with t(3)\n Variance",
       x="Index",
       y="Variance of estimator of variance") +
  theme_bw() 

## `geom_smooth()` using method = 'loess'

Find t3 mean and variance by important sampling with N(0,1)

set.seed(3)
v <- 3
tv <- function(x,v){
  z <- gamma((v+1)/2)/sqrt(v*pi)/gamma(v/2)/(1+(x^2)/v)^((v+1)/2)
  return(z)
}

# mean
diffn_mean <- NULL
diffn_var <- NULL
n <- seq(100,1000,by = 100)
for(j in 1:length(n)){
theta.hat <- NULL
for(i in 1:N){
x <- rnorm(n[j],0,1)
theta.hat[i] <- mean(x*tv(x,v)/dnorm(x,0,1))
}
diffn_mean <- c(diffn_mean,mean(theta.hat))
diffn_var <- c(diffn_var,var(theta.hat))
}
diffn <- data.frame(c(1:10),diffn_mean, diffn_var) %>% round(., digits = 5)
t(diffn)

##                [,1]     [,2]    [,3]    [,4]    [,5]     [,6]    [,7]
## c.1.10.     1.00000  2.00000 3.00000 4.00000 5.00000  6.00000 7.00000
## diffn_mean -0.03472 -0.07704 0.00640 0.01828 0.02777 -0.00569 0.00535
## diffn_var   4.17662 37.85109 0.14083 0.33347 7.30422  0.54173 0.35007
##               [,8]     [,9]    [,10]
## c.1.10.    8.00000  9.00000 10.00000
## diffn_mean 0.02893 -0.09281 -0.00659
## diffn_var  4.66431 21.42051  0.64627

# 繪圖

# mean
ggplot(diffn) + 
  geom_point(aes(x = diffn$c.1.10., y = diffn_mean)) +
  geom_smooth(aes(x = diffn$c.1.10., y = diffn_mean), se = FALSE) +
  labs(title="t(3) mean of important sampling with N(0,1)\n Mean",
       x="Index",
       y="mean of estimator of mean") +
  theme_bw() 

## `geom_smooth()` using method = 'loess'

# Variance
ggplot(diffn) + 
  geom_point(aes(x = diffn$c.1.10., y = diffn_var)) +
  geom_smooth(aes(x = diffn$c.1.10., y = diffn_var), se = FALSE) +
  labs(title="t(3) mean of important sampling with N(0,1)\n Variance",
       x="Index",
       y="Variance of estimator of mean") +
  theme_bw() 

## `geom_smooth()` using method = 'loess'

# Var
diffn_mean <- NULL
diffn_var <- NULL
n <- seq(100,1000,by = 100)
for(j in 1:length(n)){
theta.hat <- NULL
for(i in 1:N){
  x <- rnorm(n[j],0,1)
  theta.hat1 <- (mean(x*tv(x,v)/dnorm(x,0,1)))^2
  theta.hat2 <- mean((x^2)*tv(x,v)/dnorm(x,0,1))
  theta.hat[i] <- theta.hat2 - theta.hat1
}
diffn_mean <- c(diffn_mean,mean(theta.hat))
diffn_var <- c(diffn_var,var(theta.hat))
}

diffn <- data.frame(c(1:10),diffn_mean, diffn_var) %>% round(., digits = 5)
t(diffn)

##                    [,1]       [,2]       [,3]      [,4]     [,5]      [,6]
## c.1.10.         1.00000    2.00000    3.00000   4.00000  5.00000   6.00000
## diffn_mean    -12.69518    0.66387    0.50941   1.10809  1.42652   1.33154
## diffn_var  663024.83381 1238.09275 4433.13618 387.92210 16.91081 272.91020
##                    [,7]        [,8]     [,9]      [,10]
## c.1.10.         7.00000     8.00000  9.00000   10.00000
## diffn_mean     -8.73801    -0.81247  1.53318    0.52701
## diffn_var  518980.65200 22146.27862 19.25389 5697.17336

# 繪圖

# mean
ggplot(diffn) + 
  geom_point(aes(x = diffn$c.1.10., y = diffn_mean)) +
  geom_smooth(aes(x = diffn$c.1.10., y = diffn_mean), se = FALSE) +
  labs(title="t(3) Variance of important sampling with N(0,1)\n Mean",
       x="Index",
       y="mean of estimator of variance") +
  theme_bw() 

## `geom_smooth()` using method = 'loess'

# Variance
ggplot(diffn) + 
  geom_point(aes(x = diffn$c.1.10., y = diffn_var)) +
  geom_smooth(aes(x = diffn$c.1.10., y = diffn_var), se = FALSE) +
  labs(title="t(3) Variance of important sampling with N(0,1)\n Variance",
       x="Index",
       y="Variance of estimator of variance") +
  theme_bw() 

## `geom_smooth()` using method = 'loess'