Homework Solutions to Problems 7--8.
===================================
[NB solutions are in Splus 6.0]

(7) > motif()
> betvec <- rbeta(1600,2,3)
  hist(betvec, nclass=40, prob=T, xlab="Simulated variate scale",
    ylab="Scaled relative frequency", main=
   "Scaled Rel Freq Histogram: Beta Data & Overlaid Density")
  lines((1:200)/200, dbeta((1:200)/200,2,3), lty=2)

(8) (a) Here 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)") 
> lines((1:400)/100 -> 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(smooth.spline(Fx(xpt), xpt,
      spar=1.e-6), runif(1000))$y
  Ytmp2 <- FYinvcondX(runif(1000), Xtmp2)

(b) Here there are simple formulas: X = R*cos(Th), Y = R*sin(Th)
where R ~ F_R(r) = (2/pi)*arc tan(r^2) and Th ~ Unif(0,2*pi)
are independent.

> Thvec <- 2*pi*runif(1000)
  Rvec <- sqrt(tan(pi*runif(1000)/2))
  Xvec <- Rvec*cos(Thvec)
  Yvec <- Rvec*sin(Thvec)