Let $X$ be a random variable uniformly distributed on the interval $[0.1]$. It has density function $1$ on the interval, and $0$ elsewhere. Let $Y=e^{-X}$. We first compute the mean and variance of the random variable $Y$.
We have $E(Y)=\int_0^1 e^{-x}\,dx=1-e^{-1}$. Call this $\mu$. It is the exact value of our integral. This is the whole point of the simulation: If we take a sort of big sample, such as $1000$, the sample mean will probably be close to the true mean, which is the true value of the integral.
To find the variance of $Y$, we use the shortcut formula $\text{Var}{Y}=E(Y^2)-(E(Y))^2$. We have $E(Y^2)=E(e^{-2X})=\int_0^1 e^{-2x}\,dx=\frac{1}{2}(1-e^{-2})$. Now we can calculate the variance of $Y$, either as an exact expression, or by giving a decimal approximation. Call it $\sigma^2$.
In your program, you have $1000$ random variables $Y_1,Y_2,\dots, Y_{1000}$, each with the same distribution as the above $Y$. We are studying the random variable
$$\bar{Y}=\frac{Y_1+Y_2+\cdots +Y_{1000}}{1000}.$$
This random variable $\bar{Y}$ also has mean $\mu$. It can be shown (this is a general phenomenon) that the variance of $\bar{Y}$ is $\frac{\sigma^2}{1000}$.
We know $\sigma^2$, so we know the variance of $\bar{Y}$.
Note that $\bar{Y}$ is the result of a single run of your program, that is, $1000$ points.
We can imagine repeating the simulation $1000$ times, that is, effectively, using $1$ million points.
The theory is much the same. The results of the simulations are random variables $\bar{Y}_1,\dots, \bar{Y}_{1000}$. If we average them, the result has variance equal to the variance of $\bar{Y}$, divided by $1000$. so it is $\frac{\sigma^2}{1000^2}$.
It appears like you are trying to do something like this (see slide 15). In that case, what the monte carlo estimator $\langle I \rangle$ is doing does not make sense if you think of it as calculating the $E[f(x)]$. What that estimator is doing is calculating $\int_{-\infty}^{\infty} f(x) dx$ using a Monte Carlo method that samples the $x_i$ not uniformly but according to the density function $p(x_i)$, where $p(x_i)$ is chosen to sample the places where $f(x_i)$ is large more often than where it is small, thereby making better use of limited computer resources - this is called importance sampling (as copper.hat pointed out). To calculate the expected value of $f(x_i)$ you first need to specify the distribution of the $x_i$ and then generate Monte Carlo values of $f(x_i)$ with x's drawn from that distribution.
To summarize, the Monte Carlo estimator for the average of a random function and for the integral of a deterministic function are completely different things. You can't necessarily mix them. I think your confusion is thinking $\langle I \rangle$ is for $E[f(x_i)]$, which is it not. Your intuition on how to calculate $E[f(x_i)]$ is correct.
Best Answer
The wikipedia pages for the control variates method and the antithetic method are a pretty good start to learn about them.
The antithetic method seems the easiest to implement. What you should do is sample random numbers from your distribution (here the Weibull distribution) and for each of those numbers create another number according to a procedure that makes the new numbers be realizations of a random variable with negative covariance with respect to the original Weibull variable that is used for generating the first sequence of numbers. This way, you reduce the total variance on your estimate of the expected value.
Practically, I think the following should do the trick (I'm using R code):
Maybe it's more efficient to generate uniformly distributed random numbers $u_i$, take their complements $1-u_i$ and from these generate the Weibull distributed variables using
qweibull. Also this might not be the best choice to minimize the variance.