Dynamic indexing in stratified mixture model

88 views
Skip to first unread message

Ethan

unread,
Aug 30, 2023, 12:38:42 PM8/30/23
to nimble-users
Hello,

I am attempting to code a stratified mixture model, where the number of components of the mixture may differ by stratum. Specifically, there may be S strata, for each s = 1, ..., S, there are K[s] mixture components. I cannot seem to get Nimble to allow me to do this no matter if I place the vector K as data or constants.

When I declare K as data, I get the error message:

Defining model
Error in getSymbolicParentNodesRecurse(x, constNames, indexNames, nimbleFunctionNames,  :
  Dynamic indexing found in a vector of indices, 1:K[l]. Only scalar indices, such as 'idx' in 'x[idx]', can be dynamic. One can instead use dynamic indexing in a vector of indices inside a nimbleFunction.


If I declare K as a constant, I get

Defining model
Error in getSymbolicParentNodesRecurse(x, constNames, indexNames, nimbleFunctionNames,  :
  getSymbolicParentNodesRecurse: dynamic indexing of constants is not allowed in K[l]. Try adding the dynamically-indexed constant as data instead (using the data argument of nimbleModel).


Below is a reproducible example.



code <- nimbleCode({
  for ( l in 1:S ) {
    gamma[s, 1:K[l]] ~ ddirch(conc[s, 1:K[l]])
    for ( k in 1:K[l] ) {
      mu[l,k] ~ dnorm(0, sd = 10)
    }
  }
  for ( i in 1:n ) {
    c[i] ~ dcat(gamma[s, 1:K[l]])
    y[i] ~ dnorm(mu[s[i], c[i]], 1)
  }
})

## putting K in data
nimconst <- list(
  n = 100, S = 2
)
nimdata <- list(
  y = rnorm(nimconst$n)
  , s = c(rep(1, floor(nimconst$n/2) ), rep(2, ceiling(nimconst$n/2) ) )
  , K = c(2, 3)
)
mu.init <- matrix(0, nrow = nimconst$S, ncol = max(nimdata$K) )
gamma.init <- matrix(0, nrow = nimconst$S, ncol = max(nimdata$K) )
for ( s in 1:nimconst$S ) {
  gamma.init[s, 1:nimdata$K[s]] <- 1/nimdata$K[s]
}
niminit <- list(gamma = gamma.init, mu = mu.init, c = rcat(nimconst$n, c(0.5, 0.5)))
model <- nimbleModel(code, nimconst, nimdata, niminit)


## putting K in constants
nimconst$K = nimdata$K
nimdata <- nimdata[-which(names(nimdata) == 'K')]
model   <- nimbleModel(code, nimconst, nimdata, niminit)



Ethan

unread,
Sep 1, 2023, 9:36:16 AM9/1/23
to nimble-users
I should mention that, although in this case one could just fit separate models, in the actual problem I have there is a hyperprior on the means, so that the overall model cannot be estimated in stages.

Chris Paciorek

unread,
Sep 2, 2023, 1:24:43 PM9/2/23
to Ethan, nimble-users
Hi Ethan,

I think you should be able to do this by specifying `K` in constants. (We definitely don't support having `K` be non-constant and then using `K` to determine the length of a model node such as `gamma[s, 1:K[l]]`.)

Before I try to dig deeper into this, I think there are some problems with some of the lines of your model code.

In:

c[i] ~ dcat(gamma[s, 1:K[l]])

1) Do you mean `s[i]` rather than `s`?
2) `l` is not defined here. 

And in:

gamma[s, 1:K[l]] ~ ddirch(conc[s, 1:K[l]])

you again have `s`, which is not a scalar.

Let me know if you still get the error after modifying the code and putting `K` in constants.

-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/99f83fcd-2080-4b0e-9bb6-1db3ee643386n%40googlegroups.com.
Reply all
Reply to author
Forward
0 new messages