Jags does not sample from dmnorm

Help
2011-05-24
2015-10-29
  • pinkisntwell

    pinkisntwell - 2011-05-24

    I'm trying to sample from a two-dimensional multivariate normal distribution
    with an omega matrix of (which is the inverse of ) but jags claims "Non
    positive definite matrix in call to logdet". Is this a bug? How am I supposed
    to work with multivariate normal distributions in jags?

     
  • Martyn Plummer

    Martyn Plummer - 2011-05-25

    Are you using JAGS 2.2.0? There was a bug in earlier versions but it is fixed
    now.

    If you are up to date then please send me a reproducible example with model
    file, data, initial values and script so I can look into it.

     
  • pinkisntwell

    pinkisntwell - 2011-05-25

    I tried in 2.2.0 and 2.0.0. I have emailed you the necessary files to
    reproduce the problem, but I also include their contents here for the benefit
    of those who might come across the problem in the future.

    ==> update_other_s.dat <==

    dims <- 2

    rf <- c(0,0)

    rsomega.unrolled <- c(1,33, -0.67, -0.67, 1.33)

    ==> update_other_s.mod <==

    model{

    for (i in 1:dims) {

    for (j in 1:dims) {

    rsomega <- rsomega.unrolled

    }

    }

    rs ~ dmnorm(rf, rsomega)

    }

    ==> update_other_s.scr <==

    model in update_other_s.mod

    data in update_other_s.dat

    compile

    initialize

     
  • Martyn Plummer

    Martyn Plummer - 2011-05-25

    Well there's your problem right here:

    rsomega.unrolled <- c(1,33, -0.67, -0.67, 1.33)
    

    You have a comma in place of a decimal point, so rsomega.unrolled is a vector
    of length 5 and rsomega is asymmetric and definitely not positive definite.

    In the development version of JAGS I suppressed the symmetry test (It is
    O(N^2) and slows down the updating considerably with large multivariate
    normal) and all the functions that use symmetric matrices just use the lower
    triangle, assuming symmetry. This means your example does not generate an
    error, which is a bit of a problem.

     
  • pinkisntwell

    pinkisntwell - 2011-05-25

    Well, this is embarassing to say the least... Sorry for wasting your time.

     
  • Elizaveta Semenova

    Hello, I have the same message "Non positive definite matrix in call to logdet" when a specify an informative prior in the JAGS-model

    phi ~ dgamma(1.1,1)

    Is works well for other parameters if I fix the value of this one via

    phi <-1

    Any idea on how to fix it?

     
    Last edit: Elizaveta Semenova 2015-10-06
  • Martyn Plummer

    Martyn Plummer - 2015-10-06

    I think we need more details: at least the full model description and preferably a complete worked example.

     
  • Elizaveta Semenova

    My complete model looks as written below. It works fine if I don't want to infer on 'phi' but assign it a value (phi <- 2) instead. The error appears when I include it into the list of parameters.

    model{
    
      lambda0 ~ dgamma(100,10) 
      phi ~ dgamma(40,20)
      tau.sq.Z ~ dgamma(40,20)
      sigma.sq.Z<-1./tau.sq.Z
    
      for ( i in 1 : n.z){
       for ( j in 1 : n.z){ 
           Sigma.Z[i,j]  <- sigma.sq.Z*exp(-phi*D.Z[i,j])
          }
      }
    
      P.Z  <- inverse(Sigma.Z)
    
      for (i in 1:n.z){
         #Z.mu[i] <- - sigma.sq.Z/2
         Z.mu[i] <- 0.0
      }
    
      Z ~ dmnorm(Z.mu , P.Z)
    
       for (i in 1:M){
               I.components[i] <-exp(Z[id[n+i]])
             }
    
       I <-sum(I.components[]) * D/M
    
       c <- 1000000
       for (i in 1:n){
    
       l[i] <- -lambda0 * I/n + log(lambda0) + Z[id[i]]
    
       zeros.mean[i] <- -l[i]+c 
       zeros[i] ~ dpois(zeros.mean[i])
     }
    }
    
     
    Last edit: Elizaveta Semenova 2015-10-27
  • Martyn Plummer

    Martyn Plummer - 2015-10-27

    So the context here is that you start with a distance matrix D.z and transform it to a variance-covariance matrix via

    Sigma.Z  <- sigma.sq.Z * exp(-phi * D.Z)
    

    (Giving the transformation in vectorized form for simplicity) Which must then be inverted to a precision matrix.

    I think the problem is that the matrix can become ill conditioned if phi is very small. Try bounding phi away from zero. Also verify that D.Z is indeed a distance matrix.

     
  • Elizaveta Semenova

    Thank you, Martyn, bounding phi away from zero has helped:

    phi_sample ~ dgamma(40,20)
    phi <- max(0.001, phi_sample)

     
    • Martyn Plummer

      Martyn Plummer - 2015-10-29

      Probably better to use truncation here:

      phi ~ dgamma(40, 20) T(0.001, )
      
       

Log in to post a comment.

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

Sign up for the SourceForge newsletter:





No, thanks