Monte Carlo – Enhancing Efficiency in R Simulations Using Monte Carlo Methods

Question

This question appears to be off-topic because it focuses on programming, debugging, or performing routine operations, or it asks about obtaining datasets. You could try the support links we maintain or the Open Data site instead.

---

If the question is actually a statistical topic disguised as a coding question, then OP should edit the question to clarify this. After the statistical content has been clarified, the question is eligible for reopening.

Closed 2 years ago.

Improve this question

An exercise displayed in the image below shows example of finding the efficiency of an estimator. I am trying to replicate this example in R using monte carlo. Y1,Y2,Y3 are random samples of normal distribution.
The steps I thought I can take is using R function would be rnorm() as this function generates random numbers from a normal distribution. I can hard code in the calculation for mu1 and mu2 using the formula given with the three random sample numbers. We can use the randomly generated numbers generated for Y1, Y2, and Y3 and input them into the mu1 and mu2 functions stated in Example 5.4.5. and then find the variance of the two mu functions. We can find the relative efficiency by divide mu1/mu2.

I started off with code with a rnorm of 3:

Y <- rnorm(3)
mu1 = 0.25*X[1] + 0.5*X[2] + 0.25*X[3]
mu2 = (1/3)*X[1] + (1/3)*X[2] + (1/3)*X[3]

I can further get random deviates with:

Y <- replicate(1000,rnorm(3))

Would my next step be trying to store each mu1,mu2 to create a 2×1000 matrix? And use that matrix to get the variance through variance function in R?

Answer 1

I am not shure, whether I understood the question correctly but maybe a wrong answer is the best basis for an improved question?

We want 1000 instances of $Y_1$, $Y_2$ and $Y_3$:

Y1 <- rnorm(1000)
Y2 <- rnorm(1000)
Y3 <- rnorm(1000)

Now we can compute the $\mu$ vectorwise:

mu1 <- 0.25*Y1 + 0.5*Y2 + 0.25*Y3
mu2 <- (1/3)*(Y1 + Y2 + Y3)

And then variances are

var(mu1)
var(mu2)
var(mu1)/var(mu2)   # most of the time, this is larger then 1.0
                    # see histogram constructed in the edit below

Edit: It's easy to compute this many times as in

hist(replicate(1000, 
               {
                 n <- 1000
                 Y1 <- rnorm(n)
                 Y2 <- rnorm(n)
                 Y3 <- rnorm(n)

                 mu1 <- 0.25*Y1 + 0.5*Y2 + 0.25*Y3
                 mu2 <- (1/3)*(Y1 + Y2 + Y3)

                 var(mu1)/var(mu2)
                 }),
     main ="", xlab = expression(mu[1]/mu[2]),
     breaks = 20)
abline(v = (3/8)/(3/9), lwd = 3, lty = 3)