Homework Problem 5 Solution
---------------------------

(a) 
> fxy = function(x,y) sqrt(1+ x*y + y^2 + x^4)
  fvals = fxy(runif(1e5),runif(1e5))
  ests= c(mean(fvals), sd(fvals)/sqrt(1e5))
  c(Integral=ests[1], CI.lo=ests[1] - qnorm(.995)*ests[2], 
                 CI.hi=ests[1] + qnorm(.995)*ests[2])
Integral    CI.lo    CI.hi
1.316910 1.315139 1.318680   ### 99% CI accuracy .0018

> mean(fxy(runif(1e6),runif(1e6)))  ## 2nd confirmatory run
[1] 1.317606

In fact, the integral of fxy over y from 0 to 1 could be evaluated
  analytically as a function of x, but we do the integral instead 
  as an iterated numerical integral:

  gy = function(y) {
         out = numeric(length(y))
         for(i in 1:length(y)) out[i] =
                  integrate(fxy,0,1,rel.tol=1e-6, y=y[i])$val
         out }
> integrate(gy,0,1, rel.tol=1.e-6)$val
[1] 1.317694


(b) Let M be the random variable value, with values 1:5.
> table(apply(array(sample(0:8, 5e6, rep=T), dim=c(1e6,5)),1,
          function(irow) length(unique(irow))))/1e6
       1        2        3        4        5
0.000169 0.018601 0.213541 0.511832 0.255857
> pvec=.Last.value


### This took 2 minutes on my office PC
### Note that the 99.6% CI half-widths are  
> round(qnorm(.998)*sqrt(pvec*(1-pvec)/1e6),5)
      1       2       3       4       5
0.00004 0.00039 0.00118 0.00144 0.00126

### So we can say that we achieved approximately these 
### accuracies in our estimated probability masses,
### with probability of ANY of the masses being outside 
### its interval equal at most to 5*(1-.996) = .02, as desired.
(This simultaneous confidence interval idea is essentially the idea of
using "Bonferroni bounds" to make simultaneous confidence statements.)
### You can see that to choose n in advance, we would have needed
sqrt(.5*.5)*qnorm(.998)/sqrt(n) < .01 ,   or  n   at least
(100*.5*qnorm(.998))^2  = 20710

### In fact, we can also check the probabilities of these various
    numbers of distinct values, by exact combinatorial calculation:

P(1 distinct value) = 9*(1/9)^5 = 0.00015
P(2 distinct values) = 
     9*8*5*(1/9)^5 + 9*8*(5 choose 2)*(1/9)^5 =  0.01829
P(5 distinct values) = 9*8*7*6*5/9^5 = 0.25606
P(4 distinct values) = 9*8*7*6*(5 choose 2)/9^5 = 0.51212
and therefore the remaining probability mass is
P(3 distinct values) = 1-0.00015-0.01829-0.25606-0.51212 = 0.21338