Simple DLM in NIMBLE

[Colin]

Hello,

I am new to NIMBLE, and also to BUGS coding.  I am trying to program a simple DLM model with one harmonic term to get up to speed with NIMBLE.  I am having some trouble with the indexing, and am getting dimension errors--e.g., "Error in var2vertexID[1:2, 1:2] : incorrect number of dimensions."  I have played around with different set ups, but can't seem to see what's wrong.  I have two latent states, and a 2x2 transition matrix, G.  Below is my code.   Perhaps there is an obvious answer to this, but I am not seeing it at the moment.  Any tips would be much appreciated.

Thank you,

Colin

testCode <- nimbleCode({

  #Identity matrix

  W<-identityMatrix(d = 2)

  

  #F matrix

  F1<-matrix(c(1,0), nrow = 1, ncol = 2)

  

  #G Transition matrix

  G<-matrix(c(cos(2*pi/365),-sin(2*pi/365),sin(2*pi/365),cos(2*pi/365)),nrow=2,ncol=2)

  

  #initial values

  m<-G%*%x0

  var0<-sigmaSquaredInvw*W

  x[1:2,1]~dmnorm(m[1:2],var0[1:2,1:2])

  y[1]~dnorm(F1%*%x[1:2,1],sigmaSquaredInvv)

  

  #Model

  for (t in 2:T){

    mu<-G%*%x[1:2,t-1]

    sigs<-sigmaSquaredInvw*W

    x[1:2,t] ~ dmnorm(mu[1:2],sigs[1:2,1:2])

    y[t] ~ dnorm(F1%*%x[1:2,t],sigmaSquaredInvv)

  }

  

  #Initial Values for States

  c0<-c(0,0)

  m0<-sigmaSquaredInvw*W

  x0~dmnorm(c0[1:2],m0[1:2,1:2])

  

  #Priors

  sigmaSquaredInvv~dgamma(5,20)

  sigmaSquaredInvw~dgamma(5,200)

})

smosModel<-nimbleModel(code=testCode,name='2harm',constants=list(T=2190,pi=pi),data=list(y=y1),

                        inits=list(sigmaSquaredInvv=1,sigmaSquaredInvw=1))

[Perry de Valpine]

Hi Colin,

Here are a couple of things I see in your code:

On Nov 6, 2016, at 1:06 PM, Colin <clewi...@gmail.com> wrote:

Hello,

I am new to NIMBLE, and also to BUGS coding.  I am trying to program a simple DLM model with one harmonic term to get up to speed with NIMBLE.  I am having some trouble with the indexing, and am getting dimension errors--e.g., "Error in var2vertexID[1:2, 1:2] : incorrect number of dimensions."  I have played around with different set ups, but can't seem to see what's wrong.  I have two latent states, and a 2x2 transition matrix, G.  Below is my code.   Perhaps there is an obvious answer to this, but I am not seeing it at the moment.  Any tips would be much appreciated.

Thank you,

Colin

testCode <- nimbleCode({

  #Identity matrix

  W<-identityMatrix(d = 2)

  

  #F matrix

  F1<-matrix(c(1,0), nrow = 1, ncol = 2)

You can’t use c() inside the model code.  You will be able to do that in the future.  For now you will need to create a variable with the c(0,1) content, which could be created as part of your model (one element at a time) or passed in as a constant when you call nimbleModel.

  

  #G Transition matrix

  G<-matrix(c(cos(2*pi/365),-sin(2*pi/365),sin(2*pi/365),cos(2*pi/365)),nrow=2,ncol=2)

There are two issues here.  One is again that you can’t use c().  The other is that in nimble matrix() will only initialize to a single value, not a provided vector.  That is another feature we’ll be providing soon as we extend the range of basic R functions that can be compiled via NIMBLE.  So, again, you’ll need to create G outside the nimbleCode and include it as a constant when you call nimbleModel.

  

  #initial values

  m<-G%*%x0

  var0<-sigmaSquaredInvw*W

  x[1:2,1]~dmnorm(m[1:2],var0[1:2,1:2])

  y[1]~dnorm(F1%*%x[1:2,1],sigmaSquaredInvv)

At the moment it can be problematic to include a linear algebra expression as an argument to produce a scalar distribution parameter.  It would be better to use inner product, inprod(F1[1,1:2], x[1:2,1])

  

  #Model

  for (t in 2:T){

    mu<-G%*%x[1:2,t-1]

Depending on how G is provided, you might need to do G[,] or G[1:2, 1:2]

    sigs<-sigmaSquaredInvw*W

    x[1:2,t] ~ dmnorm(mu[1:2],sigs[1:2,1:2])

    y[t] ~ dnorm(F1%*%x[1:2,t],sigmaSquaredInvv)

See above

  }

  

  #Initial Values for States

  c0<-c(0,0)

  m0<-sigmaSquaredInvw*W

  x0~dmnorm(c0[1:2],m0[1:2,1:2])

  

  #Priors

  sigmaSquaredInvv~dgamma(5,20)

  sigmaSquaredInvw~dgamma(5,200)

})

smosModel<-nimbleModel(code=testCode,name='2harm',constants=list(T=2190,pi=pi),data=list(y=y1),

                        inits=list(sigmaSquaredInvv=1,sigmaSquaredInvw=1))

--
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 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/f945ecce-bee9-4767-bfe9-6b3ea22c65f4%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Colin]

Hey Perry,

Thanks so much for your comments.  Very helpful.  I really appreciate it.  

Cheers,

Colin

[Daniel Turek]

Colin, I can't fully debug or fix your code this moment, but thought a really quick reply with a few pointers might be helpful.

You're mixing up some of the concepts of "BUGS code to define the model", which defines the hierarchical model structure, and is passed as an argument to the function nimbleCode(), with "NIMBLE DSL (domain specific language) code" which is used to specify statistical algorithms in the form of nimbleFunctions, and appears as an argument to the function nimbleFunction().

It looks like you're trying to define a hierarchical model.  Here are a few things that you're doing wrong:

If you use "pi", then you'll have to pass that in as a "constant" in the constants=list(pi = 3.14...).

Can't use identityMatrix(), or matrix() in BUGS model code; only in the NIMBLE DSL (to define algorithms).  If those are fixed constant matrices, (as for W), then use W[1:2, 1:2] freely in the BUGS model code, but pass in the value of that constant matrix through the constants list constants = list(W =diag(2)).  Same goes for your other matrix, G.

If the matrix elements might vary, then define them one at a time, e.g.:

A[1,1] <- [some expression of other elements]

A[1,2] <- [some expression of other elements]

etc...

No use of c(...) function allowed, not yet.  (feature soon to come).  Instead use either:

c0[1] <- 0

c0[2] <- 0

then use c0[1:2], and c0[1], and c0[2] freely in the model code.

Alternative to those two definitions above, is to pass in the constant list(c0 = c(0,0)).

Unfortunately, your "m" term will get you into trouble here.  The definition you have:

m<-G%*%x0

defines "m" as a 2x1 matrix, not as a vector, since matrix multiplication produces a matrix.  So that's ok, but instead, use "m" in your code as:

m[1:2, 1]

That's all I see offhand.  Sorry I can't try making these changes right now, but let us know how it works.

A quick favour to ask.  We'd love to make the User Manual so clear that these things are all apparent.  If you read the User Manual, and found parts of it confusing, or that led you to some wrong conclusions, can you please email us what section(s) or specific paragraph(s) were confusing, or misleading?  That would be really helpful for us to fix them up!

Thanks,

Daniel

NIMBLE development team

--
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.

[Daniel Turek]

Looks like I was beaten to the punch here.  At least no one can say that the NIMBLE team is slow to respond.

[Colin]

Hi Daniel,

I have a quick follow up question about the best way to set up the latent nodes in a NIMBLE model step.  I am at the point where I am trying to use NimCopy to copy over simulated nodes from a ModelValues object to the model.  It looks like NIMBLE only allows you to copy by row.  So if a state is multivariate the columns will correspond to additional parameters.  

I am still working on the same model we discussed in above thread.  But the way the states are set up is the columns correspond to the latent time points, rather than the rows.  

This creates an issue when I create a modelValues object, e.g.,  modvals = modelValues(smosModel, m = 1) as it creates a 2 x 365 matrix, so you can't iterate over rows when running NimCopy.   It doesn't seem like NIMBLE allows you to take the transpose of ModelValues objects.  I have also tried to use the transpose function (t) in the model code to get the latent nodes to be 365 x 2, but that doesn't seem to work either.  I was curious if you had a suggestion that would make sense in order to use NIMBLE ModelValue objects most effectively. 

Thank you,

Colin

#Create Constant W, G and F matrices to pass to NIMBLE

W<-diag(2)

F<-matrix(c(1,0), nrow = 1, ncol = 2)

G<-matrix(c(cos(2*pi/365),-sin(2*pi/365),sin(2*pi/365),cos(2*pi/365)),nrow=2,ncol=2)

smosCode <- nimbleCode({

  

#Initial Values for States: x0 ~ Normal(c(0,0),diag(1/1000)) for precision

  m0[1] <- 0

  m0[2] <- 0

  C0[1:2,1:2] <- (1/1000)*W[1:2,1:2]

  x0[1:2] ~ dmnorm(m0[1:2],C0[1:2,1:2])

  

  #initial values

  m1[1:2,1] <- G[1:2,1:2]%*%x0[1:2]

  C1[1:2,1:2] <- sigmaSquaredInvw*W[1:2,1:2]

  x[1:2,1] ~ dmnorm(m1[1:2,1],C1[1:2,1:2]) #precision matrix

  y[1] ~ dnorm(inprod(F[1,1:2],x[1:2,1]),sigmaSquaredInvv) #precision 1/var

  

  #Model

  for (t in 2:T){

    mu[1:2,t] <- G[1:2,1:2]%*%x[1:2,t-1]

    sigs[1:2,1:2] <- sigmaSquaredInvw*W[1:2,1:2]

    x[1:2,t] ~ dmnorm(mu[1:2,t],sigs[1:2,1:2])

    y[t] ~ dnorm(inprod(F[1,1:2],x[1:2,t]),sigmaSquaredInvv)

  }

  

  sigmaSquaredInvv ~ dgamma(5,20) #in terms of shape and rate; prior on precision V

  sigmaSquaredInvw ~ dgamma(5,20) #in terms of shape and rate; prior on precision W

})

smosModel<-nimbleModel(code=smosCode,constants=list(T=365,pi=pi,G=G,W=W,F=F),data=list(y=y),

                       inits=list(sigmaSquaredInvv=1,sigmaSquaredInvw=1))

To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users...@googlegroups.com.

[Perry de Valpine]

Hi Colin,

I think the confusion here is the notion of “row” in a modelValues object.  What we mean is a complete set of states, not a row of a particular variable in nimble that might happen to be a matrix.  For modelValues purposes, you can think of all model variables as being flattened, like having their dimension attribute removed in R.

So it should work to use

copy(from = modvals, row = 1, to = smosModel, nodes = c(‘mu’,’sigs’))

Note copy is a synonym for nimCopy inside a nimbleFunction.

Also note you can omit the nodes argument and by default it will copy all variables.

Let us know if it gives you more trouble.

Perry

 

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/89da077f-7bf4-4267-a817-f77373f60b9b%40googlegroups.com.

[Colin]

Hi Perry,

Thanks for you reply, and for the clarification about how NIMBLE uses rows to denote a whole set of states from time T=1....N.  That illuminates additional problems I was encountering when trying to loop over a ModelValues object as you might with a list.  This makes me think my previous step might also not be perfectly congruous with NIMBLE ModelValue Objects.  Before I copy the states into the model, I have a function that performs FFBS to generate states from time T=0...365.  

As you can see in the code below, I am creating a model values object and then storing generated states in mvTheta as I loop over the time points.  I see now with NIMBLE this creates a ModelValues object of size 366. What I would like to do is then copy those 365 state variables (dropping the initial state 0) to a model object.  As you point out, however, when I run this code I only copy the first row of mvTheta rather than the whole set of 365.  

Do you think the best way to handle this is to loop over the Model Values object in the copy statement or somehow store all 365 states into a size 1 ModelValues object in my FFBS step? Sorry for the long post. If it would be helpful to attach my full R script, please me me know.  

Thanks again for you assistance,

Colin

modvals = modelValues(smosModel, m = 1)

nimCopy(from=test2$mvTheta,to=modvals,nodes='theta',nodesTo='x',row=1)

FFBStep<-nimbleFunction(

  setup=function(model,mvKF,nodes){

    mvThetaConfig <- modelValuesConf(vars = c('ht', 'Ht', 'theta'),types=c('double','double','double'), sizes=list(ht = 2, Ht = c(2, 2), theta = c(1,2)))

    mvTheta <- modelValues(mvThetaConfig)

    nodes <- model$expandNodeNames(nodes,sort=TRUE)

    timePoints <- length(nodes)

    len <- timePoints+1

  },

  run=function(){

    resize(mvTheta,timePoints+1)

    #Simulate Theta T

    Cho_T<-chol(mvKF['Ct',timePoints+1])

    mvTheta['theta',timePoints+1] <<- rmnorm_chol(n=1,mvKF['mt',timePoints+1],Cho_T,prec_param = FALSE)

    for (t in 1:365){

      mvTheta['ht',len-t] <<- mvKF['mt',len-t]+mvKF['Ct',len-t]%*%t(model$G)%*%inverse(mvKF['Rt',len-t+1])%*%mvTheta['theta',len-t+1]-mvKF['at',len-t+1]

      mvTheta['Ht',len-t] <<- mvKF['Ct',len-t]-mvKF['Ct',len-t]%*%t(model$G)%*%inverse(mvKF['Rt',len-t+1])%*%(model$G)%*%mvKF['Ct',len-t]

      Chold <- mvTheta['Ht',len-t]

      mvTheta['theta',len-t] <<- rmnorm_chol(n=1,mvTheta['ht',len-t],Chold,prec_param = FALSE)

    }

  }

)

[nimble.stats]

I think I see what you are saying.

Although the way you’ve done it does flip the concepts a bit, it seems reasonable and valid.  (A related design issue for NIMBLE is when and how it would be useful to make a program work across many nodes that follow the same motif, such as many latent states in a state-space model, as opposed to always treating each node as distinct.  Thoughts welcome.)

I think an alternative would be the following.  If I follow, ht is a 2x365 (or 365x2) matrix.

You could assign into pieces of it using

mvTheta[‘ht’,1][1:2, i] <- (…)

and so on for the other variables.

In a sense what I think you need here is a data structure that can hold ht, Ht and theta the way an R list would, rather than necessarily a modelValues object.  Such data structures will be in NIMBLE soon.

Also I note that the argument chold to rmnorm_chol should be a Cholesky matrix.

Perry

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/af18c2f5-7203-4bd2-a633-cddb163f9622%40googlegroups.com.

[Colin]

Hi Perry,

Thanks for your suggestion, and I'm glad to hear a NIMBLE version of List is in the works.   What I ended up doing is in the modelValuesConf, define the size of theta = c(2,366).  Then as you suggest, just update the first row, and leave the other 365 rows empty.  Note: had to resize the modelValues object to 366 rows to store the other values, ht (2x1) and Ht, which is a 2x2 matrix.

Copying over the values of theta seems to work great using the code below.  One more quick question about how NIMBLE works is in regards to the copying.  The modelValues object is 2 x 365 because there is an initial state X0, which is its own variable.  However, the theta ModelValues object has 366 columns.  When I do the copying NIMBLE seems to automatically resize the modvals object and copies all 366 states.  Is there a way to drop the first state in the copying process?  Something akin to how in R you can put mvTheta[,-1] to drop the first column?

Thanks again,

Colin

           modvals = modelValues(smosModel, m = 1)

nimCopy(from=test2$mvTheta,to=modvals,nodes='theta',nodesTo='x',row=1)

FFBStep<-nimbleFunction(

  setup=function(model,mvKF,nodes){

    mvThetaConfig <- modelValuesConf(vars = c('ht', 'Ht', 'theta'),types=c('double','double','double'), sizes=list(ht = 2, Ht = c(2, 2), theta=c(2,366)))

    mvTheta <- modelValues(mvThetaConfig)

    nodes <- model$expandNodeNames(nodes,sort=TRUE)

    timePoints <- length(nodes)

    len <- timePoints+1

  },

  run=function(){

    resize(mvTheta,timePoints+1)

    #Simulate Theta T

    Cho_T<-chol(mvKF['Ct',timePoints+1])

    mvTheta['theta',1][1:2,timePoints+1] <<- rmnorm_chol(n=1,mvKF['mt',timePoints+1],Cho_T,prec_param = FALSE)

    for (t in 1:365){

      mvTheta['ht',len-t] <<- mvKF['mt',len-t]+mvKF['Ct',len-t]%*%t(model$G)%*%inverse(mvKF['Rt',len-t+1])%*%mvTheta['theta',1][1:2,len-t+1]-mvKF['at',len-t+1]

      mvTheta['Ht',len-t] <<- mvKF['Ct',len-t]-mvKF['Ct',len-t]%*%t(model$G)%*%inverse(mvKF['Rt',len-t+1])%*%(model$G)%*%mvKF['Ct',len-t]

      Chold <- chol(mvTheta['Ht',len-t])

      mvTheta['theta',1][1:2,len-t] <<- rmnorm_chol(n=1,mvTheta['ht',len-t],Chold,prec_param = FALSE)

    }

  }

)

[Perry de Valpine]

You should be able to say

nimCopy(…, nodes = ‘theta[1:2, 2:366]’, nodesTo = ‘x[1:2, 1:365]’, row = 1)

Is that what you need?

Also nimCopy is most useful when there are nodes from multiple variables to copy.  If you are customizing to particular nodes an alternative is simply:

model[[‘theta’]][1:2, 2:366] <<- mvTheta[1, ‘x’][1:2, 1:365]

No, negative indexing doesn’t work to drop a column or row.

Perry

To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/f7014a23-050e-4c34-b136-a0612c43bed5%40googlegroups.com.

[Colin]

Perry,

Yes, this is exactly what I was looking for.  Thank you!

~Colin