Homework Set 6 Solution.                     
========================

(8) (a) The point is that, in using the "bivariate probability
    integral transform" method discussed in class, you can work
    either with F_X and F_{Y|X} or with F_Y and F_{X|Y}. If the
    objective is to do everything analytically, then F_Y and F_{X|Y}
    works better, so that will be the way I do it first (Method 1).
    The other way works too, but uses a little more in the way of
    Splus trickery because F_X, while simple, cannot be analytically
    inverted. That will be Method 2, and I use the two methods to
    check each other.

## First we exhibit the marginal and conditional df's.
> Fx = function(x) 1 - 1.5/x^2 + .5/x^3
  Fy = function(y) 1-1.5/y + .5/y^2
  Finvy = function(r) .75*(1 + sqrt(1-8*(1-r)/9))/(1-r)
  FYcondX = function(x,y) ifelse(y <= x,
    (y-1)/(2*x-1), 1 - x^2/(y*(2*x-1)))
  FXcondY = function(x,y) ifelse(x <= y,
    1.5*(1-1/x)/(1.5-1/y), 1 - .5*(y/x)^3/(1.5*y-1))
  FYinvcondX = function(r,x) ifelse(r <= FYcondX(x,x),
    (2*x-1)*r+1, x^2/((2*x-1)*(1-r)))
  FXinvcondY = function(r,y) ifelse(r <= FXcondY(y,y),
    1/(1-r*(1-2/(3*y))), y*(.5/((1-r)*(1.5*y-1)))^(1/3))

## NB we can test these inverses simply, eg by generating
##   arbitrary X values on (1,infty) and R on (0,1) and
##   then making sure that FYcondX(FYinvcondX(R,X),X)=R
> Xtmp = exp(rexp(100)); Rtmp = runif(100)

> sum(abs(FYcondX(Xtmp, FYinvcondX(Rtmp,Xtmp))-Rtmp))
[1] 1.259409e-15    ### this is =0, OK
> sum(abs(FXcondY(FXinvcondY(Rtmp,Xtmp),Xtmp)-Rtmp))
[1] 5.972653e-15    ### OK again
> sum(abs(Fy(Finvy(Rtmp))-Rtmp))
[1] 2.664535e-15    ### OK again

METHOD 1.
> Ytmp = Finvy(runif(1000))
  Xtmp = FXinvcondY(runif(1000),Ytmp)
### We display the histogram of log(Xtmp) versus
###    the theoretical density exp(z)*fx(exp(z))
> hist(log(Xtmp), nclass=36, prob=T, main=
     "Scaled Histogram of log(X)")
> lxp = (1:400)/100
  lines(lxp, 3*exp(-2*lxp)-1.5*exp(-3*lxp), lty=3)
### Virtually perfect agreement
> hist(log(Ytmp), nclass=36, prob=T, main=
     "Scaled Histogram of log(Y)")
  lines((1:400)/100 -> lxp, 1.5*exp(-lxp)-exp(-2*lxp),
     lty=3)   ### Again very close agreement

METHOD 2.
> xpt = exp((0:400)/80)
> Xtmp2 = predict(smooth.spline(Fx(xpt), xpt,
      spar=1.e-6), runif(1000))$y
  Ytmp2 = FYinvcondX(runif(1000), Xtmp2)


(b) The region of interest is the set in the plane with
y-coordinates ranging from 0 to 4, and with x ranging from
(-3y/4) to  6-9y/4 . It is easy to see that the density
of Y is not uniform but rather is proportional to 
6-9y/4 - (-3y/4) = 6-3y/2, and therefore

fY(y) = (1/2)(1-y/4)  for  0<y<4 , 0 elsewhere

FXcondY(x|y) = (6-3*y/2)^(-1)  for  (-3*y/4) < x < 6-9*y/4
     and  0 elsewhere.

Therefore for  0<r<1:

 FYinv(r) = 4*(1-sqrt(1-r)) ,   and

FXinvcondY(r|y) = (-3*y/4) + r*(6-3*y/2)

> FYinv
function(r) 4*(1-sqrt(1-r))
> FXinvcondY
function(r,y) (-3*y/4) + r*(6-3*y/2)

Yvec = FYinv(runif(1e5))
Xvec = FXinvcondY(runif(1e5), Yvec)

> mean(Xvec< 0)
[1] 0.33322        ### Theoretical prob 1/3
> mean(Yvec > 2)
[1] 0.24919    ### Theoretical prob 1/4
> mean(abs(Xvec-2.5)<.5 & abs(Yvec-.77)<.5)
[1] 0.08241    ### Theoretical prob 1/12

> plot(Xvec[10*(1:1000)],Yvec[10*(1:1000)])
### Looks nice and uniform

The rejection method is easier to program, but wasteful:

## To find uniform in the right triangle with base (-3,6)
## and vertical leg from (-3,0) to (-3,4):

Xvec = Xvec0 = 6-9*runif(1e5)
Yvec = Yvec0 = 4*runif(1e5)
ind = (Yvec0 > (4/9)*(6-Xvec0))
Xvec[ind] = 3-Xvec0[ind]
Yvec[ind] = 4-Yvec0[ind]
> ind2 = (Yvec > -4*Xvec/3)
> mean(ind2)
[1] 0.66759   ## This is the fraction of points we can use !!

> plot(Xvec[ind2][1:1000],Yvec[ind2][1:1000])
> mean(abs(Xvec[ind2]-2.5)<.5 & abs(Yvec[ind2]-.77)<.5)
[1] 0.08419839   ### compare to 1/12
> mean(Xvec[ind2] < 0)
[1] 0.3334831
> mean(Yvec[ind2] > 2)
[1] 0.2490301                      ### All OK !!!

### Various methods can be used to check the uniformity in
### the triangle, but some credit will be taken off 
(e.g., 1 point) if the student does nothing beyond looking at
a scatterplot.