Errors look harmless, but model does not contain samples

[jonatho...@gmail.com]

Hello,

I'm a first time Stan user so I might be doing something obviously stupid, but I can't seem to find the problem.  I fit a large model in Rstan, the code ran for approximately 8 hours.  The errors appear pretty harmless to me but when I look at the returned object I receive the error 'Stan model 'e3e0e19d22c8ea8ec8d41b888c7c59ca' does not contain samples.'  

I found the following statement in the Rstan manual which might be relevant.

  However, if cores > 1 and there is an error for any of the chains, then the error(s) are printed. If all chains have errors and an error occurs before or during sampling, the returned object does not contain samples. But the compiled binary object for the model is still included, so we can reuse the returned object for another sampling.

But I really don't know what that means.  Is the object that I have saved useless, or can information actually be extracted?

Any advice would be greatly appreciated.

Thanks,

Jonathon O'Brien

The model statement and associated errors are below, and the full R script is attached.  

[1] "The following numerical problems occured the indicated number of times after warmup on chain 2"

                                                                                                     count

Exception thrown at line 38: multi_normal_log: Covariance matrix is not symmetric. Covariance matrix    11

Exception thrown at line 38: multi_normal_log: LDLT_Factor of covariance parameter is not positive d     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."

Chain 3, Iteration: 1800 / 2000 [ 90%]  (Sampling)

Chain 4, Iteration: 1800 / 2000 [ 90%]  (Sampling)

Chain 1, Iteration: 1600 / 2000 [ 80%]  (Sampling)

Chain 4, Iteration: 2000 / 2000 [100%]  (Sampling)

 Elapsed Time: 18845.2 seconds (Warm-up)

               7888.17 seconds (Sampling)

               26733.4 seconds (Total)

[1] "The following numerical problems occured the indicated number of times after warmup on chain 4"

                                                                                                     count

Exception thrown at line 38: multi_normal_log: Covariance matrix is not symmetric. Covariance matrix    14

Exception thrown at line 38: multi_normal_log: LDLT_Factor of covariance parameter is not positive d     4

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

Chain 3, Iteration: 2000 / 2000 [100%]  (Sampling)

 Elapsed Time: 17757.9 seconds (Warm-up)

               9347.9 seconds (Sampling)

               27105.8 seconds (Total)

Chain 1, Iteration: 1800 / 2000 [ 90%]  (Sampling)[1] "The following numerical problems occured the indicated number of times after warmup on chain 3"

                                                                                                     count

Exception thrown at line 38: multi_normal_log: Covariance matrix is not symmetric. Covariance matrix    10

Exception thrown at line 38: multi_normal_log: LDLT_Factor of covariance parameter is not positive d     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."

Chain 1, Iteration: 2000 / 2000 [100%]  (Sampling)

 Elapsed Time: 19549.6 seconds (Warm-up)

               9836.25 seconds (Sampling)

               29385.8 seconds (Total)

[1] "The following numerical problems occured the indicated number of times after warmup on chain 1"

                                                                                                     count

Exception thrown at line 38: multi_normal_log: LDLT_Factor of covariance parameter is not positive d     7

Exception thrown at line 38: multi_normal_log: Covariance matrix is not symmetric. Covariance matrix     6

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

data {

int<lower=0> npep ;

int<lower=0> N ;

int<lower=0> nprotch ;

real y[N] ;

matrix[npep,2] int_IS ;

int<lower=1,upper=npep> opep[N] ;

int<lower=1,upper=nprotch> oprotch[N] ;

}

parameters{

real beta_pep[npep] ;  // pep effect

vector[2] beta[nprotch] ;  // channel/prot averages

corr_matrix[2] G;

real<lower=0> tau ; // peptide effect var

real<lower=0> sigma ; // experimental error

vector<lower=0>[2] bsds ; // random int and slope SD

real mu_pep ;

vector[2] mu_beta ;

}

model {

beta_pep ~ normal(0,10) ;

mu_pep ~ normal(0, 5) ;

mu_beta ~ normal(0,5) ;

bsds ~ cauchy(0, 2.5) ;

beta ~ multi_normal(mu_beta, quad_form_diag(G, bsds)) ;

for (i in 1:N)

  y[i] ~ normal(beta_pep[opep[i]] + int_IS[opep[i]] * beta[oprotch[i]], sigma) ;

}

generated quantities {

  vector[nprotch] ppred;

  for (i in 1:nprotch)

  ppred[i] =  sum(beta[oprotch[i]]) ;

}

[Bob Carpenter]

I would suggest trying a small number of iterations in a single
chain to make sure everything's working. It's the only way you'll
be able to debug problems.

We need better methodology guidelines that try to tell users not
to run for 8 hours without testing the models and programs first.
I really hope they find a way to fix the bug that different chains
swallow errors when running in parallel.

You could also try running CmdStan, which is more robust for this
kind of thing in that it streams out draws as it goes and doesn't
swallow error messages.

- 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.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
> <StanModel.R>

[Jonathon O'Brien]

Hi Bob,

Thanks for the advice.  I'll give that a try this weekend.

Also, thanks for all the work and promoting of Stan that you do.  I decided to try Stan because of one your talks on Youtube.  The Gibbs samplers I work with are horrifically slow.  I'm hoping Stan is the solution.

Jonathon

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

> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
> <StanModel.R>

--
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/BHQCvk7_vbc/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

[Ben Goodrich]

On Friday, October 21, 2016 at 10:51:43 AM UTC-4, jonatho...@gmail.com wrote:

But I really don't know what that means.  Is the object that I have saved useless, or can information actually be extracted?

Pretty useless. I would guess the culprit is NaNs in the generated quantities. You can specify the sample_file argument to write the draws to the disk, which may help to figure out what went wrong.

Ben

[Bob Carpenter]

I'd strongly suggest starting with simpler models and simulated
data so you can build up confidence that the code's working, then
throw big real data sets at it.

And feel free to send your code here for comments on efficiency.
The big thing you want to do is move to the Cholesky factor
parameterization of that correlation matrix and use the _cholesky
form of the multi-normal. It's much faster and more arithmetically
stable.

- Bob

> > To unsubscribe from this group and stop receiving emails from it, 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.
> > <StanModel.R>
>
> --
> 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/BHQCvk7_vbc/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.
>
>

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

[Jonathon O'Brien]

Over the weekend I re-ran the code with a smaller dataset and only 200 iterations on a single core.  It appears to have worked (the posterior means made sense when compared to known true values) so I'm assuming the original problem had to do with the parallelization.

I am going to work on making it faster now.  The full data requires approximately 100,000 parameters and 200 iterations was not enough to obtain convergence.  I'll let you know what happens after the using the Cholesky factor parameterization and I may look into the Variational Bayes solution as well.   

Thanks again for the help, I greatly appreciate it,

Jonathon

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

> > To post to this group, send email to stan-...@googlegroups.com.
> > For more options, visit https://groups.google.com/d/optout.
> > <StanModel.R>
>
> --
> 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/BHQCvk7_vbc/unsubscribe.

> To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
>
>
> --
> 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+unsubscribe@googlegroups.com.

> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

--

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/BHQCvk7_vbc/unsubscribe.

To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

[Jonathon O'Brien]

I tried to improve the processing speed by using a Cholesky factorization on the correlation matrix.  There are now two unexpected problems that I can't quite figure out.  Somehow I'm not correctly generating N(0,1) variates and every single one of my beta random variables has identical samples.  I'm at a loss to explain what's going wrong.

I am simulating Normal(0,1) random variables to be used in the transformed parameters section.  Looking at a traceplot of these values it is clear that they are not being sampled correctly.  The full code is below but the relevant call is 

parameters{

real beta_pep[npep] ;  // pep effect

vector[2] alpha[nprotch] ;

.

.

.

model{

for (i in 1:nprotch){

  alpha[i] ~ normal(0,1) ;  // normal(0,1)s used to efficiently get beta

}

As I understand it I have created an array of 2 dimensional vectors.  Then my model statement calls each vector and assigns it as a normal(0,1).  Yet somehow in the sampling the first component is not a N(0,1).  I have attached an image showing a pair of the traceplots.  The second component of the vector looks fine.  Have I made a mistake in my model syntax?

One more quick question, I have created a cholesky_factor_corr[2] object.  In this dataset the correlation matrix has only 1 parameter.  Since there is only one parameter would it be better to just assign one parameter with a uniform prior and place it in a matrix, or do you still expect a speedup from the factorization?

Thanks,

Jonathon

data {

int<lower=0> npep ;

int<lower=0> N ;

int<lower=0> nprotch ;

real y[N] ;

matrix[npep,2] int_SSN ;

int<lower=1,upper=npep> opep[N] ;

int<lower=1,upper=nprotch> oprotch[N] ;

}

parameters{

real beta_pep[npep] ;  // pep effect

vector[2] alpha[nprotch] ;  // channel/prot averages independent

real<lower=0> tau ; // peptide effect var

real<lower=0> sigma ; // experimental error

vector<lower=0>[2] bsds ; // random int and slope SD

vector[2] mu_beta ;

cholesky_factor_corr[2] L ;

}

transformed parameters {

vector[2] beta[nprotch] ;

for (i in 1:nprotch){

  beta[i] = mu_beta + bsds .* (L * alpha[i]) ;  //int_SSN*beta gives the protch mean

}

}

model {

// First do parameters that effect the protein level

mu_beta ~ normal(0,10) ;

bsds ~ normal(0, 10) ; //This is a half-normal

L ~ lkj_corr_cholesky(1) ; 

for (i in 1:nprotch){

  alpha[i] ~ normal(0,1) ;  // normal(0,1)s used to efficiently get beta

}

//Now compute peptide deviations

tau ~ normal(0, 5) ;  //peptide prior

beta_pep ~ normal(0, tau) ;

//finally add it together

sigma ~ normal(0, 5) ;

for (i in 1:N)

  y[i] ~ normal(beta_pep[opep[i]] + int_SSN[opep[i]] * beta[oprotch[i]], sigma) ;

}  // end model statement 

generated quantities {

  vector[nprotch] ppred;

  for (i in 1:nprotch)

  ppred[i] =  beta[i][1] ;

}

[Bob Carpenter]

Can you send us the whole model where you think you're getting
an error? If there are other statements involved in the log density
that involve those variables, the posteriors of alpha won't be
normal.

- Bob

> > > To unsubscribe from this group and stop receiving emails from it, 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.
> > > <StanModel.R>
> >
> > --
> > 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/BHQCvk7_vbc/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.
> >
> >
> > --
> > 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.

> > To post to this group, send email to stan-...@googlegroups.com.
> > For more options, visit https://groups.google.com/d/optout.
>
> --
> 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/BHQCvk7_vbc/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.
>
>
>
> --
> 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.

> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

> <Rplot.png>

[Jonathon O'Brien]

Sure,

The full code is below.  

Thanks!

Jonathon

############Code for implementing the compositional proteomics model in RSTAN########

setwd('C:/Users/Thom6/Documents/TMTModeling/GoldStandard')

library(rstan)

library(rstudioapi)

rstan_options(auto_write = TRUE)

options(mc.cores = parallel::detectCores())

dat<-read.csv('sasdat_1_ms2.csv') 

#create subset approximately the size of one tenplex

uprot <- unique(dat$protid)

nprot <- length(uprot)

set.seed(777)

protlab <- sample(1:6, nprot, replace =T)

datlabels <- protlab[match(dat$protid, uprot)]

dat <- data.frame(dat, datlabels)

dat <- dat[order(dat$protid,dat$peptide,dat$treatment),]

subdat <- by(dat, dat$datlabels, function(x) x)

################Now process the subdat by group##############

dat <- subdat[[1]]   #start with group 1

nchannel <- length(unique(dat$treatment))

channel <- as.integer(dat$treatment)

N <- length(dat$protid)

y <- dat$lratio

peps <- paste(dat$protid, dat$peptide)

npep <- length(unique(peps))

opep <- rep(NA,N)

counter <- 1

id <- peps[1]

for (i in 1:N){

  if(id == peps[i]){

    opep[i] <- counter

  }else{

    counter <- counter + 1

    opep[i] <- counter

    id <- peps[i]

  }

}

protch <- paste(dat$protid, dat$treatment)

oprotch <- as.integer(factor(protch))

nprotch <- length(unique(protch))

pepIS <- aggregate(dat$IS ~ opep, FUN = mean)

int_IS <-cbind(1,pepIS[,"dat$IS"])

SSN <- 4*(dat$totsig/max(dat$totsig)) - 4

pepSSN <- aggregate(SSN ~ opep, FUN = mean)

int_SSN <-cbind(1,pepSSN[,"SSN"])

##########################Begin Stan Code##########################

modelText <- '

'  # End of STAN code

compSSN <- stan(model_code = modelText, chains = 1, iter = 1000 )  

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

> > > To post to this group, send email to stan-...@googlegroups.com.
> > > For more options, visit https://groups.google.com/d/optout.
> > > <StanModel.R>
> >
> > --
> > 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/BHQCvk7_vbc/unsubscribe.

> > To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

> > To post to this group, send email to stan-...@googlegroups.com.
> > For more options, visit https://groups.google.com/d/optout.
> >
> >
> > --
> > 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+unsubscribe@googlegroups.com.

> > To post to this group, send email to stan-...@googlegroups.com.
> > For more options, visit https://groups.google.com/d/optout.
>
> --
> 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/BHQCvk7_vbc/unsubscribe.

> To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
>
>
>
> --
> 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+unsubscribe@googlegroups.com.

> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

> <Rplot.png>

--
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/BHQCvk7_vbc/unsubscribe.

To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

[Bob Carpenter]

I'm trying to clear out my mail backlog and it looks like this one
never got answered.

I just learned that you need to run in a single chain on a single core
in order to get error messages back from RStan. I hope we can fix that
in the future, but for now, you can do that to see where things go wrong.

- Bob

> > > > To unsubscribe from this group and stop receiving emails from it, 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.
> > > > <StanModel.R>
> > >
> > > --
> > > 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/BHQCvk7_vbc/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.
> > >
> > >
> > > --
> > > 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.

> > > To post to this group, send email to stan-...@googlegroups.com.
> > > For more options, visit https://groups.google.com/d/optout.
> >
> > --
> > 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/BHQCvk7_vbc/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.
> >
> >
> >
> > --
> > 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.

> > To post to this group, send email to stan-...@googlegroups.com.
> > For more options, visit https://groups.google.com/d/optout.
> > <Rplot.png>
>
> --
> 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/BHQCvk7_vbc/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.
>
>
> --
> 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.

[Jonathon O'Brien]

Hi Bob,

Thanks for getting back to me.  I did end up running the single chain and it turned out that at some point I was exceeding my system memory.  Anyhow, I've got the program working and the results are excellent.  I think we will be building a Stan analysis into our in house analysis workflow.

Thanks,

Jonathon