############### Assignment 5.1 #################
f <- function(x){exp(-x)}
F <- function(x){1-exp(-x)}
F_inv <- function(U){-log(U)}
x <- seq(0,10,length.out=1000)
U <- seq(0,1,length.out=1000)
par(mfrow=c(1,3))
plot(x,f(x),type="l",main="f(x)")
plot(x,F(x),type="l",main="CDF of f(x)")
plot(U,F_inv(U),type="l",main="Inverse CDF of f(x)")
X <- runif(100000,0,1) # random sample from U[0,1]
Z1 <- F_inv(X)
par(mfrow=c(1,1))
hist(Z1, freq=FALSE, breaks=c(seq(0,10,length=30),Inf), xlim=c(0,10))
lines(x,f(x),type="l",main="Density function", col="red",lty=2)
################# Assignment 5.2 ######################
g <- function(y){1/(pi * (1 + y^2))}
G <- function(y){atan(y)/pi + 0.5}
G_inv <- function(U){tan((U-0.5)*pi)}
y <- seq(-10,10,length.out=1000)
U <- seq(0,1,length.out=100)
U <- seq(0.1,0.9,length.out=100)
par(mfrow=c(1,3))
plot(y,g(y),type="l",main="g(y)")
plot(y,G(y),type="l",main="CDF of g(y)")
plot(U,G_inv(U),type="l",main="Inverse CDF of g(y)")
# random sample from U[0,1]
Y <- runif(100000,0.1,0.9)
Z2 <- G_inv(Y)
par(mfrow=c(1,1))
hist(Z2, freq=FALSE, breaks=c(seq(-10,10,length=60),Inf), xlim=c(-10,10))
lines(y,g(y),type="l",main="Density function", col="red",lty=2)
################# Assignment 5.3 ######################
######### Assignment 5.3.a Possion(4) ########
m <- 100000
X = c()
for(i in 1:m){
u <- runif(1000)
j=1
prod = prod(u[1])
while(prod>exp(-4))
{
j=j+1
prod=prod(u[1:j])
}
X[i] = j-1
}
mean(X)
## [1] 3.9901
median(X)
## [1] 4
hist(X, freq=FALSE, breaks=c(seq(0,15,0.5),Inf), xlim=c(0,15))
######### Assignment 5.3.b Gamma(3,2) #########
m <- 100000
X = c()
for(i in 1:m){
Gamma_inv <- function(U){
r = 3
lamda = 2
U <- runif(r,0,1)
s = 0
s = s - sum(log(U))/lamda
return(s)
}
X[i] <- Gamma_inv(U)
}
mean(X)
## [1] 1.50413
median(X)
## [1] 1.341912
par(mfrow=c(1,1))
s = seq(0,10,length.out=10000)
hist(X, freq=FALSE, breaks=c(seq(0,10,0.2),Inf), xlim=c(0,10))
lines(s,dgamma(s,3,2),typ='l',main="PDF of Gamma(3,2)", col="red",lty=2)
######### Assignment 5.3.c Kaisq(k=5) #########
m <- 100000
X = c()
for(i in 1:m){
k = 5
Y <- rnorm(k,0,1)
Chisq_inv <- function(Y){sum(Y^2)}
X[i] <- Chisq_inv(Y)
}
mean(X)
## [1] 5.001617
median(X)
## [1] 4.350573
par(mfrow=c(1,1))
s = seq(0,20,length.out=10000)
hist(X, freq=FALSE, breaks=c(seq(0,20,0.5),Inf), xlim=c(0,20))
lines(s,dchisq(s,5),typ='l',main="PDF of Chisq(5)", col="red",lty=2)
######### Assignment 5.3.d Student's (m=4) #########
m = 4
Y <- rnorm(100000,0,1)
Z <- rchisq(100000,m)
Student_inv <- function(Y,Z){Y/(sqrt(Z/m))}
X <- Student_inv(Y,Z)
mean(X)
## [1] 0.003242462
median(X)
## [1] -0.0007616571
par(mfrow=c(1,1))
s <- seq(-40,40,length.out=1000)
hist(X, freq=FALSE, breaks=c(seq(-40,40,0.5),Inf), xlim=c(-40,40))
lines(s,dt(s,4),type="l",main="PDF of Student's t(4)", col="red",lty=2)
######### Assignment 5.3.e Student's F(k=3, m=5) #########
k = 3
m = 5
Z <- rchisq(100000,k)
W <- rchisq(100000,m)
StudF_inv <- function(Z,W){(Z/k)/(W/m)}
X <- StudF_inv(Z,W)
mean(X)
## [1] 1.673037
median(X)
## [1] 0.8991664
par(mfrow=c(1,1))
s <- seq(0,15,length.out=1000)
hist(X, freq=FALSE, breaks=c(seq(0,15,0.15),Inf), xlim=c(0,15))
lines(s,df(s,3,5),type="l",main="PDF of Student's F(3,5)", col="red",lty=2)
######### Assignment 5.3.f Beta(5,3) #########
a = 5
b = 3
Y <- rgamma(100000,a,1)
Z <- rgamma(100000,b,1)
Beta_inv <- function(Y,Z){Y/(Y+Z)}
X <- Beta_inv(Y,Z)
mean(X)
## [1] 0.6254832
median(X)
## [1] 0.6365635
par(mfrow=c(1,1))
s <- seq(0,1,length.out=1000)
hist(X, freq=FALSE, breaks=c(seq(0,1,0.02),Inf), xlim=c(0,1))
lines(s,dbeta(s,5,3),type="l",main="PDF of Beta(5,3)", col="red",lty=2)
######### Assignment 5.3.g Binomial(10,0.4) #########
m <- 100000
Y = c()
for(i in 1:m){
X = c()
for(j in 1:10){
u <- runif(1,0,1)
p = 0.4
if (u < p)
X[j]=1
else if(u>=p)
X[j]=0
}
Y[i]= sum(X)
}
mean(Y)
## [1] 4.00237
median(Y)
## [1] 4
hist(Y, freq=FALSE, breaks=c(seq(0,15,0.25),Inf), xlim=c(0,15))
################# Assignment 5.4 ######################
######### Assignment 5.4.a #########
# Read in data
setwd('D:\\R Code')
normdata = read.table("norm.txt");
data <- c(normdata$V1)
data
## [1] 2.983 1.309 0.957 2.160 0.801 1.747 -0.274 1.071 2.094 2.215
## [11] 2.255 3.366 1.028 3.572 2.236 4.009 1.619 1.354 1.415 1.937
X <- mean(data)
nsim = 100000; # Total Number of Samples
sim = 0; # Indicator
samp = rep(NA,nsim); # Total Set of Samples
while(sim < nsim){
cand = rcauchy(1,0,1) # Candidate Generation
u = runif(1,0,1) # Comparison Gener
ratio = prod(dnorm(data,cand))/prod(dnorm(data,X)) # Ratio Calculation
if(u < ratio){ # Acceptance
sim = sim + 1
samp[sim] = cand
}
}
hist(samp,nclass=100)
plot(density(samp))
############ Assignment 5.4.b ##########
##### Generate Cauchy Distribution #####
library(VGAM)
## Loading required package: stats4
## Loading required package: splines
m = 100000; # Number of Candidates
n = 5000; # Number of Resample
y = rcauchy(m); # Random Number Generation from Cauchy Distribution
v = rnorm(m); # Random Number Generation from Normal Distribution
w.norm = function(x){ # Weight Calculation Function for a Normal Envelope
weight = dcauchy(x)/dnorm(x) # Importance Weights
weight = weight/sum(weight) # Standardization
return(weight)
}
weights.norm = w.norm(v) # Calculate Weights with a Normal Envelope
samp = sample(v,n,replace=TRUE,prob=weights.norm) # Resample with a Normal Envelope
hist(samp,freq=FALSE,breaks=seq(-6,6,by=.2),main="Histogram of Draws",ylab="Cauchy Density")
points(seq(-10,10,by=.01),dcauchy(seq(-10,10,by=.01)),type="l")
################ Assignment 5.5 ####################
######### Assignment 5.5.a #########
library(VGAM)
f <- function(x){dtriangle(x,0.25,0,1)}
invtri <- function(u){
if(u>=0 && u<0.25){
s = sqrt(u)/2
}
else{
s = 1-sqrt(3*(1-u))/2
}
return(s)
}
u <- seq(0,1,length.out=1000)
x <- seq(0,1,length.out=1000)
y = rep(0,length(u))
for (i in 1:length(u)){
y[i]=invtri(u[i])
}
X <- runif(100000,0,1) #random sample from U[0,1]
z = rep(0,length(X))
for (i in 1:length(X)){
z[i]=invtri(X[i])
}
plot(u,y,type="l",main="Inverse CDF of f(x)")
par(mfrow=c(1,1))
hist(z, freq=FALSE, breaks=c(seq(0,1,0.02),Inf), xlim=c(0,1))
lines(x,f(x),type="l", col="red")
######### Assignment 5.5.b #########
## g ~ Uniform(0,1), e~3g
f <- function(x){dtriangle(x,0.25,0,1)}
g <- function(x){1}
e <- function(x){3}
sample.x = runif(10000,0,1)
accept = c()
for(i in 1:length(sample.x)){
U = runif(1, 0, 1)
if(e(sample.x[i])*U <= f(sample.x[i])) {
accept[i] = 'Yes'
}
else if(e(sample.x[i])*U > f(sample.x[i])){
accept[i] = 'No'
}
}
x <- seq(0,1,length.out=1000)
T = data.frame(sample.x, accept = factor(accept, levels= c('Yes','No')))
hist(T[,1][T$accept=='Yes'], breaks = seq(0,1,0.02), freq = FALSE, main = 'Histogram of X', xlab = 'X')
lines(x,f(x),type="l", col="red")
################ Assignment 5.6 ####################
nsim = 10000
a = rnorm(nsim,20,1)
prob = length(a[a>=20])/nsim
x <- seq(15,25,0.1)
hist(a, breaks = seq(15,25,0.1), freq = FALSE, main = 'Histogram of X', xlab = 'X')
lines(x, dnorm(x,20,1),type="l", col="red")
