HW Problem 9 Solution
===================== 10/7/09

Xarr = array(rlogis(2e6), c(1e5,20))
Yarr = Xarr*array(runif(2e6,.75,1.25), c(1e5,20))
X10 = apply(Xarr, 1, function(xrow) sort(xrow)[10])
Y10 = apply(Yarr, 1, function(xrow) sort(xrow)[10])
cval = integrate(function(x) dbeta(x,10,11)*(qlogis(x))^2,0,1)$val
> round(c(theor=cval, empir=mean(X10^2)) ,5)
  theor   empir 
0.21033 0.20953    ### OK, these are close !

Val1 = mean((Y10-X10)^2)
> c(naive.est=Val1, naive.SE = sqrt(var((Y10-X10)^2)/1e5))
   naive.est     naive.SE 
3.297885e-03 2.289494e-05 

tmp = lm(I(Y10-X10)^2 ~ I(X10^2-cval))
> summary(tmp)$coef[1,1:2]
    Estimate   Std. Error 
3.309335e-03 1.825161e-05        ### Improvement in SE
   #### shows that with control variates we would need
### for the same CI width a number of replicates R=
> 1e5*(1.825/2.289)^2
[1] 63567.36             ### a significant saving, 36% !!