#Question 1
#Part a, the tight upper bound was found for when x=1/5=0.2
#Part c, code with the correct tight upper bound
betarej<-function(n){
#Function to sample from Beta(2,5) distribution via rejection sampling from a uniform normal
#n the number of random samples from Beta(2,5)
z<-rep(NA,n) #This will store our Beta samples
M<-0.2*(1-0.2)^4 #The correct optial tight upper bound
for(k in 1:n){
z[k]<-runif(1) #The sample from the proposal distribution
u<- runif(1) #For the accept with probability part
while(u>(z[k]*(1-z[k])^4)/M){ #With the correct acceptance probability
z[k]<-runif(1) #If u>acceptance prob, reject the proposed sample and sample again from the proposal distribution
u<-runif(1)
}
}
return(z)
}
#Part d
z<-betarej(n=1000)
hist(z, prob=T)
curve(dbeta(seq(0,1, len=50), 2,5), add=T)
#Part e
beta(2,5)/(0.2*0.8^4)
betarej.acc<-function(n){
#Function to sample from Beta(2,5) distribution via rejection sampling from a uniform normal
#n the number of random samples from Beta(2,5)
c<-rep(NA,n) #Vector where each element denotes how many samples from the proposal were used to obtain sample from the target
z<-rep(NA,n) #This will store our Beta samples
M<-0.2*(1-0.2)^4 #The correct optial tight upper bound
for(k in 1:n){
z[k]<-runif(1) #The sample from the proposal distribution
u<- runif(1) #For the accept with probability part
c[k]<-1
while(u>(z[k]*(1-z[k])^4)/M){ #With the correct acceptance probability
z[k]<-runif(1) #If u>acceptance prob, reject the proposed sample and sample again from the proposal distribution
c[k]<-c[k]+1 #As we resampled from the proposal distribution
u<-runif(1)
}
}
return(1/mean(c)) #To return the average acceptance rate
}
betarej.acc(n=1000)
#Question 2
#Part a
pT<-exp(-10) # Theoretical mean of the estimator
1-pexp(10, rate=10) #What it should be according CDF
var.pT<-exp(-10)*(1-exp(-10)) /10^4
cv.pT<-sqrt(var.pT)/pT
#Part b
p_direct<-function(k,N,T){
#k estimates of the probability of interest p(T)
#N samples used from the exponential distribution
data<-matrix(NA,k,N)
p<-rep(NA,k)
for(i in 1:k){
data[i,]<-rexp(N) #Sample from the distribution
p[i]<-sum(data[i,]>T)/N #Proportion,estimate, of samples which satisfy the event of interest.
}
return(p)
}
direct.est<-p_direct(k=100, N=10^4, T=10) #Estimation by sampling directly
mean(direct.est)
sd(direct.est)
sd(direct.est)/mean(direct.est)
#Part c
p_importance<-function(k,N,T){
#Function to estimate via importance sampling
#k= number of estimates
#N = number of samples per estimate from importance distribution
#T = the event threshold
data<-matrix(NA, N,2)
p<-rep(NA,k)
for(i in 1:k){
data[,1]<-rexp(N, rate=1/8) #Sample from the importance distribution
data[,2]<-8*exp(-data[,1]*7/8) #The corresponding importance weight
p[i] <-sum((data[,1]>T)*data[,2]/N) #The estimate with respect to the event and the importance weight
}
p
}
p_est<-p_importance(1, 10^4, 10)
p_est
#Part d
p_est<-p_importance(100, 10^4, 10)
p_mean<-mean(p_est)
p_se<-sd(p_est)
cv<-p_se/p_mean