[R] Gamma Simulation from Exponential Distribution: find sample size given power and sig.Level

Black Yellow bl@cki@h952 @ending from gm@il@com
Sun Oct 28 21:26:27 CET 2018


I built my function to simulate two gamma distributions X and Y based on the
sum of i.i.d exponential distributions. Assume my code is correct about this
simulation, I am interested in finding an equal sample size n for X and Y
such that n can be determined given 90% power and 5% significance level.

My thought process is that I create a sequence of n's value and use a for
loop to reiterate each n value in the sequence along with an "if" statement
saying that if given the n value, 95% of the time when I replicate the
process, the probability that I reject the null hypothesis given the
alternative is true is 90%.

The problem I am facing is that n value is always the maximum within my
n_seq. I am doubting there must be something wrong when I generate my two
Gamma distributions.

Any input is appreciated.

    set.seed(3701)
    n_seq = seq(from = 1, to = 200, by = 1)
    Z.list = numeric(reps)
    var.list = numeric(reps)

    gamma.sim = function(n, mu){
      for(j in 1:5e5){
          u.list = runif(n = n)
          x.list = -mu*log(1-u.list)
          Z.list[j] = sum(x.list)
      }
      return(Z.list)
    }


    # Generate n generations for X.gamma and Y.gamma
    gen_p_vals <- function(n = reps){
      p_vec <- numeric(reps)
      for(ii in 1:reps){
        x <- gamma.sim(n = 79, mu = 1)
        y <- gamma.sim(n = 80, mu = 1)
        t_stat_numerator <- mean(x) - mean(y) - 0 # Test Statistic under
null hypothesis
        s_p <- sqrt(sum((x-mean(x))^2) + sum((y-mean(y))^2)) /
(sqrt(2*reps-2))
        t_stat <- t_stat_numerator /(s_p*(sqrt((1/reps)+(1/reps))))
        p_vec[ii] <- 2* pt(-abs(t_stat), 2*reps-2)
      }
      return(p_vec)
    }
    power_pvals <- gen_p_vals(n = 10)
    mean(power_pvals < 0.05)

    for(j in 2:100){
      power_pvals = gen_p_vals(n = n_seq[j])
      if (mean(power_pvals < 0.05) == 0.90){
        print(j)
      }
    }

	[[alternative HTML version deleted]]



More information about the R-help mailing list