num.reps <- 1000
v <- 1000
thetaHat <- rep(0,num.reps)
thetaHati <- rep(0,num.reps)

for (i in 1:num.reps) {
  thetaHat[i] <- mean(sqrt(rgamma(v,3,2)));
  temp <- rexp(v,2/3);
  thetaHati[i] <- mean(sqrt(temp) * dgamma(temp,3,2) / dexp(temp,2/3))
}

mean(thetaHat)
sqrt(var(thetaHat))
mean(thetaHati)
sqrt(var(thetaHati))