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]] ______________________________________________ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.