Can I obtain first chain of the posterior samples if the program stopped before the second chain?

342 views
Skip to first unread message

Irene Zeng

unread,
Mar 6, 2013, 5:04:48 PM3/6/13
to stan-...@googlegroups.com
Dear Mailing users,
I am starting to use Rstan only recently. I have run 4 chains sampling in my laptop, after 4 days, it finish the first chain of sampling, but it showed up a error message " cannot allocate vector of size 8kb". I still want to have the results of the first chain, is there a way to get the first chain samples ? How would I know how much memeory I need if I need to run 4 chains ?
Attached my script in the email.
thank you, Irene  
script and results for mailing list.R

Bob Carpenter

unread,
Mar 6, 2013, 5:10:38 PM3/6/13
to stan-...@googlegroups.com
You should have the samples.csv file if you ran through
the command line, but if you're running through R, I don't think
there's a way to recover the data.

We've been suggesting for users that they run fewer iterations
first to make sure the model works and to get a look at output.
There's no way to debug if things go wrong after a 4-day run.

If you were using multivariate models, the problem
you may have been having is a memory leak we patched
for the 1.2 release which is due out any second now.
One reason it might've taken a long time is if your computer
went to swapping memory out for virtual memory on file.

- 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.
> For more options, visit https://groups.google.com/groups/opt_out.
>
>

Andrew Gelman

unread,
Mar 6, 2013, 5:30:08 PM3/6/13
to stan-...@googlegroups.com
Just in general, I'd recommend making sure everything is going ok before you start running things for days!
A

--
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.
For more options, visit https://groups.google.com/groups/opt_out.
 
 
<script and results for mailing list.R>

Ben Goodrich

unread,
Mar 6, 2013, 5:34:18 PM3/6/13
to stan-...@googlegroups.com
On Wednesday, March 6, 2013 5:04:48 PM UTC-5, Irene Zeng wrote:
Dear Mailing users,
I am starting to use Rstan only recently. I have run 4 chains sampling in my laptop, after 4 days, it finish the first chain of sampling, but it showed up a error message " cannot allocate vector of size 8kb". I still want to have the results of the first chain, is there a way to get the first chain samples ? How would I know how much memeory I need if I need to run 4 chains ?

There is an option for stan() called sample_file that can be passed a character string to write the samples to a file. In addition, it is possible to call stan() in parallel using any of the parallel frameworks that R supports. The main problem with memory in your model is that you define mu as a transformed parameter that has length 78400! Even if you wrote the samples to a disk, you could easily crash R when you loaded them into RAM. I personally would define mu as a local vector in the model section so that its samples are used but not stored. In addition, it is wasteful to embark on a 4 day run without optimizing your code. The main areas where your code would go faster are using matrix algebra rather than scalar algebra, using the vectorized forms of the distributions rather than manually looping over all the observations, and using bernoulli_logit() instead of bernoulli(inv_logit(...)).

Ben

Irene Zeng

unread,
Mar 7, 2013, 5:07:51 AM3/7/13
to stan-...@googlegroups.com
Thanks for all  your suggestions for optimizing codes.I did not expect it run so long and do not want to interrupt it in the middle :d) It looks like the model section can be similar as the bug language.
I have used multi_normal() function for the sampling , but I had a error message as below, so I change to univariate normal for the trial. Do you know why I will get such message although I just copy the code for multinormal function ? thanks! Irene
g~multi_normal(mn,invprec);
 ^-- here

DIAGNOSTIC(S) FROM PARSER:
no matches for function name="multi_normal_log"
    arg 0 type=real[1]
    arg 1 type=real[1]
    arg 2 type=matrix
unknown distribution=multi_normal
Parser expecting: "}"
> print(fit)
Stan model 'protein_code' does not contain samples.
> plot(fit)
Stan model 'protein_code' does not contain samples.

Bob Carpenter

unread,
Mar 7, 2013, 11:57:18 AM3/7/13
to stan-...@googlegroups.com
multi_normal requires a vector for the variate and location
arguments (vector). You're giving it arrays (real[1]).

So you must have something declared as

real g[N];

which should be declared as

vector[N] g;

- Bob

Irene

unread,
Mar 8, 2013, 12:05:52 AM3/8/13
to stan-...@googlegroups.com
Yes. That is what I have been used for the multiviariate model. But I still have the same error message, Irene

Bob Carpenter

unread,
Mar 8, 2013, 12:39:02 AM3/8/13
to stan-...@googlegroups.com
Could you send us the whole model?

The following message should only pop up if
g and mn are declared as arrays rather than vectors:

> > g~multi_normal(mn,invprec);
> > ^-- here
> >
> > DIAGNOSTIC(S) FROM PARSER:
> > no matches for function name="multi_normal_log"
> > arg 0 type=real[1]
> > arg 1 type=real[1]
> > arg 2 type=matrix

- Bob
> stan-users+...@googlegroups.com <javascript:>.
> > For more options, visit https://groups.google.com/groups/opt_out <https://groups.google.com/groups/opt_out>.

Irene

unread,
Mar 8, 2013, 4:23:17 AM3/8/13
to stan-...@googlegroups.com
thanks, you are right. I change the declaration for g and mn as vector[] , the message gone. Now I also want to sample 
the row vector of the matrix U[200,3] from multi_normal distribution, I have tried several version. this is the last one,  
but not work..Irene
 
protein_code <-'
data {
 int ns;          //number of records
      int nobs;           //number of observation
      int nmiss;                           //number of missing
      int ncensor;                         //number of censored
 int<lower=0> subject[ns];    //subject id
 int<lower=0> proteinid[ns];    //protein id
 int<lower=0,upper=1> classno[ns];   //normal comparing with control
 int<lower=0,upper=1> run1[ns];   //run 1
 int<lower=0,upper=1> run2[ns];
 int<lower=0,upper=1> run3[ns];
 int<lower=0,upper=1> run4[ns];
 int<lower=0,upper=1> run5[ns];
 int<lower=0,upper=1> run6[ns];
 int<lower=0,upper=1> run7[ns];
 int<lower=0,upper=1> run8[ns];
 int<lower=0,upper=1> tlab113[ns];   //label 113
 int<lower=0,upper=1> tlab114[ns];
 int<lower=0,upper=1> tlab115[ns];
 int<lower=0,upper=1> tlab116[ns];
 int<lower=0,upper=1> tlab117[ns];
 int<lower=0,upper=1> tlab118[ns];
 int<lower=0,upper=1> tlab119[ns];
 int<lower=0,upper=1> tlab121[ns];
 real logofAUC_obs[nobs];
      real logofAUC_censor[ncensor];
 real m_z_centered[ns];     //logofAUC
      corr_matrix[3] R;
      vector[3] Sigma_sd;
 cov_matrix[3] invprec;
      vector[3] mn;
 int<lower=0,upper=1> miss[ns];
      }
transformed data
   {
 matrix[3,3] T;
 T<-diag_matrix(Sigma_sd) * R * diag_matrix(Sigma_sd);
   }
parameters {
     matrix[200,3] U;     
     vector[3] g;    
     real<lower=0> ita;
     real beta2[64];
     real beta3_1;
     real beta3_2;
     real beta3_3;
     real beta3_4;
     real beta3_5;
     real beta3_6;
     real beta3_7;
     real beta3_8;
     real beta4_113;
     real beta4_114;
     real beta4_115;
     real beta4_116;
     real beta4_117;
     real beta4_118;
     real beta4_119;
     real beta4_121;
     real logofAUC_miss[nmiss];
     real alpha;
     real alpha1;
     real alpha2;
     }
 
model
  {    
       real mu[ns];
 
     for (pep in 1:ns)
      mu[pep]<-beta2[subject[pep]]+U[proteinid[pep],1]+U[proteinid[pep],2]*m_z_centered[pep]
  +U[proteinid[pep],3]*classno[pep]+beta3_1*run1[pep]+beta3_2*run2[pep]+beta3_3*run3[pep]+beta3_4*run4[pep]+beta3_5*run5[pep]+
  beta3_6*run6[pep]+beta3_7*run7[pep]+beta4_113*tlab113[pep]+beta4_114*tlab114[pep]+beta4_115*tlab115[pep]+beta4_116*tlab116[pep]+beta4_117*tlab117[pep]+beta4_118*tlab118[pep]+beta4_119*tlab119[pep]+beta4_121*tlab121[pep];
   
 beta3_1~normal(0,10);     //prior for fixed coefficient
 beta3_2~normal(0,10);
 beta3_3~normal(0,10);
 beta3_4~normal(0,10);
 beta3_5~normal(0,10);
 beta3_6~normal(0,10);
 beta3_7~normal(0,10);
 beta3_8~normal(0,10);
 beta4_113~normal(0,10);
 beta4_114~normal(0,10);
 beta4_115~normal(0,10);
 beta4_116~normal(0,10);
 beta4_117~normal(0,10);
 beta4_118~normal(0,10);
 beta4_119~normal(0,10);
 beta4_121~normal(0,10);
  ita~gamma(200,10);
     
   #hyper prior
 #for (i in 1:3)
 #g[i]~normal(mn[i],invprec[i,i]);
 
 g~multi_normal(mn,invprec);
 #for (prot in 1:200)
 #{U[prot,1]~normal(g[1],T[1,1]);
 #U[prot,2]~normal(g[2],T[2,2]);
      #U[prot,3]~normal(g[3],T[3,3]);
      #}
 for (prot in 1:200)
 U[prot,1:3]~multi_normal(g,T);
 for (sub in 1:64)
 beta2[sub]~normal(0,1);
     
 for (pep in 1:nobs)        //nobs includes number of censored
     {logofAUC_obs[pep]~normal(mu[pep],ita);  //lp__ <- lp__ + ncensor*log1m(normal_p(-9.2,mu[pep],ita))
 }
  for (pep in 1:ncensor)
 logofAUC_censor[pep]~normal(mu[pep],ita);
      for (pep in 1:nobs)
 miss[pep]~bernoulli(inv_logit(alpha + alpha1*m_z_centered[pep]+alpha2*logofAUC_obs[pep]));
 
      for (pep in 1:nmiss)
 {logofAUC_miss[pep]~normal(mu[pep],ita);
  miss[pep]~bernoulli_logit(alpha + alpha1*m_z_centered[pep]+alpha2*logofAUC_miss[pep]);
 }
   }
'

Ben Goodrich

unread,
Mar 8, 2013, 10:02:21 AM3/8/13
to stan-...@googlegroups.com
On Friday, March 8, 2013 4:23:17 AM UTC-5, Irene wrote:
thanks, you are right. I change the declaration for g and mn as vector[] , the message gone. Now I also want to sample 
the row vector of the matrix U[200,3] from multi_normal distribution, I have tried several version. this is the last one,  
but not work..Irene

I don't think multi_normal is defined for row vectors. Only column vectors and matrices.

Ben
 

Bob Carpenter

unread,
Mar 8, 2013, 11:53:02 AM3/8/13
to stan-...@googlegroups.com
Correct. You can find a convenient list of all the
signatures for all the functions supported by Stan in
the index.

So what you need to do is apply a transpose. So if x
is the column vector, instead of writing "x ~ multi_normal(...)"
write "x' ~ multi_normal(...)". Same for the location vector,
which also has to be a vector.

We'll probably generalize the multivariate normal when
we get around to vectorizing it.

- Bob

Irene

unread,
Mar 9, 2013, 11:02:52 PM3/9/13
to stan-...@googlegroups.com
Yes, thanks. I found the optimized code section in the manual, and it is very useful. Irene
Reply all
Reply to author
Forward
0 new messages