Using control variates & antithetic method with Monte Carlo

Question

Supposing $g(x)=\sqrt[3]{x}$, I want to calculate the expected value of g, $E(\sqrt[3]{x})$, using Monte Carlo method, by generating $x_i$ from a Weibull distribution with parameters $(1,5)$.

After that, I want to use the control variates method and the antithetic method in order to to reduce the variance of my estimator, which I found with the simple Monte Carlo. And here is my problem, I do not know how to do these methods.

I would appreciate if someone could help me do that or give any tip/help.Thank you very much for your concern, in advance.

What I have done so far

Supposing $S$ is our estimator, then we know that $S=(\sum \limits_{i=1}^{N} g(x_i))/N$.

Using Matlab, I found the expected value $S$ by generating 1000 random numbers from the weibull(1,5) distribution and calculate the sum. Here is my algorithm:

N=1000
sum=0;
for i=1:N;
  X = wblrnd(1,5);  
  res(i)=X.^(1/3); 
  sum=sum+res(i);
end
S=sum/N

Answer 1

There is no one way to implement either control or antithetic variates, however, a couple of examples may help.

Antithetic variables: Imagine instead of the random number generator you actually used, you generated a $U(0,1)$ variate, call it $u$, and ran it through the inverse CDF of the Weibull(1,5) distribution, thereby generating a Weibull(1,5) variate. For the next random number, use $1-u$ instead of generating a new $u$. For subsequent random numbers, alternate generating a new $u$ and using $1-u$. This helps to "balance" high and low values from your random number stream, thus reducing variability of your final estimates.

Control variates: These are "extra" variables that are correlated with the result, enabling you to do something like run a regression on your results against the control variates to get a more accurate estimate. In your case, for example, you know the true mean of the Weibull dist'n (5), so you could use the $x_i$ as a control variate. You would calculate the improved estimate:

$S^* = S - \frac{\widehat{cov}(x,res)}{\widehat{var}(x)} * (\bar{x} - 5)$

where the covariance and variance terms are estimated from the data. This helps correct the estimate for random number streams that are not, in some relevant way, totally representative of the underlying distribution.

Both methods, esp. control variates, are more general than these two examples might lead you to believe. The wikipedia links are at best rough introductions; plenty of books and papers covering both techniques are out there if you want to go more into depth.

Answer 2

As correctly explained by jbowman, you have to create a negative correlation between your simulations to implement antithetic variables: with a generic wblrndfunction this is not possible. You thus have to get back to the definition of a Weibull $\mathcal{W}(\lambda,k)$, which is a scale transform of the $k$th power of an exponential variate, i.e. $$ X\sim \mathcal{W}(\lambda,k) $$ is equivalent to $$ (X/\lambda)^k \sim \mathcal{E}(1) $$ This can be reworded in terms of simulation as $$ X=\lambda(-\log U)^{1/k}\,,\qquad U\sim\mathcal{U}(0,1)\,. $$ Therefore, you can implement the antithetic method by using a sample of uniforms, $U_1,\ldots,U_n$ and its complent $1-U_1,\ldots,1-U_n$ and compare the variance of the estimator of $\mathbb{E}[X^{1/3}]$ with the corresponding estimator based on $U_1,\ldots,U_{2n}$. To exhibit the improvement, you have to run a Monte Carlo experiment repeating the computation of those variances on many samples of size $n$. (You cannot see the impact of antithetic simulation on a single run.)

The control variate is implemented in your case by choosing a known moment of the Weibull, for instance as suggested by jbowman, $$ \mathbb{E}[X] = \lambda \Gamma(1+1/k) $$ and using $X$ as the control variate. This means you compute the average of the simulated $X_i$'s along the average of the $X_i^{1/3}$ and the empirical covariance between the $X_i$'s and the $X_i^{1/3}$ as well as the empirical variance of the $X_i$'s to use jbowman formula. Again, checking the improvement brought by the control variate requires a Monte Carlo experiment with several runs.