Range constraints for MVN nodes - a solution?

Help
2014-03-10
2014-03-19
  • Richard James

    Richard James - 2014-03-10

    Hello fellow JAGS users,

    The model I am running (see below) has a multivariate normal node "S" which I require to be constrainted to a 0-100 interval (i.e. no negative values). I've tried using the "I(0,100) and "T(0,100)" terms to no avail, and I note that the manual it states that "I" cannot be used to constrain MVN nodes. Has anyone come across a way around this issue, as currently I'm generating negative S values that are illogical?

    Thanks for any assistance given.

       model {
      # target data likelihood
      for(i in 1:N) {
        Y[i,1:J] ~ dmnorm(mu[i,],OmegaresZ)
        mu[i,1:J] <- p[i,1:K]%*%s[1:K,1:J,i]
      }
    
      # s
      for(i in 1:N) {
        for(k in 1:K) {
          s[k,1:J,i] ~ dmnorm(muS[,k],OmegaS[,,k])
        }
      }
    
      # p
      for(i in 1:N) {
        p[i,1:K] <- expphiV[i,1:K]/sum(expphiV[i,1:K])
        phiV[i,1:K] <- phi[i,1:(K-1)]%*%tV
        for(k in 1:K) {
          expphiV[i,k] <- exp(phiV[i,k])
        }
      }
    
      # phi
      for(i in 1:N) {
        for(k in 1:(K-1)) {
          phi[i,k] ~ dnorm(muphi[k],tauphi[k])
        }
        philast[i] <- -sum(phi[i,1:(K-1)])
      }
    
      # hyperparameters for phi prior
      for(k in 1:(K-1)) {
      muphi[k] ~ dnorm(0,1000)    
             tauphi[k] ~ dgamma(2,1)
        sigma2phi[k] <- 1/tauphi[k]
      }
    
      # Prior on OmegaresZ
      # combined precision matrix of residual + measurement error
      OmegaresZ ~ dwish(Rres,kres)
      SigmaresZ <- inverse(OmegaresZ) # for output
    
    }
    
     
  • Martyn Plummer

    Martyn Plummer - 2014-03-19

    The multivariate normal distribution can't be truncated. However, in your example it appears that the mean and variance of s are observed, because they do not appear on the left-hand side of any relation. If this is the case you can use the dinterval distribution to apply the constraint. Normally dinterval represents censoring, not truncation (see the manual) but when the parameters of a distribution are fixed they are the same.

    To do this you would need to add a data step, or supply the equivalent data:

    data {
       # s
       for(i in 1:N) {
          for(k in 1:K) {
             for (j in 1:J) {
                positive[i,j,k] <- 1
          }
       }
       zero <- 0
    }
    

    Then in your model

       # s
       for(i in 1:N) {
          for(k in 1:K) {
             for (j in 1:J) {
                positive[i,j,k] ~ dinterval(s[i,j,k], zero)
          }
       }
    

    This encodes the information that every element of the array s is positive.

    JAGS will fall back on a rather simple sampling method for s - the random walk Metropolis sampler so monitor the mixing of these parameters.

     
    Last edit: Martyn Plummer 2014-03-19

Log in to post a comment.

Get latest updates about Open Source Projects, Conferences and News.

Sign up for the SourceForge newsletter:

JavaScript is required for this form.





No, thanks