Antithetic method for monte carlo when bounds of the integral are infinite I wanted to apply Monte Carlo with antithetic variables to estimate $\int_{0}^{\infty} e^{-x} \,dx$ (equal to 1). I used this R code.
a <- runif(1000)
f <- function(x) exp(-x)
x <- f(a)
y <- f(1-a)
mean((x+y)/2)

However, it returned 0.632577961884066, which is not 1. I am really confused.

The point to the method of antithetic variates is to improve on such direct sampling methods.  What can make it work well is to re-express the integral as an expectation of something with respect to a distribution that (a) you can efficiently sample from and (b) samples the largest absolute values of the integrand with relatively high probability.
This distribution automatically finesses the problem of an infinite integration interval by assigning low probabilities to high values.
Criterion (b) means the density function of your chosen distribution should look a little like the integrand.
In mathematical notation, this means finding a distribution with density $f$ and re-expressing the integral as
$$\int_0^\infty e^{-x}\,\mathrm{d}x = \int_0^\infty \frac{e^{-x}}{f(x)}\, f(x)\mathrm{d}x.$$
This is the expectation of $e^{-x}/f(x)$ with respect to the new distribution.
Finding such a density $f$ takes some creativity but it can really pay off.  In the present instance, suppose you didn't know how to integrate exponential functions, so that "$e^{-x}$" is just some kind of black box, perhaps an expensive one to consult at that.  But let's suppose you do have some elementary command of integration; say, the ability to integrate powers (the very first thing learned in integral Calculus).  We might then try something like
$$f(x) \ \propto \ x^r$$
for some power $r.$  Its integral is
$$F(x) = \frac{x^{r+1}}{r+1}.$$
For this to remain finite for large $x$ we will need $r+1 \lt 0;$ that is, $r \lt 1.$  But there's a problem at $x=0:$ both $f$ and $F$ will blow up.  But why start integrating at $0,$ since that's a problem?  Let's start, say, at $1.$  What this amounts to is choosing
$$f(x) \ \propto \ \frac{1}{(1+x)^p}$$
for some power $p \gt 1$ (a Pareto distribution).  The integral of this $f$ is
$$F(x) = 1 - \frac{1}{(1+x)^{p-1}}.$$
Consequently
$$\frac{e^{-x}}{f(x)} = \frac{e^{-x} (1 + x)^p}{p-1}\tag{*}$$
and you can also easily solve the equation $F(x) = q$ (for $0\le  q \lt 1$) to give
$$x = F^{-1}(q) = \left(1-q\right)^{1/(1-p)} - 1\tag{**}.$$
This means you can generate random values $x_i$ from $F$ by generating $q$ from a uniform distribution and applying $(**)$ to it.  The antithetic variable is obtained by applying $(**)$ to $1-q,$ as shown in the question.
Let's compare the original integrand to this density $f:$

They are qualitatively similar,  implying we ought to gain something for our efforts.
To sum up this analysis, you will draw uniform independent values $q_1, q_2, \ldots, q_n$ from which you will compute $x_i$ and antithetic values $x_i^\prime$ using $(**)$ and, from those, apply $(*)$ to obtain values $y_i$ and $y_i^\prime,$ averaging them as $(y_i+y^\prime_i)/2$ to obtain one estimate of the integral.  By repeating this you can estimate the integral with their average
$$\frac{1}{n}\sum_{i=1}^n \frac{y_i + y^\prime_i}{2}.$$
We can do even better: we can use tiny $n$ to estimate the integral and study how variable the resulting estimate is when we alter the parameter $p.$  When I did this, I found that $p\approx 3$ gives estimates with the lowest variance.
Here, to illustrate, are the results of 50,000 estimates of the integral based on just $n=10$ draws.  On the left is the solution for $p=3$ and on the right is the solution based on drawing $x$ uniformly from the interval $[0, \log 2^{52}]$ (on which $1 + e^{-x}$ is distinguishable from $1$).

Although both are centered around the correct value of $1,$ the Pareto solution tends to be much closer to this value than the Uniform solution.  Indeed, because the ratio of their variances is $235:1,$ the Uniform solution requires $235$ times as many draws as the Pareto distribution in order to estimate the integral with the same precision.
Specifically, this example indicates you can reliably estimate the integral to within 20% of its true value with just $2n=20$ evaluations of the exponential function, while using the uniform method would require $235\times 20 = 4700$ evaluations to achieve the same precision.  This is what a decent choice of distribution can accomplish.

Code to reproduce the example
#
# The percentage point function for a Pareto distribution.
#
h <- function(q, pow) (1-q) ^ (1 / (1 - pow)) - 1

xmax <- 52 * log(2) # An effectively infinite (log) upper limit.
pow <- 3
#
# Study two estimation procedures by applying them to many independent samples.
# (Takes less one second for 50,000 iterations.)
#
# set.seed(17) # Use for reproducibility
z <- replicate(5e4, {
  q <- runif(1e1)   # Random uniform values
  #
  # The Pareto sampling estimate.
  #
  x <- h(q, pow)    # A random draw
  x. <- h(1-q, pow) # Its antithetic variable
  y <- exp(-x) * (1 + x)^pow / (pow - 1)
  y. <- exp(-x.) * (1 + x.)^pow / (pow - 1)
  z1 <- mean(c(y, y.))
  #
  # The Uniform sampling estimate.
  #
  y <- exp(-q*xmax)
  y. <- exp(-(1-q)*xmax)
  z2 <- xmax * mean(c(y, y.))
  c(Pareto=z1, Uniform=z2)
})
#
# Plot the sampling distributions of the two estimation procedures.
#
par(mfrow = c(1,2))
hist(z[1,], breaks=50, freq=FALSE, main="Pareto Solution", xlab="Z")
hist(z[2,], breaks=25, freq=FALSE, main="Uniform Solution", xlab="Z")
par(mfrow = c(1,1))
#
# How much more precise is the first procedure compared to the second?
#
(var(z[2,]) / var(z[1,]))

