3

我试图估计当动物被移除并且检测时间和空间发生变化时，在多个观察期内从 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] 在每个 jth 观察中都相同（即 p[i,1] = p[i,2] = p[i,3] 等），它将是：

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

任何建议将不胜感激。

这是 BUGS 语言代码。我在 R 中工作并通过 rjags 包将代码发送到 JAGS。BUGS、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])
}

谢谢你的帮助。担

4

3 回答 3

4

这确实看起来像一个非常适合 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
于 2013-11-22T00:30:44.907 回答
3

@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
于 2013-11-22T00:00:05.373 回答
2

这是一种方法：

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

在我的电脑上，它们几乎匹配：

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

我不知道 BUGS/JAGS，但我的想法是先将 1-p 矩阵的累积乘积穿过其列，然后再取 p*result。

于 2013-11-21T23:40:35.477 回答