[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