JAGS检测仿真误差/生态学家Bungle高级统计

Question

我对贝叶斯世界很陌生，很高兴能把它应用到我的研究中。作为第一步，我试图建立一个模拟模型，估计在10公里半径范围内的动物的平均数量和探测概率作为噪声水平的函数。基本上，我是问在两个地方是否有1k的动物，但噪音的特性是不同的，模型能告诉我“真”的意思(1k动物)在两个研究区域之间是否相同吗？希望是的。

这是一个三步走的过程，它将检测到的呼叫数(kk1)和检测概率(Pdet)与信噪比(信噪比，在一个日志尺度上)联系起来，最终与区域内动物的平均数量联系起来。

我模拟了我的动物位置，源水平，传输损耗(取决于动物位置)，信噪比和Pdet功能如下。

##################################################################
library(runjags)
#library(plotrix)
library(rjags)

# Initializations, thanks Matt
### Gamma distribution specification function, thanks Matt! ###
gammaPr<-function(mu, sd)
{
  shape<-mu^2/sd^2
  rate<-mu/sd^2
  x<-seq(max(0, mu-4*sd),mu+4*sd, length.out=100)
  #plot(x, dgamma(x, shape=shape, rate=rate), main=paste("Prior distribution with mean ",mu," and sd ", sd), type="l")
  #print(paste("The two parameters of the Gamma distribution are shape=",shape," and rate=",rate))
  return(c(shape,rate))
}

##################
# Simulated Data #
##################

# Number of animals within the detection radius (5km)
N=1000
#x=seq(0,5,by=1/N)



 # Idiot checking:
        # Since this is point sampling and (for the moment) we assume that the number
        # of animals available to be detected increases linearly with radius
    #x1 =    rbeta(N, shape1 = 2, shape2 = 1)
    #theta = runif(N, min = 0, max=2*pi)
    #radial.plot(lengths = x1, radial.pos = theta, rp.type = 's')

## We know the parameters for Source Levels values from the literature
nlpr=gammaPr(120, 30)

# Time series
tt=seq(1, 300, by=1)

# Number of animals at each site available to be sampled
N1=numeric(length=length(tt))
N2=N1
NL1=N1 # Noise Level Distribution
NL2=N1
TL1=N1 # Estimated Transmission loss
TL2=N1
kk1=N1 # Number of animals detected by the sensor 
kk2=N1

for (ii in 1:length(tt)){
  #N1[ii]=round(rnorm(1, mean = 700, sd=100)) # Screw it-there are enough problems to deal with
  #N2[ii]=round(rnorm(1, mean = 700, sd=100)) 
  N1[ii]=1000
  N2[ii]=1000

  # Random whale distances

  r1=rbeta(N1[ii], shape1 = 2, shape2 = 1)*5 # assumed (known)
  r2=rbeta(N2[ii], shape1 = 2, shape2 = 1)*5 


  # Random source levels for each whale
  slpr=gammaPr(131, 20)                     # assumed (known)
  SL1=rgamma(N1[ii], slpr[1], slpr[2])
  SL2=rgamma(N2[ii], slpr[1], slpr[2])

  # Noise Levels for each sensor
  nlpr=gammaPr(120, 30)
  NL1[ii]=rgamma(1, nlpr[1], nlpr[2])      # measured (known)

  nlpr=gammaPr(90, 30)
  NL2[ii]=rgamma(1, nlpr[1], nlpr[2])

  # TL Transmission loss from each animal to sensor
  TL1=20*log10(r1*1000)                   # this one is tricky, we can assume TL but it's based on multiple captures so how this make it into the probability density function escapes me 
  TL2=20*log10(r2*1000)                   

  # Signal to noise ratio for each animal call
  SNR1=SL1-(1/(1000*r1)^2)-NL1[ii]        # calculated
  SNR2=SL2-(1/(1000*r2)^2)-NL2[ii]

  # Detection probability for a given SNR threshold
  Pdet1=pbeta(SNR1, shape1  = 1.5, shape2 = 2.5)  # unknown, we will need to estimate shape1 and shape from the JAGS model
  Pdet2=pbeta(SNR2, shape1  = 1.2, shape2 = 2.5)
  ## Idiot checking again
    #yy=pbeta(seq(-40,40, by=.01), 1.5, 2.5)
    #plot(seq(-40,40, by=.01),yy, xlim=c(-2,2))


  # Number of animals detected is binomial based on the detection function
  kk1[ii]=sum(rbinom(N1[ii], size=1, prob = Pdet1))
  kk2[ii]=sum(rbinom(N2[ii], size=1, prob =  Pdet2))
}

data=data.frame(cbind(tt, N1,N2, NL1, NL2,kk1, kk2))
##################################################################

这很好，但我没有成功地建立一个JAGS模型，它可以预测一段时间内动物的平均数量，更不用说在两个不同的研究区域了。这个模型的最新版本，在逻辑上对我来说是有意义的，它会抛出关于重新定义节点Pdet1的错误。但是，由于这些是我试图估计的值，这是相当令人困惑的。如能提供任何援助，将不胜感激。

也许这是定义临时可变资产的问题？E.g.？但这似乎与我的问题有些不同，因为我只处理一个循环。

提前谢谢你抽出时间。

#############################################################################
# Model Fitting#
#############################################################################

# Build the JAGS model to determine if the true number 
# of animals at each location is different


# We are trying to estimate the mean population size (mu)
# Number of trials (detection periods)
N=nrow(data)

# We guess that the median range of animals within the detection 
# area is 2500 m (will be updated later with propationg modelling)
# This is stupid. There needs to be stochasticity in both values.
MedTL=20*log10(2500)
MedSL=rgamma(1, slpr[1], slpr[2])


# We are trying to estimate the population mean and the 
# detection function parameters
modeltext <- '
model{

  #### liklihood ### 
  for(ii in 1:N){

    # Number of animals detected is distributed by a normal function the total available (mu)
    # and the probability that they are detected(Pdet)
    kk1[ii] ~ dbin(Pdet, mu)

    # Pdet is a beta distribution with parameters estimated from above
    Pdet<-pbeta(SNR[ii], betaPar1, betaPar2)

    # SNR (signal to noise ratio) is the Source Level-Noise Levl-Transmission Loss
    SNR[ii]<-MedSL-MedTL-NL1[ii]

    # SL and RL need stochasticity.....
    }

  ## Priors ###
    # Mu is dependent on how the animals are distributed in space. In this case, 1000. Because that's what I told it to do. 
    mu~dpois(lambda)

  # We know mu must be greater than 0 and poison
  lambda~dgamma(11, .01) # approximately what they should be based on mu 1000, sd 300

  # We will assume we know a bit about the detection function, we will be helpful
  betaPar1 ~ dnorm(1,.3)
  betaPar2 ~ dnorm(2,.3)
}
  #data# kk1, NL1, N, slpr, MedTL, MedSL
  #monitor# betaPar1, betaPar2, mu, pdet
'

# We could just run the model like this:
results1 <- run.jags(modeltext, n.chains=2)

Answer 1

在JAGS模型中，Pdet是一个标量概率

kk1[ii] ~ dbin(Pdet, mu)

但是，您将Pdet指定为一个beta概率分布函数：

Pdet<-pbeta(SNR[ii], betaPar1, betaPar2)

大概你的意思是，Pdet是从beta发行版中提取出来的。这会是

Pdet ~ pbeta(SNR[ii], betaPar1, betaPar2)

或者(与您模拟的数据更一致)，您的意思是kk1ii是1000次伯努利试验的总和，每个个体的检测概率都是Pdet。所以，你需要做的是让kk1ii是一堆个别的伯努利试验的总和。在JAGS中，所观察到的数据必须是某些分布的随机实现，因此不能编写，例如

for(ii in 1:300){
for(jj in 1:mu){
kk11[jj] ~ dbern(p[jj])
}
kk1[ii] <- sum(kk11[jj])
}

对您来说，幸运的是，JAGS已经预见到了您的困境，并为您的目的实现了特殊的可观察函数dsum()。查一下雅格斯手册。

总之，，，除非我误解了你的模拟，你真正需要的是用kk1ii在时间步骤中循环，但是在这个循环中，你需要在每个个体上循环，每个人都需要在每一个时间步骤中都有自己的模型检测概率。在这一点上，你可能会想问，“但是他们可能有足够的数据来模拟跨时间步骤的个体检测概率吗？”答案是，您真正需要建模的是跨年单个检测概率的分布。这将与距离传感器的单个距离的分布、系统的(潜在的时变)噪声特性等有关。

好消息是(以第一次近似)，，只要您可以在R中模拟数据，就可以在JAGS中建立相应的模型。您现在面临的挑战是确保模拟生成您希望它生成的数据，然后确保JAGS模型实际上与仿真相对应。

Postscript

为此，我想指出一个奇怪的区别，在模拟和模型之间，你可能没有打算。在模拟中，您可以修复N=1000。在模型中，你的先验是伽马-泊松混合物。你有理由这么做吗？一个原因是，如果真正感兴趣的数量是Poisson参数lambda，它对应于系统中动物的密度，从中提取N。在这种情况下，，，您也应该将它作为数据模拟的一部分：

N <- rpois(1000)

而且，您应该监视lambda (而不是mu)作为JAGS输出中的主要关注量。您还需要在这里考虑N是绘制一次，然后固定，还是在不同的时间步骤重复绘制。，否则，，您可能会考虑一个更直观的优先于mu。