\[f(x)=\frac{x}{\sigma^2}e^{-x^2/(2\sigma^2)}\] \(x\geq0\), \(\sigma \geq 0\)
Implement a function to generate samples from a Rayleigh(\(\sigma\)) distribution and compare the variance reduction between \(\frac{X_1+X_2}{2}\) and \(\frac{X+X'}{2}\):
Let \(Y=U[0,1]\) then to use the antithetic approach:
\(Y_j=g(F_X^{-1}(U_1^{(j)},...,F_X^{-1}(U_n^{(j)}))\)
\(Y'_j=g(F_X^{-1}(1-U_1^{(j)}),...,F_X^{-1}(1-U_n^{(j)}))\)
for \(U_i^{(j)}\) uniform for \(i=1,...,n\), \(j=1,...,\frac{m}{2}\)
Then \(\hat{\theta} = \frac{2}{m}\sum_{j=1}^{m/2}\Big(\frac{Y_j+Y'_j}{2}\Big)\)
Given that the CDF of the Rayleigh distribution is:
\[F(x)=1-e^{-x^2/(2\sigma^2)}\] for the scale \(\sigma >0\). Setting this to \(U\), the uniform distribution between \([0,1]\), I get:
\[X_{Rayleigh}=\sigma\sqrt{-2\ln{U}}\]
Where \(X\) has a Rayleigh distribution with paramter \(\sigma\) per this Wikipedia page.
First I create a function to generate a \(X_{Rayleigh}\) across some vector of Uniform RVs:
my_raylei <- function(x, sig){
return(sig * sqrt(-2 * log(x)) )
}Next I create two datasets on which to apply this function. One is going to be a list of uniform RV’s representing \(X_1+X_2\) and one is going to be a list of \(X + X'\) where the half the \(X_i\) are \(U\) and the other half are \(X'_i\) are \(1-U\) and then I’m going to apply the above function across all of the values and compare the variances:
set.seed(1060)
antithetic_compare <- function(sig, n){
X1 <- runif(n/2)
X2 <- runif(n/2)
X <- X1
X_anti <- 1-X1
partA <- append(X1, X2)
partB <- append(X, X_anti)
A <- var(my_raylei(partA, sig))
B <- var(my_raylei(partB, sig))
# return the variance for each and
# the % difference between them
return(data.frame(`X1+X2`=A,
`X+Xanti`=B,
`pdiff` =round((B-A)/A * 100,4)))
}
x <- antithetic_compare(3,100)
x## X1.X2 X.Xanti pdiff
## 1 4.15897 4.358541 4.7986Above, we see a rougly 4.8% improvement in variance. That improvement required, in theory, half the number of calculations which is significant.