Homework Problem 7 solutions,     10/6/2016.
--------------------------------------------

(1) 10,000 random points from the region in the plane containing (0,0) 
and bounded by the lines

     y = 3 + 2x,     y = 2 -7x   and  y= 2,    y=-5.

SOLUTION: 

this region is {(x,y):  -5 <= y <= 2,  (y-3)/2 <= x <= (2-y)/7}
so that 

f_{X|Y}(x) = 1/((2-y)/7 - (y-3)/2) = 14/(25-9*y)  for  (y-3)/2 <= x <= (2-y)/7

and  (since  f(x,y) = 1/C on the trapezoid and 0 elsewhere, 
with C=(1/2)(7)(1/2 + 5) = 19.25, we have 

f_Y(y) = (25-9*y)/(14*19.25)    for  -5 < y < 2

and the cdf is

F_Y(y) = (1/269.5)*(25*(y+5) - (9/2)*(y^2-25)) ,   -5 < y < 2

and its inverse, by the quadratic formula, is

F_Y^{-1}(u) = (1/9)*(25 - sqrt(25^2 + 4*(9/2)*((9/2)*25 + 125 -269.5*u)))
            = (1/9)*(25 - sqrt(4900-4851*u))

Then, to generate a pair (X,Y) with the required density, take two independent
Unif[0,1] variates U_1 and U_2, and put

Y = F_Y^{-1}(U_1) = (1/9)*(25 - sqrt(4900-4851*U_1)
X = (Y-3)/2 + ((25-9*Y)/14)*U_2

or, in R code

> FY = function(y) (25*(y+5)- 4.5*(y^2-25))/269.5
  FYinv = function(u) (1/9)*(25 - sqrt(4900-4851*u))
> sum(abs(FY(FYinv(x))-x))
[1] 5.377643e-15
  #### Verifies inverse property !!

> Yvec = (25  - sqrt(4900-4851*runif(1e4)))/9
  Xvec = (Yvec-3)/2 + ((25-9*Yvec)/14)*runif(1e4)
> plot(Xvec,Yvec)

YOU COULD ALSO HAVE DONE THIS PROBLEM (WITHOUT TOO MUCH WASTE BY AN 
ACCEPT-REJECT  METHOD BY EMBEDDING THE TRAPEzOID within the rectangle
[-4,1] x [-5,2] .

#-------------------

(2) 10,000 independent random variable variates on (0,5) with 
univariate density 

     g(x) = C * (2 + 0.2*sin(3*x^2))   ,   for  x in (0,5)
          = 0                             elsewhere

for a constant C which makes this a proper density.

HINT: use an accept-reject method.

SOLUTION:

## Here the constant C is 
> Const = 1/(10+0.2*integrate(function(x) sin(3*x^2),0,5)$val)
[1] 0.09934203

## Method is to generate pairs W~Unif[0,5], U ~ Unif[0,2.2*Const]
## and accept X=W whenever    U <= g(W)

## NB The probability of this happening is
> g = function(x) Const*(2+0.2*sin(3*x^2))
> integrate(function(w) g(w)^2, 0,5)$val/(2.2*Const)
[1] 0.9193678

## So to get 10,000, it ought to be enough to generate 1.1*1e4/0.91937=11965
## (W,U) pairs

  W = runif(11965,0,5)
  U = runif(11965,0,2.2*Const)
  ind = (U < g(W))
> sum(ind)
[1] 10959           ### so we have more than enough
> X = W[ind][1:1e4]
> hist(X, nclass=50, prob=T)
> curve(g(x), 0, 5, add=T, col="green")  ### reasonably close !

##----------------------------------

(3) 1,000 independent random pairs (X[i],Y[i]) on (0,Inf) x (0,Inf)
with joint density

     f(x,y) = K y/(1+x^2+y^2)^4      for x, y > 0

where K is a constant which makes this a proper density.

SOLUTION:  

## Note that by subsitution integration (z = 1+x^2+y^2) in the dy integral,
##   the analytical f(x) density = (K/6)/(1+x^2)^3   for x>0, 
##   which happens to make  X  equal to abs(W)/sqrt(5) for W ~ t_5, i.e.,
##   W a student's t random variable with 5 degrees of freedom.

Then we do a partial integral in the y variable to find

integral from 0 to t of f(y|x) dy = F(y|x) = 1 - ((1+x^2)/(1+x^2+y^2))^3

To invert this, we solve  u = F(y|x) to give

        y = F^{-1}(u|x) = sqrt(  (1+x^2)*((1-u)^(-1/3)-1) )

R code for generating the pairs is then

X = abs(rt(1e4,5))/sqrt(5)
Y = sqrt((1+X^2)*((1-runif(1e4))^(-1/3)-1))

### We can check this by doing a couple of probabilities both 
###    empirically and `theoretically', by integration.

> K = 6/integrate(function(x) 1/(1+x^2)^3,0,Inf)$val
[1] 10.18592

## (You can check that the constant factor K/6 of the density
##   is exactly equal to the value at 0 of the density of abs(W)/sqrt(5),
##   ie K/6 = 2*dt(0,5)*sqrt(5).)

## First, P(X<1, Y<1):

> mean( X<1 & Y<1 ) 
[1] 0.7812

> (K/6)*integrate(function(x) 1/(1+x^2)^3 - 1/(1+x^2+1)^3,0,1)$val
[1] 0.7785172

### A second check is P(X<0.5, Y<2)
> mean( X<0.5 & Y<2)
[1] 0.6793
> (K/6)*integrate(function(x) 1/(1+x^2)^3 - 1/(1+x^2+2^2)^3,0,.5)$val
[1] 0.6791571

### Many more such empirical versus theoretical probabilities can 
###  be produced. Another test would be to look at the theoretical 
###  versus empirical curve for P(Y>=X>=x)  which is related to the t_5
###  distribution.