Solutions to HW Set 6
---------------------

### (a) Here   f_X(x) = (2-(x-1)^2) / integral of 2-(x-1)^2 on [0,2]
##                    = 0.3*(1+2*x-x^2)

##          F_X(x) = 0.3*(-x^3/3 + x^2 + x)    on [0,2]    and  
##    f_{Y|X}(y|x) = 1/( 2-(x-1)^2 ) ~ Unif(0, 2-(x-1)^2)

> unlist(integrate(function(x)  2-(x-1)^2,0,2)[1:2])
       value    abs.error 
3.333333e+00 3.700743e-14 

## so cdf F_X =  -x^3/3 + x^2 + x    on [0,2]
## Simplest method is Accept-reject !!

Xbig = runif(2e5,0,2)
Ybig = runif(2e5,0,2)
inds = which( Ybig <= 1+2*Xbig-Xbig^2 )[1:1e4]
Mout = cbind(Xbig[inds], Ybig[inds])
> plot(Mout)         #### shows that the desired region is filled !!

## Second method: inversion of the df requires solving a cubic via uniroot.

FXptA = function(x,u) 0.3*(x + x^2-x^3/3) - u

PairA = function(N) { 
           uvec = runif(N)
           out = numeric(N)
           for(i in 1:N)          ### very slow !!
               out[i] = uniroot(FXptA, c(0,2), u=uvec[i])$root
           cbind(out, runif(N, 0, 2-(out-1)^2)) }

> plot(PairA(1e3))     ### again fills the desired region visually uniformly


### (b) Here   F(x) = (1/24)*((1+x^2)^2 - 1) ,  so the inverse function
###     is     F^(-1)(r) = sqrt( sqrt(24*r+1) - 1 )     for   0<=r<=1

PtB = function(N) sqrt( sqrt(24*runif(N)+1) - 1 )

### Here is a small check: we want the fraction of r.v. values between 
###  0.5 and 1 to be  F(1)-F(.5) = (1/24)*(4 - 1.25^2) = 0.1015625 and 
###  the fraction between 1.5 and 1.9   to be 
###       (1/24)*((1+1.9^2)^2 - (1+2.25)^2) = .4454

tmp = PtB(1e6)
> mean( (tmp>0.5)*(tmp<=1) )
[1] 0.100975
> mean((tmp>1.5)*(tmp<=1.9))
[1] 0.445823


 ### (c) Here a quick change of variable Y = sqrt(X) gives f_Y(y) = 2C y exp(-y)
###     So Y ~ Gamma(2,1)    and  C = 0.5

### One method is  :   PtC = function(N) rgamma(N,2,1)^2

### Second method : integration (and linear interpolation between points)

> M = 1e3; tol = 1e-6; N=1e4
  ztmp = seq(tol,1, by=(1-tol)/M)
  xtmp = rev(log( ztmp )^2)
  Ftmp = numeric(M+1)
  tmpf = function(x) exp(-sqrt(x))
  for(i in 1:M) Ftmp[i+1] = integrate(tmpf,xtmp[i],xtmp[i+1])$val
  Ftmp = 0.5*cumsum(Ftmp)
### now get inverse at all points in [0,1] by linear interpolation between
### the points with x-coords  Ftmp  and  y-coords  ztmp = seq(tol,1, by=(1-tol)/M)
### Recall function "LinInterpB" from the solution to the first part of HW4:
  PtC = LinInterpB(cbind(x=Ftmp, y=xtmp), runif(N))[,2]
> hist(PtC, nclass=50, prob=T)
  lines(xtmp, 0.5*exp(-sqrt(xtmp)), col="red", lty=2)
### Low-resolution picture, but clearly a good fit !! Same as for 1st method.


## (d)     f(x, y) = C exp(-x -y -2xy)

 ##     f_X(x) = C/(1+2x) exp(-x)
   and   f_{Y|X}(y|x) ~ Expon(1+2x),   but

> Cons = 1/integrate(function(x) exp(-x)/(1+2*x),0,Inf)$val
> Cons
[1] 2.167057
> xt = rev(-log(seq(1e-6,1,length=1001)))

### xt values are generated so that the exp(-xt) points are equally spaced

  Ft = numeric(1001)
  for(i in 1:1000) Ft[i+1] = integrate(function(x)
           Cons*exp(-x)/(1+2*x),xt[i],xt[i+1])$val
  Ft = cumsum(Ft)

### The main point here is to understand that since  (xt[i], Ft[i])
### are points on the cdf curve, the reversed points (Ft[i], xt[i])
### are points on the inverse cdf curve !! This is a trick we will use again!

  PtD = LinInterpB(cbind(x=Ft,y=xt), runif(1e5))

### This vector  PtD  is a column of generated X's according to f_X density.
### Then the matrix output of X's and associate Y's would be

 XYmat = cbind(newx= PtD, newy = rexp(1e5, 1+2*xt))


## In (a) the Accept-Reject method is easy to do without for-loop; 
##  the other method is not. In (b) the function naturally vectorizes.
## In (c) the method by change of variable and rgamma naturally vectorizes,
##  but the other does not unless the integration is programmed to vectorize.
## Similarly, in (d), the repeated integration could vectorize if numerical 
##  quadrature  were newly programmed, but repeated integration and uniroot does not!