error in applying Stan to multivariate Poisson-log normal model

[Julian]

Hello,

I started to learn Stan recently and I am applying it to obtaining inference of some count data via a multivariate Poisson-log normal model with an inverse Wishart prior.

 

My Stan code seemed good after checking, but when I plug in data in R, the fit won't work.

 

Error: Error in compileCode(f, code, language = language, verbose = verbose) : Compilation ERROR, function(s)/method(s) not created!

---

Stan code:

data {

  int<lower=0> J; // types of count data

  int<lower=0> y[J]; /

  cov_matrix[J] S; // covariance matrix in inverse wishart prior

}

parameters {

  vector<lower=0>[J] lambda;

  vector[J] mu; 

  cov_matrix[J] sig_mvn;// covariance matrix of multivariate normal distribution

}

 

model {

  sig_mvn~inv_wishart(J,S);

 

  for (j in 1:J)

   mu[j]~normal(0,0.0001);

 

  log(lambda)~multi_normal(mu,sig_mvn);

 

  target += -log(fabs(lambda));

 

  for (j in 1:J)

   y[j]~poisson(lambda[j]);

}

 

----

R code:

library(rstan)

outages_dat <- list(J = 1,y = c(280,380),S=matrix(c(1,0,0,1),nrow=2,ncol=2))

fit <- stan(file = 'outages.stan', data = outages_dat,

            iter = 2000, chains = 4)

print(fit)

 

I am sorry if the problem seems too obvious to you. Any advice on improving the code would be deeply appreciated! Thank you very much!

 

Julian Yu

----------------------------

multivariate_Poisson_log_normal.stan

data {
  int<lower=0> J; // types of count data

  int<lower=0> y[J]; /

  cov_matrix[J] S; // covariance matrix in inverse wishart prio 
}
parameters {

  vector<lower=0>[J] lambda;

  vector[J] mu; 

  cov_matrix[J] sig_mvn;// covariance matrix of multivariate normal distribution
}
model {

  sig_mvn~inv_wishart(J,S);

  for (j in 1:J)
   mu[j]~normal(0,0.0001);

  log(lambda)~multi_normal(mu,sig_mvn);

  target += -log(fabs(lambda));

  for (j in 1:J)

   y[j]~poisson(lambda[j]);
}

[Ben Goodrich]

On Monday, April 10, 2017 at 12:01:44 PM UTC-4, Julian wrote:

 My Stan code seemed good after checking, but when I plug in data in R, the fit won't work.

 

Error: Error in compileCode(f, code, language = language, verbose = verbose) : Compilation ERROR, function(s)/method(s) not created!

We need the entire error message, but it is likely the case that you do not have a C++ toolchain setup (correctly).

Ben

[Jinzhu]

Hi Ben,

Thank you for your reply. I think you are right about the reason for the error. I reinstall Rtools and ran the code again. This time the error meassage showed first was:

DIAGNOSTIC(S) FROM PARSER:
Warning (non-fatal):
Left-hand side of sampling statement (~) may contain a non-linear transform of a parameter or local variable.
If it does, you need to include a target += statement with the log absolute determinant of the Jacobian of the transform.
Left-hand-side of sampling statement:
    log(lambda) ~ multi_normal(...)

which I had already included in my Stan code, so it continued. Then there was another error meassage after I ran the fit code:

Chain 4, Iteration: 2000 / 2000 [100%]  (Sampling)
 Elapsed Time: 2.334 seconds (Warm-up)
               0.966 seconds (Sampling)
               3.3 seconds (Total)

[1] "The following numerical problems occured the indicated number of times on chain 4"
                                                                                                     count
Exception thrown at line 15: inv_wishart_log: LDLT_Factor of random variable is not positive definit     3
[1] "When a numerical problem occurs, the Hamiltonian proposal gets rejected."
[1] "See http://mc-stan.org/misc/warnings.html#exception-hamiltonian-proposal-rejected"
[1] "If the number in the 'count' column is small,  do not ask about this message on stan-users."
Warning messages:
1: In readLines(file, warn = TRUE) :
  incomplete final line found on 'C:\Users\Jinju\Desktop\UQ_termPaper\outages.stan'
2: There were 12 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup 
3: Examine the pairs() plot to diagnose sampling problems

Thanks!

Julian

--
Stan code:
data {
  int<lower=0> J; // types of count data

  int<lower=0> y[J]; /

  cov_matrix[J] S; // covariance matrix in inverse wishart prio
}
parameters {

  vector<lower=0>[J] lambda;

  vector[J] mu;

  cov_matrix[J] sig_mvn;// covariance matrix of multivariate normal distribution
}
model {

  sig_mvn~inv_wishart(J,S);

  for (j in 1:J)
   mu[j]~normal(0,0.0001);

  log(lambda)~multi_normal(mu,sig_mvn);

  target += -log(fabs(lambda));

  for (j in 1:J)

   y[j]~poisson(lambda[j]);
}

R code:

outages_dat <- list(J = 2,y = c(280,380),S=matrix(c(1,0,0,1),nrow=2,ncol=2))

fit <- stan(file = 'outages.stan', data = outages_dat,
            iter = 2000, chains = 4)
print(fit)

plot(fit)

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/hJh0fQzK6R4/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.
To post to this group, send email to stan-...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Ben Goodrich]

On Monday, April 10, 2017 at 1:38:49 PM UTC-4, Julian wrote:

Thank you for your reply. I think you are right about the reason for the error. I reinstall Rtools and ran the code again. This time the error meassage showed first was:

None of these are error messages.
 

DIAGNOSTIC(S) FROM PARSER:
Warning (non-fatal):
Left-hand side of sampling statement (~) may contain a non-linear transform of a parameter or local variable.
If it does, you need to include a target += statement with the log absolute determinant of the Jacobian of the transform.
Left-hand-side of sampling statement:
    log(lambda) ~ multi_normal(...)

which I had already included in my Stan code, so it continued.

Yes, although this can (and should) be avoided by declaring

vector[J] log_lambda;

in the parameters block and

vector[J] lambda = exp(log_lambda);

in the transformed parameters block and

target += multi_normal(log_lambda | mu, sig_mvn);

in the model block.

Although it is better still to declare a vector of standard deviations and a Cholesky factor of a correlation matrix in the parameters block rather than a covariance matrix, git rid of the inverse Wishart prior in favor of some prior on the standard deviations and a LKJ prior on the correlation matrix. Also, it is often better to compose log_lambda in the transformed parameters block as

vector[J] log_lambda = mu + diag_pre_multiply(L, sds) * z;
vector<lower=0>[J] lambda = exp(log_lambda);

where z is declared in the parameters block as

vector[J] z;

and has a standard normal prior in the model block

target += normal_lpdf(z | 0, 1);
 

[1] "The following numerical problems occured the indicated number of times on chain 4"
                                                                                                     count
Exception thrown at line 15: inv_wishart_log: LDLT_Factor of random variable is not positive definit     3

That is what happens when you put inverse Wishart priors on covariance matrices, and you were lucky it only happened on 3 iterations.
 

[1] "When a numerical problem occurs, the Hamiltonian proposal gets rejected."
[1] "See http://mc-stan.org/misc/warnings.html#exception-hamiltonian-proposal-rejected"

Read that link
 

[1] "If the number in the 'count' column is small,  do not ask about this message on stan-users."
Warning messages:
1: In readLines(file, warn = TRUE) :
  incomplete final line found on 'C:\Users\Jinju\Desktop\UQ_termPaper\outages.stan'

Make your outages.stan file end in a blank line
 

2: There were 12 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup 

Read that link
 

3: Examine the pairs() plot to diagnose sampling problems

Look at that pairs plot, especially the diagonal elements of sig_mvn

Ben

[Jinzhu]

Hi Ben,

That's really helpful! Thank you so much.

Julian

[Jinzhu]

Good morning Ben,

In the last email, you suggested me to compose log_lambda in the transformed parameters block as

     vector[J] log_lambda = mu + diag_pre_multiply(L, sds) * z;

Could you tell me more about how to specify L and sds?

Also, since I know lambda usually belongs to (1,100), is it appropriate to define vector<lower=0>mu in my stan code? Or how can we formulate a better hyper prior for mu given the information？

Thank you very much!

Julian

---

在 2017/4/10 12:58, Ben Goodrich 写道:

[Bob Carpenter]

Since you talk about log_lambda, then presumably lambda has to be positive.
But given that log_lambda is a derived quantity from mu, L, sds, and z,
you want to put priros on those, not on log_lambda. If you try to put
priors on log_lambda, you have a very messy Jacovian to account for.

You can check out the manual for a description of this non-centered
parameterization for multivariate data.

- Bob

> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.

[Jinzhu]

Hi Bob,

Thanks for the pointer. I studied the multivariate parameterization section of the reference manual (Version 2.14.0) and modified my code following the second example in that section, but there was also an error. Some in the group have also encounted such error, but  I found the reasons for theirs bugs were different, so I write to you for help. Thanks.

Julian

SYNTAX ERROR, MESSAGE(S) FROM PARSER:

variable "cholesky_decompose" does not exist.

ERROR at line 10

  9:    transformed parameters{
 10:      matrix[J,J]L=cholesky_decompose[S];
                                         ^
 11:      //L*L'=sigma of mnv reference manual 2.14 p327

Error in stanc(model_code = paste(program, collapse = "\n"), model_name = model_cppname,  : 
  failed to parse Stan model 'outages -removeWishart' due to the above error.
-----
Stan code:

data {

  int<lower=0> J; // number of types of unavailabilities
  int<lower=0> N; // number of data points of each type in a site
  int<lower=0> y[N,J]; //
  //vector[N]<lower=0> y[J]; // J1:unscheduled, J2: unavailabilities for N docks. saved for later use.
  cov_matrix[J] S; //
}

transformed parameters{

  matrix[J,J]L=cholesky_decompose[S];
  //L*L'=sigma of mnv. reference manual 2.14, page 327
}


parameters {

  vector[J] log_lambda;
  vector[J] mu;  
  // matrix[J,J] sig_mvn;
  // precision matrix Omega(the inverse of covariance) in multi_normal_preci//sion   //distri.
  vector[J] z;
} 

transformed parameters{ 

  vector[J] log_lambda;
  log_lambda= mu + L*z;
  //implies: log_lambda ~mvn(mu, sigma) where sigma=L*z
  // to put prior on mu L,z
  //otherwise, a very messy Jacovian to account for
  vector<lower=0>[J] lambda = exp(log_lambda);
}

model {
  target += normal_lpdf(z|0,1);
  
  for (j in 1:J){
  target+=normal_lpdf(mu[j]|0,1e5); 
  }
  
  for (i in 1:N){
    for (j in 1:J){
     target += poisson_lpmf(y[i,j]|lambda[j]);    
    }
  }
}

----------------------

在 2017/4/17 13:55, Bob Carpenter 写道:

Since you talk about log_lambda, then presumably lambda has to be positive.
But given that log_lambda is a derived quantity from mu, L, sds, and z,
you want to put priros on those, not on log_lambda.  If you try to put
priors on log_lambda, you have a very messy Jacovian to account for.

You can check out the manual for a description of this multivariate
parameterization for multivariate data.

- Bob

[Bob Carpenter]

You need to use parens for function arguments, as in

cholesky_decompose(S)

I'd also advise adding some space between symbols---it's free!

- Bob

> On Apr 19, 2017, at 12:27 AM, Jinzhu <yujin...@gmail.com> wrote:
>
> ...