Problem 1

for the first problem we have a polinomy of high power, which means that finding the inverse function for the cumulative function is really complex, for that reason we’ll be using the rejection method

In this case \(f(x) = 30\cdot(x^2-2x^3+x^4)\)

Since we are working in the [0,1] interval the candidate density g is going to be the uniform distribution between 0 and 1.

\[ \left\{ \begin{array}{ll} 1 & 0 \leq x \leq 1\\ 0 & otherwise \\ \end{array} \right. \]

rRejP1 = function(size = 1000){
  
  
  x = rep(0,size)
  
  for(i in 1:size){
    
    u = runif(1)
    y = runif(1)
    c = 1.875 # This is the upper bound of the likelihood ratio f(x)/g(x)
    
    while (u > (30*(y^2-2*y^3+y^4)/(c*1))){ # This equation is f(x)/cg(x)
      
      u = runif(1)
      y = runif(1)
      
      
    }
    
    x[i] = y
  }
  
  return(x)
}

Next we obtain the 1000 samples using the rejection method and draw the histogram with the actual density of the random variable

samp = rRejP1(1000)

hist(samp,freq = F,main = "Histogram with 1000 samples")
curve(30*(x^2-2*x^3+x^4),add = T,col = "red")

The red line corresponds to the actual density of the given function which seems to fit fairly well with our random sample.

Problem 2

Now, for this case, we want to use both the rejection and the inverse transform method.

Inverse transform

First we will proceed with the inverse transform

We know that the given \(f(x) = 2xe^{-x^2}\) for x > 0

Firstly, we want to find the cumulative function by integrating f(x) between 0 and t.

This results on \(F(t) = 1 - e^{-t^2}\).

Next we calculate the inverse function of F.

Resulting in \(F^{-1}(u) = \sqrt{-ln(1-u)}\)

rInvTP2 = function(size = 1000){
  
  u = runif(size)
  
  y = sqrt(-log(1-u))
  
  return(y)
  
}

Next we obtain the 1000 samples using the inverse transform method and draw the histogram with the actual density of the random variable.

samp2 = rInvTP2(1000)

hist(samp2,freq = F,main = "Histogram with 1000 samples\n (Inverse transform)",ylim = c(0,0.9))
curve(2*x*exp(-x^2),add = T,col = "red")

The red line corresponds to the actual density of the given function which seems to fit fairly well with our random sample.

Rejection method

For this method we first choose a candidate density g

Since the given f(x), as we saw before, is defined in x > 0, then a good candidate density would be the exponential distribution (with \(\lambda = 1\) in this case), so \(g(x) = e^{-x}\)

The upper limit c for the likelihood ratio \(\frac{f(x)}{g(x)} = 2xe^{x(1-x)}\) is 2

the resulting \(\frac{f(x)}{c\cdot g(x)} = x\cdot e^{x(1-x)}\)

rRejP2 = function(size = 1000){
  
  x = rep(0,size)
  
  for (i in 1:size){
    
    u = runif(1)
    y = rexp(1)
    
    
    while(u > (y*exp(y*(1-y)))){
      
      u = runif(1)
      y = rexp(1)
      
    }
    x[i] = y
  }
  
  return(x)
  
}

Next we obtain the 1000 samples using the rejection method and draw the histogram with the actual density of the random variable.

samp3 = rRejP2(1000)
hist(samp3,freq = F,main = "Histogram with 1000 samples\n (Rejection)",
     ylim = c(0,0.9))
curve(2*x*exp(-x^2),add = T,col = "red")

The red line corresponds to the actual density of the given function which seems to fit fairly well with our random sample.