Can I optimize my complex, and slow to fit model?

[Josef Fruehwald]

Hello everyone,

I have to apologize ahead of time that I am still very much a beginner at graphical models, bayesian statistics, and obviously Stan. My problem is that I have a Stan model that's taking forever to fit. It's basically taking a full working day for it to run through 3 chains of 1000 iterations. It might just the fact that for this data for this model, that's how long it takes, but I'd appreciate any advice on optimizing the model.

I've attached my data as an R workspace. The response variable is F1.n, which is an acoustic measurement of a specific vowel. I have observations from 308 speakers, who spoke this vowel in 429 different words. I am interested in the rate of change of F1.n in two different linguistic contexts over dates of birth. Here's what a basic plot of speaker means looks like.

I'm specifically interested in whether the rates of change are ever the same between NContext==1 and NContext==2, i.e. are they ever parallel? In this case it seems obvious that they aren't, but I'll be looking at other cases where this isn't so obvious. Also, in this case, men and women look to be basically the same, there are other cases I'll look at where they're definitely not.

I've decided to model the rate of change using b-splines, because I understand how to implement them. So, I have a parameter in the model delta[i], which represents the amount of change from DOB==-i-1 to DOB==1, which is modeled with a b-spline. I also have a parameter Delta[i], which is equal to sum(delta[1:i]), which represents the total displacement from DOB==1. Here's my Stan model.

model_code <- "

data{

  int<lower = 0> N;           // number of observations

  real y[N];                  // y
  int<lower = 0> speaker[N];  // speaker index

  int<lower = 0> word[N];     // word index

  int<lower = 0> context[N];  // context index

  int<lower = 0> max_context; // max context

  int<lower = 0> max_speaker;      // number of speakers

  int<lower = 0> dob[max_speaker]; // by speaker dob

  int<lower = 0> sex[max_speaker]; // by speaker sex

  int<lower = 0> max_sex;          // max sex

  int<lower = 0> n_bspline;        // number of b-spline knots

  int<lower = 0> max_dob;          // max_dob

 

  matrix[max_dob, n_bspline] b_spline;   // bsplines

  int<lower = 0> max_word;        // max word

}

parameters{

  vector[n_bspline] b_coefs;           // bspline coefs for sex==1, context==1

  vector[n_bspline] male_b_coefs;      // bspline coefs for sex==2, context==1

  vector[n_bspline] c2_b_coefs;        // bspline coefs for sex==1, context==2

  vector[n_bspline] c2_male_b_coefs;   // bspline coefs for sex==2, context==2

  matrix[max_speaker, max_context] speaker_mu;   // by-speaker mu

  real<lower = 0> speaker_sigma[max_speaker];    // by-speaker sigma

  real<lower = 0> sigma;              // sigma for distribution of speaker_mu

  real intercept;                    // intercept for speaker_mu

  real male_effect;                   // male effect for speaker_mu

  real context_effect;                // context effect for speaker_mu

  real word_mu[max_word];            // by-word mu

  real<lower = 0> words_sigma;       // sigma for distribution of word_mu

}

transformed parameters{

  // delta[i] = diff in y from dob==i-1 to dob==i

  vector[max_dob] delta;           // delta for for sex==1, context==1     

  vector[max_dob] male_delta;      // delta for for sex==2, context==1

  vector[max_dob] c2_delta;        // delta for for sex==1, context==2

  vector[max_dob] c2_male_delta;   // delta for for sex==2, context==2

  // Delta = total displacement from dob==1

  real Delta[max_dob, max_sex, max_context]; 

  // mu for intercept prior

  real likely_start;

  // calculate deltas

  delta <- b_spline * b_coefs;

  male_delta <- b_spline * male_b_coefs;

  c2_delta <- b_spline * c2_b_coefs;

  c2_male_delta <- b_spline * c2_male_b_coefs;

  // This whole block is basicaly cumsum(delta[1:i])

  Delta[1,1,1] <- delta[1];

  Delta[1,1,2] <- delta[1] + c2_delta[1];

  Delta[1,2,1] <- delta[1] + male_delta[1];

  Delta[1,2,2] <- delta[1] + c2_delta[1] + male_delta[1] + c2_male_delta[1];

  for(i in 2:max_dob){

    Delta[i,1,1] <- Delta[i-1,1,1] + delta[i];

    Delta[i,1,2] <- Delta[i-1,1,2] + delta[i] + c2_delta[i];

    Delta[i,2,1] <- Delta[i-1,2,1] + delta[i] + male_delta[i];

    Delta[i,2,2] <- Delta[i-1,2,2] + delta[i] + c2_delta[i] + male_delta[i] + c2_male_delta[i];
  }

  //end cumsum

  // Set mu for intercept prior to sample mean

likely_start <- mean(y);

}

model{

  // sigma prior as suggested by Stan reference

  words_sigma ~ cauchy(0, 2.5) T[0,];  // sigma for distribution of word effects

  for(i in 1:max_word){

     word_mu[i] ~ normal(0, words_sigma);  // word effect

   }

  // Calculate a mu for every speaker for each context

  for(i in 1:max_speaker){

    for(j in 1:max_context){

      speaker_mu[i,j] ~ normal(Delta[dob[i],sex[i],j] + (male_effect * (sex[i]-1)) + (context_effect * (j-1)) + intercept, sigma);

    }

    // Each speaker can have a different sigma for their data.

    speaker_sigma[i] ~ cauchy(0, 2.5) T[0,];

  }

  // priors for parameters in speaker_mu estimate

  intercept ~ normal(likely_start, 100);

  male_effect ~ normal(0, 100);

  context_effect ~ normal(0,100);

  // b_coefs

  for(i in 1:n_bspline){

    b_coefs[i] ~ normal(0, 100);

    male_b_coefs[i] ~ normal(0, 100);

    c2_b_coefs[i] ~ normal(0, 100);

    c2_male_b_coefs[i] ~ normal(0, 100);

  }

  // data loop

  // mu = speaker_mu + word_mu

  // sigma = speaker_sigma

  for(i in 1:N){

    y[i] ~ normal(speaker_mu[speaker[i], context[i]] + word_mu[word[i]], speaker_sigma[speaker[i]]);
  }

}

"

And my code for running it in rstan:

load("data.Rdata")

library(rstan)

library(splines)

library(plyr)

spinfo <- unique(data[,c("NSpeaker","DOB1","NSex")])

spinfo <- arrange(spinfo, NSpeaker)

b_spline <- cbind(1, bs(1:max(spinfo$DOB1), df = 3))

dat <- list( N = nrow(data),

             y = data$F1.n,

             speaker = data$NSpeaker,

             word = data$NWord,

             context = data$NContext,

             max_context = max(data$NContext),

             

             max_speaker = max(data$NSpeaker),

             dob = spinfo$DOB1,

             sex = spinfo$NSex,

             

             

             max_sex = max(spinfo$NSex),

             

             n_bspline = ncol(b_spline),

             max_dob = max(spinfo$DOB1),

             b_spline = b_spline,

             

             max_word = max(data$NWord)

)

model <- stan(fit = mod_init,

              data = dat, 

              iter = 1000,

              chains = 3)

To examine the difference in the rate of change between NContext==1 and NContext==2 for women, I'll be looking at c2_delta, and for men, c2_delta + c2_male_delta.

But, as I said, this is taking forever to fit. If it's just a fact that this model will take forever, given the number of parameters involved and the size of the data, then I'll have to live with it. At least after waiting a day, the result is convergence on all parameters, which is more than I can say for my JAGS implementation. Here's my R session info, if that makes a difference.

> sessionInfo()

R version 2.15.1 (2012-06-22)

Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)

locale:

[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:

[1] graphics  grDevices utils     datasets  stats     methods   base     

other attached packages:

[1] ggplot2_0.9.1

loaded via a namespace (and not attached):

 [1] colorspace_1.1-1   dichromat_1.2-4    digest_0.5.2       grid_2.15.1       

 [5] labeling_0.1       MASS_7.3-18        memoise_0.1        munsell_0.3       

 [9] plyr_1.7.1         proto_0.3-9.2      RColorBrewer_1.0-5 reshape2_1.2.1    

[13] scales_0.2.1       stringr_0.6       

Thanks!

-Joe

[Josef Fruehwald]

Apologies! My description of delta[i] got all messed up. It should be:

... I have a parameter in the model delta[i], which represents the amount of change from DOB==i-1 to DOB==i, which is modeled with a b-spline...

[Marcus Brubaker]

Hi Josef,

There are a number of things you can do to speed things up with your model.  One big thing you can do is use the vectorized form of the normal distribution.

E.G.,  instead of

  for(i in 1:max_word){

     word_mu[i] ~ normal(0, words_sigma);  // word effect

   }

and

  // b_coefs

  for(i in 1:n_bspline){

    b_coefs[i] ~ normal(0, 100);

    male_b_coefs[i] ~ normal(0, 100);

    c2_b_coefs[i] ~ normal(0, 100);

    c2_male_b_coefs[i] ~ normal(0, 100);

  }

you can simply write

word_mu ~ normal(0,words_sigma);
b_coefs ~ normal(0,100);
male_b_coefs ~ normal(0,100);
c2_b_coefs ~ normal(0,100);
c2_male_b_coefs ~ normal(0,100);

which should be a fair bit faster.

Your cumsum would probably be faster done as follows:

Instead of

  // This whole block is basicaly cumsum(delta[1:i])

  Delta[1,1,1] <- delta[1];

  Delta[1,1,2] <- delta[1] + c2_delta[1];

  Delta[1,2,1] <- delta[1] + male_delta[1];

  Delta[1,2,2] <- delta[1] + c2_delta[1] + male_delta[1] + c2_male_delta[1];

  for(i in 2:max_dob){

    Delta[i,1,1] <- Delta[i-1,1,1] + delta[i];

    Delta[i,1,2] <- Delta[i-1,1,2] + delta[i] + c2_delta[i];

    Delta[i,2,1] <- Delta[i-1,2,1] + delta[i] + male_delta[i];

    Delta[i,2,2] <- Delta[i-1,2,2] + delta[i] + c2_delta[i] + male_delta[i] + c2_male_delta[i];
  }

  //end cumsum

Use

  for(i in 1:max_dob){
    t0 <- sum(delta[1:i]);
    t1 <- sum(c2_delta[1:i]);
    t2 <- sum(male_delta[1:i]);

    Delta[i,1,1] <- t0

    Delta[i,1,2] <- t0 + t1;

    Delta[i,2,1] <- t0 + t2;

    Delta[i,2,2] <- Delta[i,1,2] + t2 + sum(c2_male_delta[1:i]);
  }

There may be a few others that can be done now but that's a start.  Right now indexing limitations will prevent making other big gains (e.g., vectorizing y or speaker_mu) but become available in the future.

Hope this helps.

Cheers,
Marcus

[Josef Fruehwald]

Thanks! The vectorized distributions worked, but this didn't.

  for(i in 1:max_dob){
    t0 <- sum(delta[1:i]);
    t1 <- sum(c2_delta[1:i]);
    t2 <- sum(male_delta[1:i]);

    Delta[i,1,1] <- t0;

    Delta[i,1,2] <- t0 + t1;

    Delta[i,2,1] <- t0 + t2;

    Delta[i,2,2] <- Delta[i,1,2] + t2 + sum(c2_male_delta[1:i]);
  }

I get the following error:

> stanc(model_code = model_code)

TRANSLATING MODEL 'model_code' FROM Stan CODE TO C++ CODE NOW.
Error in stanc(model_code = model_code) : c++ exception (unknown reason)

 

I've tried declaring t0, t1, t2 inside of the for loop and higher up in the block with the other declarations, and neither works. Is it possible that there's something wrong with the indexing, or with sum()? I've declared these deltas as vectors, if that matters.

-Joe

[Josef Fruehwald]

Even if I relax how I want to sum up my deltas, I still get the c++ error:

transformed parameters{

  ...

  vector[max_dob] delta;

  vector[max_dob] male_delta;

  

  vector[max_dob] c2_delta;

  vector[max_dob] c2_male_delta;

  real Delta[max_dob, max_sex, max_context];  

  ...

  delta <- b_spline * b_coefs;

  male_delta <- b_spline * male_b_coefs;

  c2_delta <- b_spline * c2_b_coefs;

  c2_male_delta <- b_spline * c2_male_b_coefs;

  ...

  for(i in 1:max_dob){

    Delta[i,1,1] <- sum(delta[1:i]);

    Delta[i,1,2] <- sum(c2_delta[1:i]);

    Delta[i,2,1] <- sum(male_delta[1:i]);

    Delta[i,2,2] <- sum(c2_male_delta[1:i]);

  }

...

}

[Bob Carpenter]

Can you send me the rest of the model so I can
track down the error?

My guess is it's the same bug I just patched.
The next release (1.0.2) will fix a bug in accessing
arrays of matrices and vectors. We hope to have
it out in a couple of days.

In 1.0.1, this won't work:

matrix[M,N] x[L];

x[l,m,n] <- y;

because 1.0.1 has a bug in type inference that
type inference for x[l].

Thanks,

- Bob

[Josef Fruehwald]

Here's the full model that's causing the error:

model_code <- "

data{

    int<lower = 0> N;           // number of observations

    real y[N];                  // y

    int<lower = 0> speaker[N];  // speaker index

    int<lower = 0> word[N];     // word index

    int<lower = 0> context[N];  // context index

    int<lower = 0> max_context;    // max context

    int<lower = 0> max_speaker;      // number of speakers

    int<lower = 0> dob[max_speaker]; // by speaker dob

    int<lower = 0> sex[max_speaker]; // by speaker sex

    

        

    int<lower = 0> max_sex;          // max sex

    

    int<lower = 0> n_bspline;        // number of b-spline knots

    int<lower = 0> max_dob;          // dob index

     

    matrix[max_dob, n_bspline] b_spline;   // bsplines

    

    int<lower = 0> max_word;

}

parameters{

    vector[n_bspline] b_coefs;

    vector[n_bspline] male_b_coefs;

    vector[n_bspline] c2_b_coefs;

    vector[n_bspline] c2_male_b_coefs;

    

    real<lower = 0> sigma;

    real intercept;    

    real male_effect;

    matrix[max_speaker, max_context] speaker_mu;

    real<lower = 0> speaker_sigma[max_speaker];

    

    real word_mu[max_word];

    real<lower = 0> words_sigma;

    

    real context_effect;

}

transformed parameters{

    vector[max_dob] delta;

    vector[max_dob] male_delta;

    

    vector[max_dob] c2_delta;

    vector[max_dob] c2_male_delta;

    real Delta[max_dob, max_sex, max_context];    

    real likely_start;

    

    delta <- b_spline * b_coefs;

      male_delta <- b_spline * male_b_coefs;

    c2_delta <- b_spline * c2_b_coefs;

    c2_male_delta <- b_spline * c2_male_b_coefs;

    for(i in 1:max_dob){        

        Delta[i,1,1] <- sum(delta[1:i]);

        Delta[i,1,2] <- sum(c2_delta[1:i]);    

        Delta[i,2,1] <- sum(male_delta[1:i]);

        Delta[i,2,2] <- sum(c2_male_delta[1:i]);

    }

    //Delta[1,1,1] <- delta[1];

    //Delta[1,1,2] <- delta[1] + c2_delta[1];

    //Delta[1,2,1] <- delta[1] + male_delta[1];

    //Delta[1,2,2] <- delta[1] + c2_delta[1] + male_delta[1] + c2_male_delta[1];

    //

    //for(i in 2:max_dob){

    //    Delta[i,1,1] <- Delta[i-1,1,1] + delta[i];

    //    Delta[i,1,2] <- Delta[i-1,1,2] + delta[i] + c2_delta[i];        

    //    Delta[i,2,1] <- Delta[i-1,2,1] + delta[i] + male_delta[i];        

    //    Delta[i,2,2] <- Delta[i-1,2,2] + delta[i] + c2_delta[i] + male_delta[i] + c2_male_delta[i];        

    //}

    likely_start <- mean(y);

}

model{

    words_sigma ~ cauchy(0, 2.5) T[0,];

    word_mu ~ normal(0, words_sigma);

        

    for(i in 1:max_speaker){

        for(j in 1:max_context){

            speaker_mu[i,j] ~ normal(Delta[dob[i],sex[i],j] + (male_effect * (sex[i]-1)) + (context_effect * (j-1)) + intercept, sigma);

        }

        speaker_sigma[i] ~ cauchy(0, 2.5) T[0,];

    }

    intercept ~ normal(likely_start, 100);

    male_effect ~ normal(0, 100);    

    context_effect ~ normal(0,100);

    b_coefs ~ normal(0, 100);    

    male_b_coefs ~ normal(0, 100);

    c2_b_coefs ~ normal(0, 100);    

    c2_male_b_coefs ~ normal(0, 100);        

    

    for(i in 1:N){   

        y[i] ~ normal(speaker_mu[speaker[i], context[i]] + word_mu[word[i]], speaker_sigma[speaker[i]]);

    }

}

"

-Joe

[Bob Carpenter]

I strongly recommend using the

file=<filepath>

argument to stan() in RStan to specify a model
from a file rather than the

model_code=<string>

argument to specify the model as a string.

Whatever you do, don't cut and paste from a file
into a string. Not only is it tedious, it's
error prone (c.f., Fitts's law).

I think people may not realize the file option is
there because the first example in the RStan getting
started guide uses a string rather than a file, as does
the example in the draft appendix of the 3rd edition
of BDA.

- Bob

[Josef Fruehwald]

Can I expect different behavior from the 

    file = <filepath>

approach vs the

    model_code = <string>

approach? I realized there was the file option, but opted for the string approach because I preferred to keep as much of my code for this analysis in one place as possible. 

-Joe

[Bob Carpenter]

On 9/21/12 2:04 PM, Josef Fruehwald wrote:> |
> data{
> int<lower = 0> N; // number of observations
> real y[N]; // y

...

Thanks for sending the model. In the version you sent me,
I'm not getting an issue with C++, but rather with parsing:

============================================================

--- Translating Stan graphical model to C++ code ---
bin/stanc /Users/carp/temp2/fruehwald.stan --o=/Users/carp/temp2/fruehwald.cpp
Model name=fruehwald_model
Input file=/Users/carp/temp2/fruehwald.stan
Output file=/Users/carp/temp2/fruehwald.cpp

ERROR PARSING
LOCATION: file=input; line=61, column=36

Delta[i,1,1] <- sum(delta[1:i]);

^-- here

DIAGNOSTIC(S) FROM PARSER:

Parser expecting: "]"

=============================================================

Stan doesn't support an R-style : operation inside
of indexes. We're thinking about adding that option
some day, but it's not among our top priorities yet.

The loop is also problematic in terms of efficiency.
It winds up being quadratic because it does summations
for vectors of size 1:max_dob, which takes on
order max_dob * max_dob operations.
So I commented it out:

// for(i in 1:max_dob){
// Delta[i,1,1] <- sum(delta[1:i]);
// Delta[i,1,2] <- sum(c2_delta[1:i]);
// Delta[i,2,1] <- sum(male_delta[1:i]);
// Delta[i,2,2] <- sum(c2_male_delta[1:i]);
// }

and went with the way you had commented out:

Delta[1,1,1] <- delta[1];

Delta[1,1,2] <- delta[1] + c2_delta[1];

Delta[1,2,1] <- delta[1] + male_delta[1];

Delta[1,2,2] <- delta[1] + c2_delta[1] + male_delta[1] + c2_male_delta[1];

for(i in 2:max_dob){

Delta[i,1,1] <- Delta[i-1,1,1] + delta[i];

Delta[i,1,2] <- Delta[i-1,1,2] + delta[i] + c2_delta[i];

Delta[i,2,1] <- Delta[i-1,2,1] + delta[i] + male_delta[i];

Delta[i,2,2] <- Delta[i-1,2,2] + delta[i] + c2_delta[i] + male_delta[i] + c2_male_delta[i];
}

That version works on my machine, though your mileage
may vary.

There may be a problem with the version of C++ you're
using. Older versions have trouble with some of the
template metaprograms, so could raise errors we're not
anticipating.

- Bob

[Bob Carpenter]

On 9/21/12 2:50 PM, Josef Fruehwald wrote:
> Can I expect different behavior from the
>
> file = <filepath>
>
> approach vs the
>
> model_code = <string>
>
> approach?

No. It's exactly the same.

> I realized there was the file option, but opted for the string approach because I preferred to keep as much of
> my code for this analysis in one place as possible.

Ah. My motivation was to cut down cut-and-paste errors.
If you have an entire R script in a file along with
the model, that arrangement makes sense, too.

I tend to run models with different data, though, and
like to keep the following concerns separated:

-- model
-- R code (or other code) to munge data
-- R code to call fit code (including inits)
-- R code to draw plots and analyze
-- one script to run everything for an analysis

The main reason is that I often run the same model over
data from different sources with different shapes and then
run different output analyses over them.

The style in R is to have one big function do everything,
which gives you nice encapsulation but lower flexibility,
whereas my arrangement is more unix-like, which is more
modular, but less encapsulated.

- Bob

[Josef Fruehwald]

Ok, thanks. I think my approach now will be to do some initial aggregation. I'm almost certain it's taking so long because of the volume of data (

79,246 observations). Words are nested within speakers, and the context effects vary between words, so I think what I'll do is summarise mean(y) for words within speakers, giving me 12,844 observations. 

-Joe

[Bob Carpenter]

I was only looking for where the error in parsing/compiling
came from.

You can get a whole lot more efficiency by vectorizing
the call to

speaker_mu[i,j] ~ normal(...);

Allocate a vector for the means, fill it with the mean
values in a loop, then call vectorized normal all at
once. That means replacing this:

for (i in 1:max_speaker) {
for (j in 1:max_context) {

speaker_mu[i,j] ~ normal(Delta[dob[i],sex[i],j]
+ (male_effect * (sex[i]-1))
+ (context_effect * (j-1))
+ intercept,
sigma);
}
speaker_sigma[i] ~ cauchy(0, 2.5) T[0,];
}

with this:


for (i in 1:max_speaker) {
real predix[max_context];
for (j in 1:max_context)
predix[j] <- Delta[dob[i],sex[i],j]

+ (male_effect * (sex[i]-1))
+ (context_effect * (j-1))

+ intercept;
speaker_mu[i] ~ normal(predix,sigma);
}

for (i in 1:I)
speaker_sigma[i] ~ cauchy(0, 2.5);

Note you don't need truncation because sigma's constrained
to (0,inf) and we only need computation to a proportion (and
those cumulatives are expensive!).

You could also put it all in one big array, which is even
faster, but requires loop fiddling. We'll eventually get around
to vectorizing to matrices, but we're doing vectors for
everything first.

And then the same trick with N, which is even more important
if N is large relative to everything else:

for(i in 1:N){
y[i] ~ normal(speaker_mu[speaker[i], context[i]] + word_mu[word[i]], speaker_sigma[speaker[i]]);
}

becomes

real predix2[N];
real spkr_sigma_spkr[N];
for (i in 1:N){
predix2[i] <- speaker_mu[speaker[i], context[i]] + word_mu[word[i]];
spkr_sigma_spkr[i] <- speaker_sigma[speaker[i]]);
}
y ~ normal(predix2,spkr_sigma_spkr);

This is MUCH faster. Basically, you can think of loops and assignments
as basically free in Stan (they're compiled to efficient C++, unlike
R or BUGS loops). What's expensive is each ~ used and each complicated
function used in an assignment.

- Bob

[Josef Fruehwald]

Just so I understand the basics, this is what I understand from your advice:

1) I won't be eliminating any loops, but will be changing what's done inside of a loop.

2) It's less expensive to create a vector and assign values into it through a loop and then do the vectorized sampling, than it is to do the sampling inside of the loop.

Is that right?

-Joe

[Josef Fruehwald]

I'm getting an "Error in stanc(file = "delta_stan.R") : c++ exception (unknown reason)" again. I've only made the change from

  for(i in 1:N){
    y[i] ~ normal(speaker_mu[speaker[i], context[i]] + word_mu[word[i]], speaker_sigma[speaker[i]]);

  }

to

  real predix2[N];

  real speaker_sigma_v[N];    

  for(i in 1:N){

      predix2[i] <- speaker_mu[speaker[i], context[i]]  + word_mu[word[i]];

      speaker_sigma_v[i] <- speaker_sigma[speaker[i]];

        

      //y[i] ~ normal(speaker_mu[speaker[i], context[i]] + word_mu[word[i]], speaker_sigma[speaker[i]]);

  }

    

  y ~ normal(predix2, speaker_sigma_v);

I'd tell you my C++ version, but I'm afraid I'm a bit too amateurish to do that. I'm running OSX 10.8, and I have whatever XCode 4.4.1 installs.

-Joe

[Bob Carpenter]

On 9/21/12 4:44 PM, Josef Fruehwald wrote:
> Just so I understand the basics, this is what I understand from your advice:
>
> 1) I won't be eliminating any loops, but will be changing what's done inside of a loop.
> 2) It's less expensive to create a vector and assign values into it through a loop and then do the vectorized sampling,
> than it is to do the sampling inside of the loop.
>
> Is that right?

Right.

It's not the loops that are expensive, it's the sampling
statements executed inside them. Assignments and vector
allocations (implied by local variable declarations) are
relatively inexpensive (probably not even measurable as
a fraction of total work in cases like yours).

For reasons having to do with how automatic differentation
of the log probability function is done, it's MUCH MUCH
more efficient to call one vectorized sampling statement than
to call the individual sampling statements in a loop.
It saves substantial work in allocation, reduces memory
usage, and greatly reduces the number of virtual function
calls (the big bottleneck for modern CPUs is figuring
out where a program will branch to and a bunch of
very unpredictable virtual function calls is about as bad
as it gets for that).

For what it's worth, you can open the hood on JAGS and
add vectorized distributions there. This provides a pretty
big speedup for similar reasons (they interpret loops and
function calls -- we interpret expressions during auto-dif).

- Bob

[Ben Goodrich]

On Friday, September 21, 2012 5:13:56 PM UTC-4, Josef Fruehwald wrote:

I'd tell you my C++ version, but I'm afraid I'm a bit too amateurish to do that. I'm running OSX 10.8, and I have whatever XCode 4.4.1 installs.

You can go open the Terminal program in the Applications folder and execute

g++ --version

which in your case is most likely 4.2.x . That is the oldest compiler that is known to basically work well with Stan.

Ben

[Jiqiang Guo]

Could you please providing all the relative R code and stan code that generated the error again?  In general, the more information you provide, the easier for us to replicate your error.  From the code you provided once, it seems you have this error, but Bob only got an error like this

LOCATION:  file=input; line=61, column=36

        Delta[i,1,1] <- sum(delta[1:i]);

                                   ^-- here

DIAGNOSTIC(S) FROM PARSER:

Parser expecting: "]"

By the way (to stan-dev), this is an error generated from stanc function in RStan, not the generated c++ code 

cannot be compiled.  To give an idea how it is generated, function stanc in c++ is kind of implemented like this

SEXP stanc(SEXP) {

   try {

     // parse the stan code 

   } catch(std::exception& e) {

      // report error
   } catch(...) {

      Rf_error("c++ exception (unknown reason)")

   }

} 

I cannot think of what issue would trigger this type of error. 

Jiqiang

[Josef Fruehwald]

I've attached the relevant R workspace, Stan model, and R code here. I've been getting the following error a lot as I've been writing my models:

Error in stanc(file = file, model_code = model_code, model_name = model_name,  : 

  c++ exception (unknown reason)

Which is why I started running them through stanc() first. Relevant info my my version of rstan and C++ are included in the R file.

-Joe

[Jiqiang Guo]

Well, I still cannot replicate the error you had.  Instead, I have

> stanc(file = "delt_vector.stan")

TRANSLATING MODEL 'delt_vector' FROM Stan CODE TO C++ CODE NOW.

Error in stanc(file = "delt_vector.stan") : 

  failed to parse Stan model 'delt_vector' and error message provided as:

LOCATION:  file=input; line=116, column=5

    real predix[N];

    ^-- here

DIAGNOSTIC(S) FROM PARSER:

variable "real" does not exist.

Parser expecting: "}"

So I move variable declaration "real predix[N];" and "real sp_sigma_v[N];" to the 

beginning of model block and add semicolon to the end of line 120.  Then at least 

stanc can parse it for me. 

116     real predix[N];

117     real sp_sigma_v[N];

118       for(i in 1:N){

119         predix[i] <- speaker_mu[speaker[i], context[i]] + word_mu[word[i]];

120         sp_sigma_v[i] <- speaker_sigma[speaker[i]]

121     }

122     y ~ normal(predix, sp_sigma_v);

  

For now, the only thing that I can think of is to ask what is your version of Rcpp: the most recent one is Rcpp_0.9.13.

Jiqiang

[Josef Fruehwald]

I've managed to vectorize the data loop, which has sped things up, but still no success with the loop defining speaker_mu. My version of Rcpp is 0.9.13.

Should I generally be expecting informative error messages from stanc() like the ones you've been posting? I've only ever been getting "c++ exception (unknown reason)" which has been  a little frustrating with this largish model, since I'm never sure if I've left out a ; or if I failed to declare a variable, or what.

-Joe

[Jiqiang Guo]

In general, everyone should have the an informative error message if stanc found problem with the model.  I am not sure why you always get 

the 'c++ exception ....' message.  It is a bit weird though there has got to be a reason. 

Could you please try the following R code and let's see what error message you will get?

library(Rcpp)

library(inline)

src <- 'throw std::domain_error("domain error with error message: a fake one");' 

fx <- cxxfunction(body = src, plugin = 'Rcpp', verbose = TRUE)

fx()

The following is what you are supposed to get:

> fx()
Error in fx() : domain error with error message: a fake one

--

Jiqiang

[Josef Fruehwald]

I ran this in a clean R session, and got the unknown reason error. Here's everything:

> sessionInfo()

R version 2.15.1 (2012-06-22)

Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)

locale:

[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:

[1] graphics  grDevices utils     datasets  stats     methods   base     

other attached packages:

[1] ggplot2_0.9.1

loaded via a namespace (and not attached):

 [1] colorspace_1.1-1   dichromat_1.2-4    digest_0.5.2       grid_2.15.1        labeling_0.1       MASS_7.3-18        memoise_0.1       

 [8] munsell_0.3        plyr_1.7.1         proto_0.3-9.2      RColorBrewer_1.0-5 reshape2_1.2.1     scales_0.2.1       stringr_0.6       

> library(Rcpp)

> library(inline)

> 

> src <- 'throw std::domain_error("domain error with error message: a fake one");' 

> fx <- cxxfunction(body = src, plugin = 'Rcpp', verbose = TRUE)

 >> setting environment variables: 

PKG_LIBS =  /Library/Frameworks/R.framework/Versions/2.15/Resources/library/Rcpp/lib/x86_64/libRcpp.a

 >> LinkingTo : Rcpp

CLINK_CPPFLAGS =  -I"/Library/Frameworks/R.framework/Versions/2.15/Resources/library/Rcpp/include" 

 >> Program source :

   1 : 

   2 : // includes from the plugin

   3 : 

   4 : #include <Rcpp.h>

   5 : 

   6 : 

   7 : #ifndef BEGIN_RCPP

   8 : #define BEGIN_RCPP

   9 : #endif

  10 : 

  11 : #ifndef END_RCPP

  12 : #define END_RCPP

  13 : #endif

  14 : 

  15 : using namespace Rcpp;

  16 : 

  17 : 

  18 : // user includes

  19 : 

  20 : 

  21 : // declarations

  22 : extern "C" {

  23 : SEXP file739f60f567f5( ) ;

  24 : }

  25 : 

  26 : // definition

  27 : 

  28 : SEXP file739f60f567f5(  ){

  29 : BEGIN_RCPP

  30 : throw std::domain_error("domain error with error message: a fake one");

  31 : END_RCPP

  32 : }

  33 : 

  34 : 

Compilation argument:

 /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB file739f60f567f5.cpp 2> file739f60f567f5.cpp.err.txt 

g++ -arch x86_64 -I/Library/Frameworks/R.framework/Resources/include -I/Library/Frameworks/R.framework/Resources/include/x86_64 -DNDEBUG  -I/usr/local/include -I"/Library/Frameworks/R.framework/Versions/2.15/Resources/library/Rcpp/include"   -fPIC  -g -O3  -c file739f60f567f5.cpp -o file739f60f567f5.o

g++ -arch x86_64 -dynamiclib -Wl,-headerpad_max_install_names -undefined dynamic_lookup -single_module -multiply_defined suppress -L/usr/local/lib -o file739f60f567f5.so file739f60f567f5.o /Library/Frameworks/R.framework/Versions/2.15/Resources/library/Rcpp/lib/x86_64/libRcpp.a -F/Library/Frameworks/R.framework/.. -framework R -Wl,-framework -Wl,CoreFoundation

> fx()

Error in fx() : c++ exception (unknown reason)

> packageVersion("Rcpp")
[1] ‘0.9.13’

I'm guessing this is some sort of problem either with my installation of Rcpp, or with C++, then.

-Joe

[Jiqiang Guo]

I do not know what goes wrong with your C++ or Rcpp.  You could try the following simple C++ code to rule out it is because of Rcpp.  Then C++ experts out there (I am not) might have some idea.  If it is a problem of Rcpp, you may ask people on the rcpp mail list. 

First copy the following to a file named a.cpp (or whatever you want). 

#include <stdexcept>

#include <iostream>

int main() {

  try {

    throw std::domain_error("a fake error"); 

  } catch (std::exception& e) {

    std::cerr << e.what() << std::endl; 

  } catch (...) {

    std::cerr << "unknow c++ error" << std::endl; 

  } 

} 

Suppose the file name is a.cpp, open a Terminal and change to the directory where you put file a.cpp, compile a.cpp using g++, and execute it as follows. 

sh-3.2$ g++ a.cpp

sh-3.2$ ./a.out 

a fake error

sh-3.2$ 

--

Jiqiang