## R code for last problem of STAT 512 Homework 8, Spring 2022:

# Generating data to get things started
#
# generating lots of samples of size 20 from a Unif(0.5*theta, 1.5*theta) distribution:

num_samps <- 100000
my_theta <- 100
samp_size <- 20

Y_mat <- matrix(runif(n=num_samps*samp_size, min=0.5*my_theta, max=1.5*my_theta),nrow=num_samps, ncol=samp_size, byrow=T)

##
## Code for part (b):
##
## Rao-Blackwellized estimator that is the sample midrange here:
##

theta_hat_RB <- function(vec){
est_RB <- 0.5*(max(vec) + min(vec))
return(est_RB)
}
# apply the theta_hat_RB function to each row of Y_mat
# Will return 'num_samps' values of theta_hat_RB:
theta_hat_RB_values <- apply(Y_mat,1,theta_hat_RB)

# Approximating the expected value of theta_hat_RB:
E_theta_hat_RB <- mean(theta_hat_RB_values)
print(paste("Approximate expected value of theta_hat_RB is: ", round(E_theta_hat_RB,4) ))

# Approximating the variance of theta_hat_RB:
V_theta_hat_RB <- var(theta_hat_RB_values)
print(paste("Approximate variance of theta_hat_RB is: ", round(V_theta_hat_RB,4) ))

##
## Code for part (c):
##
## LV estimator:
##

theta_hat_LV <- function(vec){
est_LV <-  (1/(0.25*((samp_size-1)/(samp_size+1))+1))*(0.5*(0.5*min(vec) + 1.5*max(vec)))
return(est_LV)
}
# apply the theta_hat_LV function to each row of Y_mat
# Will return 'num_samps' values of theta_hat_LV:
theta_hat_LV_values <- apply(Y_mat,1,theta_hat_LV)

# Approximating the expected value of theta_hat_LV:
E_theta_hat_LV <- mean(theta_hat_LV_values)
print(paste("Approximate expected value of theta_hat_LV is: ", round(E_theta_hat_LV,4) ))

# Approximating the variance of theta_hat_LV:
V_theta_hat_LV <- var(theta_hat_LV_values)
print(paste("Approximate variance of theta_hat_LV is: ", round(V_theta_hat_LV,4) ))

##
## Code for part (d):
##
## Bayes estimator:
##

theta_hat_Bayes <- function(vec){
numer<- (1.5*min(vec))/(0.5*max(vec))-1
denom<-( (1.5*min(vec)) / (0.5*max(vec)) )^(samp_size+1) - 1
est_Bayes <- ((samp_size+1)/(samp_size)) * (1- numer/denom) * (max(vec)/1.5)
return(est_Bayes)
}
# apply the theta_hat_Bayes function to each row of Y_mat
# Will return 'num_samps' values of theta_hat_Bayes:
theta_hat_Bayes_values <- apply(Y_mat,1,theta_hat_Bayes)

# Approximating the expected value of theta_hat_Bayes:
E_theta_hat_Bayes <- mean(theta_hat_Bayes_values)
print(paste("Approximate expected value of theta_hat_Bayes is: ", round(E_theta_hat_Bayes,4) ))

# Approximating the variance of theta_hat_Bayes:
V_theta_hat_Bayes <- var(theta_hat_Bayes_values)
print(paste("Approximate variance of theta_hat_Bayes is: ", round(V_theta_hat_Bayes,4) ))