Гамма-симуляция из экспоненциального распределения: найти размер выборки с учетом мощности и сигн. Уровня - PullRequest
Новая
       7

Гамма-симуляция из экспоненциального распределения: найти размер выборки с учетом мощности и сигн. Уровня

0 голосов
/ 28 октября 2018

Я построил свою функцию для моделирования двух гамма-распределений X и Y на основе суммы iid-экспоненциальных распределений.Предположим, что мой код верен в этой симуляции, я заинтересован в том, чтобы найти равный размер выборки n для X и Y, чтобы n можно было определить, учитывая мощность 90% и уровень значимости 5%.

Мой мыслительный процесс заключается в том, чтоЯ создаю последовательность значений n и использую цикл for, чтобы повторить каждое значение n в последовательности вместе с оператором «if», говорящим, что если дано значение n, в 95% случаев, когда я повторяю процесс, вероятность того, что яотвергнуть нулевую гипотезу, учитывая, что альтернатива истинна и составляет 90%.

Проблема, с которой я сталкиваюсь, заключается в том, что значение n всегда является максимальным значением в моем n_seq.Я сомневаюсь, что должно быть что-то не так, когда я генерирую свои два гамма-распределения.

Любой вклад приветствуется.

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)
  }
}
...