Как написать цикл for, когда функция увеличивается с каждой итерацией?

Я пытаюсь оценить вероятность обнаружения животных из n.sites за несколько периодов наблюдения, когда животные удаляются, а обнаружение меняется во времени и пространстве. Это работает, если я делаю что-то подобное для 5 периодов наблюдения:

for(i in 1:nsites){
  mu[i,1] <- p[i,1]
  mu[i,2] <- p[i,2]*(1-p[i,1])
  mu[i,3] <- p[i,3]*(1-p[i,1])*(1-p[i,2])
  mu[i,4] <- p[i,4]*(1-p[i,1])*(1-p[i,2])*(1-p[i,3])
  mu[i,5] <- p[i,5]*(1-p[i,1])*(1-p[i,2])*(1-p[i,3])*(1-p[i,4])
}

Вероятность в момент времени 2 зависит от вероятности в момент времени 1, а вероятность в момент времени 3 зависит от вероятностей в моменты времени 1 и 2. Если бы я делал это только для 5 периодов времени, было бы несложно написать это вне. Но поскольку я получаю 10, 15, 20+ периодов времени, записывать их довольно запутанно. Я чувствую, что должен быть способ написать этот цикл без ввода каждого шага, но я просто не могу придумать, как это сделать. Может быть, дополнительная индексация или другой оператор управления или функция мощности. Если бы p[i] было одинаковым в каждом j-м наблюдении (т.е. p[i,1] = p[i,2] = p[i,3] и т. д.), то это было бы:

p[i]*(1-p[i])^5

Любые предложения будут ценны.

Это код языка BUGS. Я работаю в R и отправил код в JAGS через пакет rjags. ОШИБКИ, R или псевдокод подойдут для моих целей.

Вот код R, который имитирует проблему:

set.seed(123)
testp <- matrix(runif(108, 0.1, 0.5), 108, 5)
testmu <- matrix(NA, 108, 5)

for(i in 1:nsites){
  testmu[i,1] <- testp[i,1]
  testmu[i,2] <- testp[i,2]*(1-testp[i,1])
  testmu[i,3] <- testp[i,3]*(1-testp[i,1])*(1-testp[i,2])
  testmu[i,4] <- testp[i,4]*(1-testp[i,1])*(1-testp[i,2])*(1-testp[i,3])
  testmu[i,5] <- testp[i,5]*(1-testp[i,1])*(1-testp[i,2])*(1-testp[i,3])*(1-testp[i,4])
}

Спасибо за любую помощь. Дэн


r for-loop jags winbugs
person djhocking    schedule 21.11.2013    source источник

Ответы (3)


arrow_upward
3
arrow_downward

Ответ @Frank чище (и, вероятно, быстрее), но это также сработает и может быть немного проще для понимания.

testmu2 <- matrix(NA, 108, 5)
nsites = 108
np = 5

for (i in 1:nsites) {
  fac <- 1
  testmu2[i,1] <- testp[i,1]
  for (j in 2:np) {
    fac <- fac * (1-testp[i,j-1])
    testmu2[i,j] <- testp[i,j] * fac
  }
}
max(abs(testmu2-testmu))
[1] 2.775558e-17
person jlhoward    schedule 22.11.2013
comment
+1. Я не знаю, легче ли это читать, но я мог бы написать тело внутреннего цикла как testmu2[i,j] <- testp[i,j]*(fac <- fac * (1-testp[i,j-1])) - person Frank; 22.11.2013
comment
@Frank - я даже не знал, что можно использовать два задания в одной строке. Спасибо #многоузнать - person djhocking; 22.11.2013
comment
Вы были правы в своих оценках скорости. Метод Фрэнка был в десять раз быстрее, чем ваш или мой. - person IRTFM; 22.11.2013
comment
Я выбрал это, потому что оно более общее (не использует специфические функции R). У него все еще может быть проблема с языком BUGS, потому что fac переназначается, но я надеюсь, что смогу обойти это. - person djhocking; 23.11.2013

arrow_upward
4
arrow_downward

Это действительно выглядит как задача, хорошо подходящая для R Reduce:

testmu3 <- matrix(NA, 108, 5)
nsites = 108
np = 5

for (i in 1:nsites) {
   testmu3[ i, ] <- Reduce( function(x,y) x*(1-y), testp[i, ], 
                                                   accumulate=TRUE)
}
max(abs(testmu3-testmu))
[1] 0

Параметр накопления создает растущий вектор промежуточных результатов.

> testp[1, ]
[1] 0.215031 0.215031 0.215031 0.215031 0.215031

> Reduce( function(x,y) x*(1-y), testp[1, ],  accumulate=TRUE)
[1] 0.21503101 0.16879267 0.13249701 0.10400605 0.08164152
person IRTFM    schedule 22.11.2013
comment
+1, Круто. Я не знал об этой опции accumulate, хотя я довольно часто использую Reduce. Для справки OP... цикл можно было бы устранить и здесь, используя t(apply(testp,1,function(z) Reduce( function(x,y) x*(1-y),z,accumulate=TRUE))) - person Frank; 22.11.2013
comment
На самом деле я никогда раньше не слышал о функции Reduce. Очень аккуратный. Спасибо за ответ, а также за новую функцию! - person djhocking; 22.11.2013

arrow_upward
2
arrow_downward

Вот один из способов:

testmu2 <- testp*t(apply(cbind(1,1-testp[,-5]),1,cumprod))

На моем компьютере они почти совпадают:

> max(abs(testmu2-testmu))
[1] 2.775558e-17

Я не знаю о BUGS/JAGS, но идея состоит в том, чтобы сначала взять кумулятивный продукт вашей 1-p матрицы по ее столбцам, а затем взять p * результат.

person Frank    schedule 21.11.2013
comment
Очень ловко, спасибо. У меня все еще есть проблемы с обдумыванием функции применения. Циклы for для меня гораздо более естественны, но функция применения намного быстрее. - person djhocking; 22.11.2013