Ones (zeros) trick

658 views
Skip to first unread message

Rahel Sollmann

unread,
Jan 14, 2018, 12:57:47 PM1/14/18
to nimble-users
Dear all,
I am new to Nimble and I am currently trying to implement an SCR model with an irregular state space S. I am trying to do that using the Ones Trick (used a bunch in BUGS/JAGS and described for the purpose of constraining S by Mike Meredith: http://www.mikemeredith.net/blog/1309_SECR_in_JAGS_patchy_habitat.htm). 
Below is the model code; "habmat" is a matrix indicating whether a cell in S is suitable habitat (1) or not (0), and it's passed to NIMBLE as a constant. I can't get this to work (the SCR model with rectangular S works fine for the same data) and I was wondering if there are any restrictions regarding use of the Ones Trick in Nimble. 
I know I could re-code this and have the activity centers, s, be drawn from a multinomial distribution (rather than a Uniform, and then truncate to an integer); at least in JAGS, though, that would slow the model down considerably and I am dealing with a large data set, so trying to avoid making things slower. But maybe that's less of an issue in NIMBLE?
Any thoughts would be much appreciated. Happy to provide the full example but figured I'd check first whether there might be a general issue with the Ones Trick.
Thanks!
Rahel

PS: This is a message I get when I execute nimbleModel(): ‘habmat’ is not a valid field or method name for reference class “bird2_modelClass_UID_39_UID_40”. 
       Also, I do provide initial values for "OK" (see code below) and I generate starting values for s that are all within the irregular S.


birdcode2<-nimbleCode({
lam0~dunif(0,5)
sigma~dunif(0, 100)
psi~dunif(0,1)

for(i in 1:M){
 s[i,1]~dunif(1,Xu)
 s[i,2]~dunif(1,Yu)
 ind1[i]<-trunc(s[i,1]+0.5)
 ind2[i]<-trunc(s[i,2]+0.5)
 pOK[i] <- habmat[ind1[i], ind2[i]] # habitat check --> pOK[i]=1 if s falls within suitable habitat 
 OK[i] ~ dbern(pOK[i]) # OK[i] = 1 if s falls within suitable habitat, 0 if not, the ones trick

for(j in 1:J){
d[i,j]<- pow(pow(s[i,1]-X[j,1],2) + pow(s[i,2]-X[j,2],2),0.5)
y[i,j] ~ dpois(p[i,j])
p[i,j]<- z[i]*lam0*exp(-(d[i,j]*d[i,j])/(2*sigma*sigma))
}
}
N<-sum(z[1:M])

})
 

Jose Jimenez

unread,
Jan 14, 2018, 1:25:27 PM1/14/18
to nimble-users
Dear Rahel,

I did no use this aproach, but if you want to try the zero's trick from Chandler (https://groups.google.com/forum/#!searchin/spatialcapturerecapture/zero$20trick%7Csort:date/spatialcapturerecapture/NzqUovn8jF0/4Z4OzGRAAgAJ) it works very well in Nimble to reduce the computation time.
You have to supply buffer as a data, not a constant

Regards,
José

#####################################################################
# IN JAGS:
.............
for (i in 1:M) {
  z[i] ~ dbern(psi)
  s[i,1] ~ dunif(xlim[1], xlim[2])
  s[i,2] ~ dunif(ylim[1], ylim[2])

  for(j in 1:J) {
    d2[i,j] <- (X[j,1] - s[i,1])^2 + (X[j,2] - s[i,2])^2
    outj[i,j] <- sqrt(d2[i,j]) > buffer # is d(i,j) > buffer??
    lambda[i,j] <- lambda0*exp(-d2[i,j]/(2*sigma^2))
    # No need to loop over occs
    lzk[i,j] <- lambda[i,j]*z[i]*K[j]
  }
  out[i] <- equals(sum(outj[i,]), J) # Is s(i) outside buffer?
  zeros[i] ~ dbern(out[i])
}
................

#####################################################################
# IN NIMBLE:
................
for (i in 1:M) {
  z[i] ~ dbern(psi)
  s[i, 1] ~ dunif(xlim[1], xlim[2])
  s[i, 2] ~ dunif(ylim[1], ylim[2])
  d2[i,1:J] <- (s[i,1] - X[1:J,1])^2 + (s[i,2] - X[1:J,2])^2
  p[i,1:J] <- z[i] * p0[1:J] * exp(-d2[i,1:J]/(2*sigma^2))

  for (j in 1:J) {
    outj[i,j] <- sqrt(d2[i,j]) > buffer
    y[i,j] ~ dpois(p[i,j]*K[j])
  }
  out[i] <- sum(outj[i,1:J])==J    # Is s(i) outside buffer?
  zeros[i] ~ dbern(out[i])
}
................

Rahel Sollmann

unread,
Jan 14, 2018, 1:32:18 PM1/14/18
to nimble-users
Thanks so much Jose, I will try this and report back. 
Best
Rahel

Perry de Valpine

unread,
Jan 14, 2018, 2:19:39 PM1/14/18
to Rahel Sollmann, nimble-users
Hi Rahel,

I played with your code a bit.  There are a couple of issues that are on the list of things to improve, but let me give you some work-arounds.

1. I think you need to provide habmat as data, even though it would make sense as constants.

2. I think you also need to provide a habmat entry in the dimensions argument (e.g. "dimensions = list(habmat = c(10, 10))", if habmat will be 10x10), even though you'd think the dimensions would be determined by its value as data. (The problem is that in nimble, data can be changed.)

-Perry


--
You received this message because you are subscribed to the Google Groups "nimble-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users+unsubscribe@googlegroups.com.
To post to this group, send email to nimble...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/bb325b81-c047-4385-bfef-05afdc02bcdd%40googlegroups.com.

For more options, visit https://groups.google.com/d/optout.

Rahel Sollmann

unread,
Jan 14, 2018, 2:42:41 PM1/14/18
to nimble-users
Based on a shirt test run this appears to be working. Thanks again!


On Sunday, January 14, 2018 at 10:25:27 AM UTC-8, Jose Jimenez wrote:

Rahel Sollmann

unread,
Jan 14, 2018, 2:44:53 PM1/14/18
to Perry de Valpine, nimble-users
Thanks Perry;
Jose provided an alternative version of the ones trick and that works fine. I actually think I figured out what I did wrong in the original version, which includes not providing habmat as data.
Best
Rahel

Perry de Valpine

unread,
Jan 14, 2018, 2:59:30 PM1/14/18
to Rahel Sollmann, nimble-users
Great!
Reply all
Reply to author
Forward
0 new messages