Using R, I am trying to simulate how the power of a Monte Carlo two-sample test for central tendency changes with sample size. However, my simulation results does not show power increasing with sample size which is clearly wrong.

Could someone please advise where I am going wrong? Below is the basic set-up.

The null hypothesis is:


The alternative hypothesis is:


The test statistic is:


For my simulation, I use two samples $X\sim N(\mu_X, \sigma_y)$, $Y\sim N(\mu_Y, \sigma_Y)$ of size $n_X$ and $n_Y$ respectively. I run 100 simulations for every increase in the sample size of X.

Below is my coded simulation with comments:

 mc.pvalue<-function(nx, ny, mu.x, mu.y, sig.x, sig.y){
    x<-rnorm(nx, mu.y, sig.x)            # Generate samples under H1: mu_x > mu_y
    y<-rnorm(ny, mu.x, sig.y)            
    obs.diff<-mean(x)-mean(y)            # Test-stat is observed diff. in means 
    se.x<-sd(x)/sqrt((length(x)))        # Calculate standard errors 
    count=0                              # Set counter equal to zero
    for(i in 1:1000){                    # Generate 1000 pairs of samples and calculate
      x1<-rnorm(nx, mu.y, se.x)          # the difference in means between each pair.
      y1<-rnorm(ny, mu.y, se.y)          # The samples have been generated under
      sim.diff<-mean(x1)-mean(y1)        # H0: mu_x = mu_y
      if(sim.diff<= obs.diff){count=count+1} 
    count/1000                           # Estimated p-value for test statistic

# Calculate 100 estimated p-values for the test statistic. Then find the proportion of times
# reject hypothesis at level of significance 0.05

calc.pvalues<-function(nx, ny=10, mu.x=21, mu.y=20, sig.x=0.5, sig.y=0.25, 
  pvalues<-replicate(n.sim, mc.pvalue(nx, ny, sig.x, sig.y, mu.x, mu.y))

chge.sample.size<-seq(2, 100, by=0.05)          # Sample size x increasing from 1 to 100
result<-sapply(chge.sample.size, calc.pvalues)  # Apply the calc.pvalues function to
                                                # estimate p-value for each sample size
                                                # of x    
plot(chge.sample.size, result)                  # Plot increasing sample size of X 
                                                # against test power

When I an edit the code to calculate the power of a t-test, it works perfectly:

calc.pvalue.t<-function(nx, ny, sig.x, sig.y, mu.x, mu.y){    
  x<-rnorm(nx, mu.x, sig.x)
  y<-rnorm(ny, mu.y, sig.y)
  t.test(x, y, alternative="greater", 
         paired=FALSE, conf.level=0.95)$p.value

calc.t<-function(nx, ny=10, mu.x=21, mu.y=20, sig.x=0.5, sig.y=0.25, 
                    n.sim=100, alpha=0.05){
  p.values<-replicate(n.sim, calc.pvalue.t(nx, ny, sig.x, 
                                            sig.y, mu.x, mu.y))

chge.sample.size<-seq(2, 10, by=1)
tpower1<-sapply(chge.sample.size, calc.t)

plot(chge.sample.size, tpower1)

Something seems to be going wrong with the first part of my code when I try to do the Monte Carlo test

  • $\begingroup$ I thought that if you assumed H1 was true when generating the data, the above calculations give power. Assuming H0 is true gives size? $\endgroup$ Commented Nov 14, 2014 at 20:10
  • 1
    $\begingroup$ I don't think your mc.pvalue function does anything meaningful. (By the way, the variable cnt is not defined -- you mean count, right?) $\endgroup$
    – Adrian
    Commented Nov 14, 2014 at 20:12
  • $\begingroup$ As written, mc.pvalue(100, 100, 10, 50, 1, 1) will always return 1. What does that mean? $\endgroup$
    – Adrian
    Commented Nov 14, 2014 at 20:13
  • $\begingroup$ Yes, sorry I have changed it to count. $\endgroup$ Commented Nov 14, 2014 at 20:15
  • $\begingroup$ That the probability of getting a test statistic at least as large as the one observed is 1. $\endgroup$ Commented Nov 14, 2014 at 20:17

1 Answer 1


I think maybe the reason is that the empirical distribution is changed with every observation in the function "mc.pvalue". So the 100 p-values are calculated with different empirical distribution.

  • $\begingroup$ Could you expand on this and perhaps provide an alternative code. As it stands this is not a complete answer. Please see these guidelines for a good answer: stats.stackexchange.com/help/how-to-answer $\endgroup$
    – André.B
    Commented Oct 17, 2019 at 2:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.