
# program to estimate E[2^(Y+Z/2)] where Y,Z ~ Unif{1,2,...,10}

# choose simulation size
n = 10000

# make a vector of values
Y = sample(10, n, replace=TRUE)
Z = sample(10, n, replace=TRUE)
X = 2^(Y+Z/2)

# compute and output the mean, sd, and confidence interval
m = mean(X)
se = sd(X) / sqrt(n)
cat("MC:  ", m, " +- ", se, "  (n=", n, ")", "\n", sep='')  # \n is "newline"
cat("  95% C.I.:  (", m-1.96*se, ",", m+1.96*se, ")\n", sep='')

# compare to exact answer
thesum = 0
for (j in 1:10)
  for (k in 1:10)
    thesum = thesum + 2^(j+k/2)
print(thesum/100)

