Homework Problems 6, STAT 705, Fall 2015.

Assigned 9/30/2015,  Due Wednesday 10/7/2015

Give R code to generate (and then test the distribution in order 
to check your work) 10,000 pseudo-random variates from each of 
the two following densities

(1) the density on  [0,Infty)  proportional to  exp(-x) + exp(-x^2/2)

#---------- Solution to part (1): ----------------------------------------

The density is K*(exp(-x)+exp(-x^2/2)), where K = 1/(1+ sqrt(2*pi)*0.5) = 0.44379.
It can be viewed as a mixture with weights K and K*(sqrt(2*pi)/2) = (1-K) of the 
densities f(x) = exp(-x) and g(x) = (2/sqrt(2*pi))*exp(-x^2/2) on (0,Inf), and 
g(x) is the density of the absolute value of a standard normal r.v. So simulate by:


MixSim = function(N) {
  eps = rbinom(N,1,1/(1+ sqrt(2*pi)*0.5))
  out = numeric(N)
  m = sum(eps)
  out[eps==1] = rexp(m)
  out[eps==0] = abs(rnorm(N-m))
  out  }

K = 1/(1+ sqrt(2*pi)*0.5)

> tmp = MixSim(1e4)
> hist(tmp, nclass=40, prob=T, main=
     "Histogram of 10,000 Mixed Exp & Abs.Normal RVs")
> curve((exp(-x)+exp(-x^2/2))*K, 0,7, col="red",add=T)

Theoretical vs empirical mean and 2nd moment:

round(array(c(K*2, K*(2+sqrt(2*pi)/2), mean(tmp), mean(tmp^2)), c(2,2), 
  dimnames=list(c("Mean","2ndMoment"), c("Theoretical","Empirical"))),5)
          Theoretical Empirical
Mean          0.88758   0.88902
2ndMoment     1.44379   1.43769
> 
### Very close agreement

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


(2) the density on  [0,Infty)  proportional to  min(exp(-x), exp(-x^2/2), exp(-x^3/3)).

## ---------  Solution to part (2): -----------------------------------

The density is actually given in two pieces, since the min is 

exp(-x) for  x < sqrt(3),  exp(-x^3/3) on  (sqrt(3),Inf)

as you can check using:
> curve(exp(-x),0,7)
> curve(exp(-x^2/2),0,7, add=T, col="blue")
> curve(exp(-x^3/3),0,7, add=T, col="red")

## Normalizing constant for density is
> B = 1/integrate(function(x) pmin(exp(-x), exp(-x^3/3)),0,Inf)$val
## value 1.15066

so the desired density is a mixture with weights  B*(1-exp(-sqrt(3)) = 0.94708, 
and  B*integrate(function(x) exp(-x^3/3),sqrt(3),Inf)$val = 0.05292 of the 
respective densities exp(-x)/(1-exp(-sqrt(3)) on (0,sqrt(3)) and 
exp(-x^3/3)/0.04598772 on (sqrt(3),Inf).

## It is easy to simulate from exp(-x) on the first interval (0,sqrt(3)).

But to simulate from exp(-x^3/3) on (2,Inf), the simplest way may be accept-reject, 
simulating from x^2*exp(sqrt(3)-x^3/3), and making use of the uniform bound 
on (sqrt(3),Inf) of

(exp(-x^3/3)/0.04598772) / (x^2*exp(sqrt(3)-x^3/3)) < 

   exp(-sqrt(3))/(3*0.04598772) = 1.28238

Then the simulation method is:

MinSim = function(N) {
    K = exp(-sqrt(3))/(3*0.04598772)
    p = (1-exp(-sqrt(3)))/
            integrate(function(x) pmin(exp(-x), exp(-x^3/3)),0,Inf)$val
    eps=rbinom(N,1,p)
    out = numeric(N)
    out[eps==1] = -log(1-runif(sum(eps))*(1-exp(-sqrt(3))))
    m = N-sum(eps)
    mbound = max(round(2.5*m),1000)
    W = (-3*log(exp(-sqrt(3))*runif(mbound)))^(1/3)
    tmp = W[runif(mbound)<3/W^2] 
    if(length(tmp)<m) {
        Walt = (-3*log(exp(-sqrt(3))*runif(mbound)))^(1/3)
        tmp = c(tmp, Walt[runif(mbound)<3/Walt^2]) }
    out[eps==0] = tmp[1:m]
    out
   }

tmp = MinSim(1e4)    
### if there are any 0's in tmp, then we did not generate enough of the second
###  accept-reject set of random variables values. But the number of those 
###  that were generated is  2.5 times the expected number needed, and then if
###  that is not enough (as may rarely happen), the function generates another 
###  batch of 2.5 times the expected number. By the Central limit Theorem, when
###  N is large the probability is really vanishingly small that enough rv's 
###  are generated, but if N is small it could happen.

> sum(tmp==0)
[1] 0

> hist(tmp, nclass=40, prob=T, main="Histogram of 10,000 Min-Density RVs")
  curve(B*exp(-x),0,sqrt(3), add=T, col="green")
  curve(B*exp(-x^3/3),sqrt(3),max(tmp), add=T, col="green")
### very close agreement, check again with larger N

tmp = MinSim(1e5)  
hist(tmp, nclass=50, prob=T, main="Histogram of 10,000 Min-Density RVs")
  curve(B*exp(-x),0,sqrt(3), add=T, col="green")
  curve(B*exp(-x^3/3),sqrt(3),max(tmp), add=T, col="green")