Цикл списка и использование элементов в этом списке

Question

В приведенном ниже коде вместо того, чтобы начинать с p = 0,01 и затем увеличивать его, я хотел бы иметь возможность сделать что-то вроде p = (1:99) / 100, а затем просто перебрать p.Например, вместо p = 0,01, давайте пусть p = (1:99) / 100.Теперь я пытаюсь заменить 1:99 в цикле for на p.Однако, когда я запускаю код, у меня возникают проблемы с покрытием [i] (возвращается число (0)).Кажется, что это должно быть довольно тривиально, поэтому я надеюсь, что просто пропускаю что-то в своем коде.

Кроме того, если вы видите легкое повышение эффективности, пожалуйста, не стесняйтесь вмешиваться!Спасибо =)

w=seq(.01,.99,by=.01)
coverage=seq(1:99)
p=0.01

for (i in 1:99){
    count = 0
    for (j in 1:1000){
        x = rbinom(30,1,p)  
        se=sqrt(sum(x)/30*(1-sum(x)/30)/30)
        if( sum(x)/30-1.644854*se < p && p < sum(x)/30+1.644854*se )
        count = count + 1
    }
    coverage[i]=count/1000
    print(coverage[i])
    p=p+0.01 
}

Answer 1

Я бы работал над средней частью, а не над этим внешним контуром.

coverage <- p <- 1:99/100
z <- qnorm(0.95)

for (i in seq(along=p) ){
  # simulate all of the binomials/30 at once
  x <- rbinom(1000, 30, p[i])/30

  # ses
  se <- sqrt(x * (1-x)/30)

  # lower and upper limits
  lo <- x - z*se
  hi <- x + z*se

  # coverage
  coverage[i] <- mean(p[i] > lo & p[i] < hi)
}

Это почти мгновенно для меня.Хитрость заключается в том, чтобы векторизовать эти симуляции и вычисления.Увеличение до 100 000 повторов моделирования заняло всего 4 секунды на моем 6-летнем Mac Pro.

(Вы можете увеличить количество повторений, чтобы увидеть структуру в результатах; plot(p, coverage, las=1) с повторениями по 100 тыс. Даетследующее: это не было бы ясно только с 1000 повторениями.)

Answer 2

Чтобы ответить на исходный вопрос, установка i in p, где p - 0,01, 0,02 и т. Д., Означает, что coverage[i] пытается выполнить coverage[0.01];поскольку [] требует целое число, оно обрезает его до нуля, в результате чего получается цифра нулевой длины, numeric(0).

Одно из других решений лучше, но для справки, сделать это, используя исходный циклвы, вероятно, захотите что-то вроде

p <- seq(.01, .99, by=.01)
coverage <- numeric(length(p))
N <- 1000
n <- 30
k <- qnorm(1-0.1/2)
for (i in seq_along(p)) {
    count <- 0
    for (j in 1:N) {
        x <- rbinom(n, 1, p[i]) 
        phat <- mean(x) 
        se <- sqrt(phat * (1 - phat) / n)
        if( phat-k*se < p[i] && p[i] < phat+k*se ) {
            count <- count + 1
        }
    }
    coverage[i] <- count/N
}
coverage

Answer 3

@ У Карла Бромана есть отличное решение, которое действительно показывает, как векторизация имеет значение.Тем не менее, его можно немного улучшить (около 30%):

Я предпочитаю использовать vapply - хотя улучшение скорости здесь не заметно, поскольку цикл составляет всего 99 раз.

f <- function(n=1000) {
    z <- qnorm(0.95)/sqrt(30)

    vapply(1:99/100, function(p) {
      # simulate all of the binomials/30 at once
      x <- rbinom(n, 30, p)/30

      zse <- z*sqrt(x * (1-x))

      # coverage
      mean(x - zse < p & p < x + zse)
    }, numeric(1))
}

system.time( x <- f(50000) ) # 1.17 seconds

Вот версия исходного кода OP с использованием vapply, и она примерно на 20% быстрее, чем оригинал, но, конечно, все же на несколько медленнее, чем полностью векторизованные решения ...

g <- function(n=1000) {
    vapply(seq(.01,.99,by=.01), function(p) {
        count = 0L
        for (j in seq_len(n)) {
            x <- sum(rbinom(30,1,p))/30 

        # the constant is qnorm(0.95)/sqrt(30)
            zse <- 0.30030781175850279*sqrt(x*(1-x))
            if( x-zse < p && p < x+zse )
                count = count + 1L
        }
        count/n
    }, numeric(1))
}

system.time( x <- g(1000) ) # 1.04 seconds