Stochastic node issue

[Lu]

Hi All,

I am a Nimble user and I have a question on the stochastic nodes or dynamic indexing. Can anyone help me with that?

The data and model code I used are shown as follows:

### data simulation ###

set.seed(1007)

y <- c()

y[1:100] <- rnorm(100, mean = 1, sd = sqrt(1)) 

y[101:200] <- rnorm(100, mean = 4, sd = sqrt(1)) 

### model code ###

# self-defined function

VFun <- function(M, l, alpha) {

  if (M < 1) {M <- 1}

  V <- matrix(rbeta(M * l, 1, alpha), nrow = M, ncol = l)

  if (l == 1) {

    Vtem <- V[1:M, 1]

  } else {

    library(purrr)

    Vtem <- as.numeric(unlist(purrr::map(1:M, ~(V[.x,1]*prod(1-V[.x,2:l])))))

  }

  return(Vtem)

}

VFunNim <- nimbleRcall(function(M = double(0), l = double(0), alpha = double(0)){},

                       Rfun = "VFun", returnType = double(1))

# nimble model code

BPCode <- nimbleCode({

  for (i in 1:N) {

    y[i] ~ dnorm(mu_y[i], psi)

    mu_y[i] <- mu0 + res[i]

    res[i] <- sum(H[1:L,i])

  }

  for (l in 1:L) {

    Ci[l] ~ dpois(gamma)

    Vtem[1:Ci[l]] <- VFunNim(Ci[l], l, alpha)

    for (j in 1:Ci[l]) {

      Theta[l,j] ~ dunif(0, N)

    }

    for (i in 1:N) {

    H[l,i] <- inprod((Theta[l,1:Ci[l]] <= i), Vtem[1:Ci[l]])

    }

  }

 alpha ~ dgamma(1, 1)

  gamma ~ dunif(0, 100)

  mu0 ~ dnorm(0, 0.01)

  psi ~ dgamma(0.01, 0.01)

})

BPdata <- list(y = y)

BPconsts <- list(L = 20, N = length(BPdata$y))

BPinits <- list(psi = 100, gamma = 2, mu0 = 0, res = rnorm(BPconsts$N),

                  alpha= 1, Ci = rpois(BPconsts$L, 2), Vtem = runif(30,0,1), 

                  Theta = matrix(runif(20*30, 0, BPconsts$N), nrow = 20, ncol = 30),

                  H = matrix(runif(BPconsts$L * BPconsts$N), nrow = BPconsts$L, ncol = BPconsts$N))

BPmodel <- nimbleModel(code = BPCode, name = "BPmodel", constants = BPconsts,

                         data = BPdata, inits = BPinits)

## error message:

Error in getSymbolicParentNodesRecurse(x, constNames, indexNames, nimbleFunctionNames,  : 

  R function ':' has arguments that cannot be evaluated; either the function must be a nimbleFunction or values for the following inputs must be specified as constants in the model: Ci[l].

Thanks a lot.

Best,

Zh

[Chris Paciorek]

hi Zh,

Nimble doesn't support the use of a stochastic index (C[i] in your case) as part of a vector of indices. Unfortunately, the error message you got is not very clear in saying that. Thanks for posting as this is something we should consider as a future feature and that we should have a better error message about.

Fortunately, you can do this in nimble by using a nimbleFunction that takes the entire vector(s) and the index and does the computation within the nimbleFunction.  Here's a simple example that you could modify to use in your situation:

dynamic_inprod <- nimbleFunction(

    run = function(vec1 = double(1), vec2 = double(1), L = double(0)) {

        returnType(double(0))

        len1 <- length(vec1)

        len2 <- length(vec2)

        if(L > len1 | L > len2)

            stop("L is greater than the length of an input vector.")

        return(inprod(vec1[1:L], vec2[1:L]))

    })

code <- nimbleCode({

    L ~ dcat(p[1:m])  # notice here that 'L' can be no bigger than 'm' by definition (unlike with a Poisson)

    S <- dynamic_inprod(a[1:m], b[1:m], L)

})

m <- nimbleModel(code, constants = list(m = 4), inits = list(L = 3, p=rep(.25, 4), a=rnorm(4), b = rnorm(4)))

You will need to do something similar if you want to assign values to a subset of Vtem as you do in "Vtem[1:Ci[l]] <- VFunNim(Ci[l], l, alpha)". Basically you can't have "1:C[i]" anywhere in your model code. One possibility would be to combine the inprod calculations with the code in VFunNim all as one nimbleFunction and fill the values of Vtem inside that nimbleFunction.

One other caution is that a Poisson random variable is unbounded so you'll need to make sure that C[i] is never larger than the length of any vector it is indexing into. 

-Chris

--
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...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/41fb6051-e14c-4584-9184-33a8f5522046n%40googlegroups.com.

[Lu]

Hi Chris,

Thanks for your advice. I will have a try.

Best,

Zh