Control Variates, Monte Carlo integration Exercise: Calculate $P(N>2.5)$ where $N$~$N(0,1)$ through simple monte carlo integration, and then use control variables to reduce the variance of my estimator.
I did
#Simulate Normal Random Variables
 randnormal<- function(n){
        Z<-NULL
        for(i in 1:n){
            U1=runif(1)
            U2=runif(1)
            Z[i]=sqrt(-2*log(U1))*cos(2*U2*pi)
        }
    return(Z)}

Then I realized the simple monte carlo integration through
#Monte Carlo Integration
montecarlo<-function(n){
    Z=randnormal(n)
    sum(Z>2.5)/n}

My question is how can I choose my control variable in this case, there is a criteria for selecting a control variable?
 A: The idea behind "control variates" is to use the same Monte-Carlo simulation to estimate the known expectation of a function $g$ that is similar to the one, $f$, whose expectation you are trying to estimate.  When $g$ and $f$ are strongly correlated, you can use the errors in estimating $g$ to correct the estimate of $f$. 
Therefore there are two principal criteria governing your choice of control variate:


*

*It must be a function whose expectation you know (or can somehow estimate to a high accuracy).

*It must be correlated with the target function $f$.
In the present case, $f(z)$ is zero for all $z \lt 2.5$ and jumps to unity for all $z \gt 2.5$.  Any function that will be correlated with it should either tend to increase with increasing $z$ or it should to tend to decrease: that will take care of (2).
Finding a $g$ whose expectation you can easily compute is more of a challenge.  But in writing its formula as
$$\mathbb{E}(g) = \frac{1}{\sqrt{2\pi}}\int_\mathbb{R} g(z) \exp(-z^2/2) dz,$$
ask yourself whether there might be a relatively simple form of $g$ that enables this integral to be calculated.  Indeed there is: one of the simplest choices would be $g(z) = z$, for then the integral could easily be computed using the substitution $z^2/2 = y, z dz = dy$.
Although $g(z)=z$ does work, following criterion (1) we might elect to modify it so that it most closely emulates $f$.  I would suggest doing this by cutting $g$ off so that it's zero when $f$ is zero, whence
$$g(z) = z I(z \ge 2.5)$$
ought to do the trick.  Its expectation is easily found in closed form as
$$\mathbb{E}(g) = \mathbb{E}(z I(z \ge 2.5)) = \frac{1}{\sqrt{2\pi}} \exp(-2.5^2/2).$$
You will find that this reduces the estimation variance by a factor of almost $200$--which means that with this technique, you need only $1/200$ as many iterations.  That's a pretty good improvement!
