Search code examples
rfor-loopwinbugsjags

How to write for loop when function increases with each iteration?


I am trying to estimate the probability of detecting animals from n.sites over multiple observation periods when animals are removed and detection changes in time and space. It works if I do something like this for 5 observation periods:

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

The probability at time 2 is dependent on the probability at time 1 and the probability at time 3 is dependent on the probabilities at times 1 and 2. If I were only doing this for 5 time periods it wouldn't be a big deal to write this out. But as I get 10, 15, 20+ time periods, it's is quite messy to write out. I feel like there should be a way to write this loop without typing out each step, but I just can't think of how to do it. Maybe additional indexing or other control statement or power function. If p[i] were the same in each jth observation (i.e. p[i,1] = p[i,2] = p[i,3], etc.) it would be:

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

Any suggestions would be greatly appreciated.

This is BUGS language code. I work in R and sent the code to JAGS via the rjags package. BUGS, R, or pseudo code would suit my purposes.

Here is R code that would simulate the problem:

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

Thanks for any help. Dan


Solution

  • @Frank's answer is cleaner (and faster, probably), but this will also work and might be a little easier to understand.

    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