0
$\begingroup$

I'm trying to complete an exercise in R related to Monte Carlo estimation for the shortest path in the bridge network problem. The exercise asks me to first perform the estimation using the plain Monte Carlo method by generating the variable according to an exponential distribution, and I don't have any issues with that. I'm facing problems in the second part of the exercise where I'm asked to apply the method of antithetic variates. When using the qexp() function, I'm obviously encountering errors. Any advice on how to proceed? Here is my code:

lambda = 1
simulations= 10000
replicas = 1000
H.X = function(simulations, lambda, replicas){
  H_i = array(NA, dim = rip)
  for (i in 1:rip) {
    X = matrix(NA, nrow = sim, ncol = 5)
    for (j in 1:5) {
      X[,j] = rexp(sim, lam)
      H = pmin(X[,1] + X[,4], X[,1] + X[,3] + X[,5], X[,2] + X[,3] + X[,4], 
               X[,2] + X[,5])
      l = mean(H)
    }
    H_i[i] = mean(l)
  }
  return(H_i)
}
$\endgroup$

1 Answer 1

1
$\begingroup$

For a continuous distribution, $Y= F^{-1}\left(1-F(X)\right)$ provides the antithetic random variable to $X$, using the CDF and its inverse the quantile function, assuming they are in fact usable.

So given $X \sim \textrm{Exp}(\lambda)$, with $F(x)=1-e^{-\lambda x}$ for $x\ge 0$, and $F^{-1}(p)=- \frac{1}\lambda \log_e (1-p)$,

then $Y=- \frac{1}\lambda \log_e (1-e^{- \lambda X})$ proves the antithetic random variable, with $Y \sim \textrm{Exp}(\lambda)$ too.

So you could use

X[,j] = rexp(sim, lam)
Y[,j] = qexp(1 - pexp(sim, lam), lam)

or explicitly

Y[,j] = -(1/lam) * log(1 - exp(-lam * X[,j]))

Alternatively, for a continuous distribution, $U=F(X)$ has a uniform distribution with support on $[0,1]$.

You can generate $U$ first, letting $X=F^{-1}(U)$ and its antithetic random variable $Y=F^{-1}(1-U)$, just using the quantile function or inverse CDF, assuming it is in fact usable.

This quantile function approach also works with discrete distributions, avoiding the problem of the discontinuous CDF.

So you could use

U = runif(sim)
X[,j] = qexp(U, lam)
Y[,j] = qexp(1 - U, lam)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.