Difficulty compiling Stan

[Nick Menzies]

Hi guys,

I have having some trouble putting together a disease progression model in Stan. Essentially, the compilation fails despite the debugger thinking everything is fine. 

I have pasted the error message, session info, and compiler info below. I have also attached the full error output as a txt file in case that is useful, and the model I am trying to run (not that I want anyone to wade  though this, but to be comprehensive).

The model simulated the distribution of a population across various compartments, with rates of transition between compartments changing over time (discrete time steps) as a function of the parameters. I was not sure whether such a model could be built and fit in Stan, so have been building it up progressively, but am stymied now. One clue might be this text buried in the full error output (just after the first instance of the term "error"):

'YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY' is not a member of 'Eigen::internal::static_assertion<false>'

Any ideas? Thanks for your help.

Nick

#### ERROR MESSAGE

Error in compileCode(f, code, language = language, verbose = verbose) : 

  Compilation ERROR, function(s)/method(s) not created! In file included from C:/Users/Nick/Documents/R/win-library/3.0/rstan/include//stanlib/eigen_3.1.2/Eigen/Core:276:0,

                 from C:/Users/Nick/Documents/R/win-library/3.0/rstan/include//stanlib/eigen_3.1.2/Eigen/Dense:1,

                 from C:/Users/Nick/Documents/R/win-library/3.0/rstan/include//stansrc/stan/math/matrix/Eigen.hpp:5,

                 from C:/Users/Nick/Documents/R/win-library/3.0/rstan/include//stansrc/stan/math/matrix.hpp:4,

                 from C:/Users/Nick/Documents/R/win-library/3.0/rstan/include//stansrc/stan/io/dump.hpp:12,

                 from C:/Users/Nick/Documents/R/win-library/3.0/rstan/include/rstan/stan_fit.hpp:13,

                 from C:/Users/Nick/Documents/R/win-library/3.0/rstan/include/rstan/rstaninc.hpp:3,

                 from file24906023060.cpp:6:

C:/Users/Nick/Documents/R/win-library/3.0/rstan/include//stanlib/eigen_3.1.2/Eigen/src/Core/Assign.h: In member function 'D

In addition: Warning message:

running command 'C:/PROGRA~1/R/R-30~1.1/bin/x64/R CMD SHLIB file24906023060.cpp 2> file24906023060.cpp.err.txt' had status 1 

##### SESSION INFO

R version 3.0.1 (2013-05-16)

Platform: x86_64-w64-mingw32/x64 (64-bit)

locale:

[1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United States.1252    LC_MONETARY=English_United States.1252 LC_NUMERIC=C                          

[5] LC_TIME=English_United States.1252    

attached base packages:

[1] stats4    stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:

[1] rstan_1.3.0   inline_0.3.11 Rcpp_0.10.3  

loaded via a namespace (and not attached):

[1] codetools_0.2-8 tools_3.0.1  

#### COMPILER

G__~1.EXE (GCC) 4.6.3 20111208 (prerelease)

Copyright (C) 2011 Free Software Foundation, Inc.

This is free software; see the source for copying conditions.  There is NO

warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

[Ben Goodrich]

Anytime something parses but does not compile, it is a bug in Stan. The debug mode is to diagnose run-time errors.

As far as I can tell, the cuts vector should be declared and populated in the transformed data block, rather than the transformed parameters block (and the same for anything else that is fixed). The compilation error happens when you try to divide by 70, and integer division is a bit fragile. I always do double division (i.e. / 70.0).

Ben

[Nick Menzies]

Ben, thanks for the quick response, and the tips -- you are right, it was the integer division that was the problem.

Best - N

--
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/os-sGfUglBo/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.
For more options, visit https://groups.google.com/groups/opt_out.

--

Nick Menzies
nick.m...@gmail.com
404 217 1076

[Bob Carpenter]

On 9/18/13 5:16 PM, Ben Goodrich wrote:
> Anytime something parses but does not compile, it is a bug in Stan. The debug mode is to diagnose run-time errors.

Could someone send me (or the list) the model that parses but
doesn't compile? As Ben says, it's a bug, and we want to fix
these. They're mostly related to int vs. double these days,
and all the ones we're aware of are fixed for the forthcoming
2.0 release.

> As far as I can tell, the cuts vector should be declared and populated in the transformed data block, rather than the
> transformed parameters block (and the same for anything else that is fixed). The compilation error happens when you try
> to divide by 70, and integer division is a bit fragile. I always do double division (i.e. / 70.0).

Integer division rules are simple, and follow programming
language conventions (C, C++ and Java, anyway):

(int/int) returns int, rounds off

(int/real) and (real/int) return real, promoting
int to real before divison

- Bob

[Bob Carpenter]

On 9/19/13 2:56 PM, Bob Carpenter wrote:
>
>
> On 9/18/13 5:16 PM, Ben Goodrich wrote:
>> Anytime something parses but does not compile, it is a bug in Stan. The debug mode is to diagnose run-time errors.
>
> Could someone send me (or the list) the model that parses but
> doesn't compile?

Never mind. I worked back to the e-mail that had
the attachment.

- Bob

[Nick Menzies]

In case it is useful: I was able to solve the problem as Ben suggested, making sure I divided by real rather than int, and shifting these fixed values to the transformed data block. It was the first of these changes that corrected everything.

Of note, I now get the following error message numerous times as Stan iterates through. I have checked through and cannot see how the variable fc_all takes on a negative value (it is a linear combination of things which are defined to be non-negative, unless I am missing a coding error), but was hoping that it will go away once I have gotten to the point where I can run the model with realistic inputs.

Informational Message: The parameter state is about to be Metropolis rejected due to the following underlying, non-fatal (really) issue (and please ignore that what comes next might say 'error'): Invalid value of fc_all: Error in function validate transformed params N4stan5agrad3varE: fc_all is -1.#IND:0, but must be greater than 0

If the problem persists across multiple draws, you might have a problem with an initial state or a gradient somewhere.

 If the problem does not persist, the resulting samples will still be drawn from the posterior.

I tried adding epsilon() , and even 0.5 (i.e. below), to sum that defines fc_all , but no change.

    fc_all <- V1[4]*f0c+V1[5]*fic[1]+V1[6]*fic[2]+V1[7]*fic[3]+V1[8]*fic[4]+

              V1[9]*fic[5]+V1[10]*fic[6]+V1[11]*fic[7]+V1[12]*fic[8]+V1[13]*fic[9]+

              +0.5;    # adding small pos value to see if error code goes away... does not.

Thanks - Nick

- Bob

--
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/os-sGfUglBo/unsubscribe.

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

For more options, visit https://groups.google.com/groups/opt_out.

[Daniel Lee]

On Thu, Sep 19, 2013 at 3:15 PM, Nick Menzies <nick.m...@gmail.com> wrote:

Of note, I now get the following error message numerous times as Stan iterates through. I have checked through and cannot see how the variable fc_all takes on a negative value (it is a linear combination of things which are defined to be non-negative, unless I am missing a coding error), but was hoping that it will go away once I have gotten to the point where I can run the model with realistic inputs.

Informational Message: The parameter state is about to be Metropolis rejected due to the following underlying, non-fatal (really) issue (and please ignore that what comes next might say 'error'): Invalid value of fc_all: Error in function validate transformed params N4stan5agrad3varE: fc_all is -1.#IND:0, but must be greater than 0

If the problem persists across multiple draws, you might have a problem with an initial state or a gradient somewhere.

 If the problem does not persist, the resulting samples will still be drawn from the posterior.

The fc_all value isn't actually a negative number. In Windows, NaN is represented by 1.#IND. See: http://www.johndcook.com/IEEE_exceptions_in_cpp.html

I'm not sure why it thinks it's a negative NaN, but the bottom line is that the computation is hosed. Just FYI, the ":0" printed after the number is the adjoint value of Stan's autodiff variables.

Have you tried the print() function to see the values of:
  V1[4]*f0c+V1[5]*fic[1]+V1[6]*fic[2]+V1[7]*fic[3]+V1[8]*fic[4]+V1[9]*fic[5]+V1[10]*fic[6]+V1[11]*fic[7]+V1[12]*fic[8]+V1[13]*fic[9]
?

I'm having a trouble following the model -- I'd suggest starting with the print and working backwards. My guess is that one of the parameters above isn't initialized before it's being used.

 

I tried adding epsilon() , and even 0.5 (i.e. below), to sum that defines fc_all , but no change.

    fc_all <- V1[4]*f0c+V1[5]*fic[1]+V1[6]*fic[2]+V1[7]*fic[3]+V1[8]*fic[4]+

              V1[9]*fic[5]+V1[10]*fic[6]+V1[11]*fic[7]+V1[12]*fic[8]+V1[13]*fic[9]+

              +0.5;    # adding small pos value to see if error code goes away... does not.

I tried inspecting the code and following around fic and V1, but it looks like a lot of these variables are set and reset multiple times.

Some suggestions:

- Move constants like mortN to the "transformed data" block.

- It looks like V1 and V2 could be transformed data or possibly local variables declared in "transformed parameters." If so, move them. This will help with the clarity and efficiency. (This might be true for some of the other variables also.)

- Start simpler. It's really hard coming in from the outside trying to debug a model of this complexity.

Good luck. If you get past the NaN issue and are still stuck, keep us posted.



Daniel

[Nick Menzies]

Hi Dan, thanks for these comments,

That makes more sense about the value of fc_all (i.e. being NAN).  I will hunt around and see if I can find anything uninitialized, though the model still runs most of the time, it just seems that 1 out of 20 iterations falls apart. This should rule out something being uninitialized, correct?

You are correct that some parameters set set and reset a number of times. The model I am trying to construct is a Markov model, so there is a vector that gets updated many times.

For the variables that get updated many times over, maybe I should be initializing them as local variables inside the loop, and I will give that a try. 

I definitely agree with starting simple -- I have been trying to do so but this is the moment where things have started to break down. Sorry the model looks like a condstruction zone, if I need to send out for further advice it will definitely be a cleaned up verison.

Frankly, my big worry is around the autodiff: this model will do a lot of manipulation within the transformed parameters block, such that the values that will finally be used in the likelihood will be a <very> complex function of the model parameters. Is this something that might boggle the gradient calculations?  This has been my concern all along but things have (seemed to have) been working so far.  I saw Bob's post on Andres'w blog the other day on 'samplers for big science', and though this model will be a teeny-tiny thing compared to some of them, it might be heading in that direction.

Again, thanks for the advice. 

Nick

--

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/os-sGfUglBo/unsubscribe.

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

For more options, visit https://groups.google.com/groups/opt_out.

[Bob Carpenter]

On 9/19/13 6:04 PM, Nick Menzies wrote:
> Hi Dan, thanks for these comments,
>
> That makes more sense about the value of fc_all (i.e. being NAN). I will hunt around and see if I can find anything
> uninitialized, though the model still runs most of the time, it just seems that 1 out of 20 iterations falls apart. This
> should rule out something being uninitialized, correct?

It seems more likely that it would be a divide-by-zero issue,
which can arise due to the "zero" being underflowed.

For instance (not an example from your code),
it looks mathematically like 1 / exp(x) would never be NaN,
but if x < -800, then exp(x) = 0, and 1 / exp(x) == 1 /0 == NaN.

Extreme initializations can put you out in places where
the model underflows, especially if there's a lot of data that
tightly constrains the posterior --- values outside of that can
be ill behaved.

> You are correct that some parameters set set and reset a number of times. The model I am trying to construct is a Markov
> model, so there is a vector that gets updated many times.
>
> For the variables that get updated many times over, maybe I should be initializing them as local variables inside the
> loop, and I will give that a try.

This shouldn't change behavior, but it's usually easier
to visually spot uninitialized variables that way.

> I definitely agree with starting simple -- I have been trying to do so but this is the moment where things have started
> to break down. Sorry the model looks like a condstruction zone, if I need to send out for further advice it will
> definitely be a cleaned up verison.
>
> Frankly, my big worry is around the autodiff: this model will do a lot of manipulation within the transformed parameters
> block, such that the values that will finally be used in the likelihood will be a <very> complex function of the model
> parameters. Is this something that might boggle the gradient calculations?

No, should be fine. By very complex, I was talking about code
that took person-years to write, not something you could put in
a Stan file.

The biggest issue we face is overflow and underflow in gradients.

> This has been my concern all along but
> things have (seemed to have) been working so far. I saw Bob's post on Andres'w blog the other day on 'samplers for big
> science', and though this model will be a teeny-tiny thing compared to some of them, it might be heading in that direction.

In which case, you might want to think about pulling out some of
the calculations into C++. Stan isn't really engineered for building
large code bases --- no namespace support, no stepwise debuggers, very
weak printing and exception control, etc.

- Bob

[Nick Menzies]

Hi Bob - thanks again. Your suspicion could be dead on.  

I had a 'placeholder' likelihood included in the code, to teach myself how it worked and sit there until I got all the other parts of the model set up, but now that other parts of the model are starting to produce realistic results there was a substantial gap between the likelihood and the model estimates, which could have created the situation you describe. I have fixed this now and things work better.

Having said that, I am planning to fit the model to data which will probably constrain some parameters very tightly, so I will keep my eye out for the situation you describe. 

If this becomes a problem, what do people think about this idea:

1.  Fitting the model to a subset of the data (i.e. more diffuse likelihood)

2. Using the point estimates (or maybe random draws) from the resulting posterior as initial values. Not sure whether this training data should then be excluded from the main analysis, but I assume not.

Thanks - Nick

--
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/os-sGfUglBo/unsubscribe.

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

For more options, visit https://groups.google.com/groups/opt_out.

[Bob Carpenter]

On 9/21/13 3:41 PM, Nick Menzies wrote:
> Hi Bob - thanks again. Your suspicion could be dead on.
>
> I had a 'placeholder' likelihood included in the code, to teach myself how it worked and sit there until I got all the
> other parts of the model set up, but now that other parts of the model are starting to produce realistic results there
> was a substantial gap between the likelihood and the model estimates, which could have created the situation you
> describe. I have fixed this now and things work better.
>
> Having said that, I am planning to fit the model to data which will probably constrain some parameters very tightly, so
> I will keep my eye out for the situation you describe.
>
> If this becomes a problem, what do people think about this idea:
> 1. Fitting the model to a subset of the data (i.e. more diffuse likelihood)
> 2. Using the point estimates (or maybe random draws) from the resulting posterior as initial values.

You can even use a direct point estimate as an initial
value, like a posterior mode.

Or you can just take a guess if you know something about
the problem itself.

You have to be careful that you start in an area near
high posterior mass. There's a counterintuitive situation
for the high-dimensional normal where there is very
little probability mass around the posterior mean and mode.

So I'd think random draws would be more reliable.
Or something randomized around the posterior mode or mean.

> Not sure whether
> this training data should then be excluded from the main analysis, but I assume not.

Leave it in. You still want to sample until convergence
and then run long enough that there's no measurable dependence
on the initial point. Initializing at diffuse starting points
helps ensure you go far enough.

So even if you start from a point estimate, it can help to add
noise. You just have to make sure to add noise on the right scale
and in such a way that constraints are satisfied. (We're going to
make doing this automatic in the future, but we're not there yet.)

- Bob

[Nick Menzies]

Hi guys, an update and additional questions.  I managed to get past the underflow problems mentioned below. Part of the solution was getting the model to a point where I was not using placeholder inputs, and the other aspect was to calculate the multinomial coefficients outside of the model. This was mainly to be quicker (only needs to happen once, and doing it in R made it easier to see what was going wrong), but also showed me that I had a problem where I was calculating log(theta^y) as log(theta)*y. It is my understanding that the later is more numerically stable, but I was getting caught by log(0.0)*0.0 --> the log(0) gave me a -Inf, that then became NAN when multiplied by 0.0.

So all good on that front.

I now have a new problem: the model is intended to have a likelihood with three components (data on different aspects of smoking behavior). When I get the first likelihood going everything seems to work fine (well, close to fine), but when I include the other two (or one of the other two), I get a problem similar to that described by Ross Boylan in an earlier post, in that the model compiles, starts to sample, and then memory usage climbs smoothly up to 8GB and then the process fails.  Any ideas?  In case it is useful, I am attaching the model script (StanScriptV30.r. I can also provide the data as well if it is useful to try and replicate the behavior.  I checked through the data for these other likelihoods and things seem fine there. It usually requires killing R to get memory usage back down to normal.

Also, if anything looks like it could be sped up let me know. 

Other useful details: 

> system("g++ --version")

G__~1.EXE (GCC) 4.6.3 20111208 (prerelease)

> sessionInfo()

R version 3.0.1 (2013-05-16)

Platform: x86_64-w64-mingw32/x64 (64-bit)

locale:

[1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United States.1252    LC_MONETARY=English_United States.1252 LC_NUMERIC=C                          

[5] LC_TIME=English_United States.1252    

attached base packages:

[1] stats4    stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:

[1] rstan_1.3.0   inline_0.3.11 Rcpp_0.10.3  

loaded via a namespace (and not attached):

[1] tools_3.0.1

Thanks (again and again...)   Nick

[Bob Carpenter]

On 9/25/13 10:31 AM, Nick Menzies wrote:
> Hi guys, an update and additional questions.

Thanks for the update --- we appreciate it.

> I managed to get past the underflow problems mentioned below. Part of the
> solution was getting the model to a point where I was not using placeholder inputs, and the other aspect was to
> calculate the multinomial coefficients outside of the model. This was mainly to be quicker (only needs to happen once,
> and doing it in R made it easier to see what was going wrong), but also showed me that I had a problem where I was
> calculating log(theta^y) as log(theta)*y. It is my understanding that the later is more numerically stable,

It is. It prevents overflow or underflow in the theta^y term.

> but I was
> getting caught by log(0.0)*0.0 --> the log(0) gave me a -Inf, that then became NAN when multiplied by 0.0.

With theta^y, it's not clear what to do with 0^0. R defines
the value to be 1, so log(0^0) = 1.

You can do one-time calculations in the transformed data block
in Stan. It should be more efficient than R.

To deal with the log(0^0) issue, we have a function

multiply_log(x,y) == x * log(y)

that returns 0 for multiply_log(0,0) for just the reason you mention.
I just realized in looking it up that we don't doc the boundary
behavior.

> So all good on that front.
>
> I now have a new problem: the model is intended to have a likelihood with three components (data on different aspects of
> smoking behavior). When I get the first likelihood going everything seems to work fine (well, close to fine), but when I
> include the other two (or one of the other two), I get a problem similar to that described by Ross Boylan in an earlier
> post, in that the model compiles, starts to sample, and then memory usage climbs smoothly up to 8GB and then the process
> fails. Any ideas? In case it is useful, I am attaching the model script (StanScriptV30.r. I can also provide the data
> as well if it is useful to try and replicate the behavior. I checked through the data for these other likelihoods and

> things /seem/ fine there. It usually requires killing R to get memory usage back down to normal.

If you can send me (ca...@alias-i.com) the data, I'll look into
the memory usage. We definitely want to find any potential memory
leaks, so it's a high priority for us.

It may not be a leak so much as just using a lot of data. With
transformed parameters like:

vector[4] out1[n_states*7000] ; // to accumulate results.

you're going to be writing a whole lot of output, which R needs to
store in each iteration (not to mention causing a lot of memory allocation
in auto-diff).

You could try running from the command-line version, which streams
the output to file (there may be an option to do that in R, too).

> Also, if anything looks like it could be sped up let me know.

First, I'd suggest using more standard bracketing and indentation.
I find your version very hard to scan:

for(j in 1:9) {
ficRR[j] <- exp(BficRR[j]); }

And you don't even need the brackets here.

Your model's a bit dense to understand directly, but if your
comment is reliable (they often aren't, which is why you always
have to read the code, which doesn't lie):

# Calc parameters that dont vary
zBnc_YR[1] <- 1.0;
zBpec_YR[1] <- 1.0;
zBpcfc_YR[1] <- 1.0;
...


then you can save a lot of space and output by defining these
as transformed data rather than transformed parameters.

Also, if you have variables that you want to use for convenience
that you don't need in output, define them as local variables.

In general, when you have repeated subexpressions, such as
(1 - fic[1]) in:

V2[5 ] <- V1[5 ]*(1-fic[1])*(1-ff[2]) + V1[4]*(1-f0c)*ff[1];
V2[6 ] <- V1[5 ]*(1-fic[2])*(1-ff[3]) + V1[4]*(1-fic[1])*ff[2];

it's more efficient to define local variables rather than recomputing.
Each expression takes up space and requires time in auto-diff calc,
and re-using expressions allows the calculations and space to be shared.

When you have big sums, it's more efficient to define an array, assign
elements into the array, then all the sum() function. It saves memory
and time. That is, rewrite:

a <- b + c + d + e;

as:

{
vector[4] temp;
temp[1] <- b;
temp[2] <- c;
temp[3] <- d;
temp[4] <- e;
a <- sum(temp);
}

It seems like it'd be more expensive, but it involves fewer
operations (each of those + operations takes space and time,
whereas sum() is much more efficient). Allocation, loops and
assignment are very fast and minimally sized in R.

And please please use indentation and not something like this:

} ##### # # END OF QTR LOOP
} ##### # # END OF COHORT LOOP
} ##### # # END OF STATE LOOP (not implemented yet)

This:

for(j in 1:n_states) {
Bnc_ST[j] ~ normal( 0.0 , pri_sd ); } // state fixed effect

should be vectorized to:

Bnc_ST ~ normal(0, pri_sd);

Again, it does less memory allocation internally and is faster.

And don't ever do what you did in the same place in the code, which was:

for(j in 1:n_states) {
Bnc_ST[j] ~ normal( 0.0 , pri_sd ); } // state fixed effect
Bnc_AG1 ~ normal( 0.0 , pri_sd );

it makes it look like the Bnc_AG1 is in the loop! Indentation should
always always match the block structure.

Especially for programs this big, I'd recommend keeping the model
in a separate file. I'm constantly having to select and cut models
out of R strings to run them from the command line. And it keeps
the models outside of the code for a particular application.

This could be vectorized

for(j in 1:12) {
rrMort[j] ~ gamma( pri_mort[j,1] , pri_mort[j,2] ) ; }

by reversing order of indexing for pri_mort and writing

rrMort ~ gamma(pri_mort[1], pri_mort[2]);

- Bob

[Nick Menzies]

Hi Bob -- thanks for this. I will send you the data now.  Nick

--
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/os-sGfUglBo/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.
For more options, visit https://groups.google.com/groups/opt_out.

[Sergio Polini]

Il 25/09/2013 19:53, Bob Carpenter ha scritto:
> With theta^y, it's not clear what to do with 0^0. R defines
> the value to be 1, so log(0^0) = 1.

A very nice question. R is in good company:
http://www.math.hmc.edu/funfacts/ffiles/10005.3-5.shtml
http://www.askamathematician.com/2010/12/q-what-does-00-zero-raised-to-the-zeroth-power-equal-why-do-mathematicians-and-high-school-teachers-disagree/
http://en.wikipedia.org/wiki/Exponentiation#Zero_to_the_power_of_zero

In any case, if 0^0 = 1 then log(0^0) = 0 :-)

Sergio

[Bob Carpenter]

Absolutely -- it's all about limits.

It's even trickier when you think about x^y instead
of x^x because x and y can be approaching 0 at different
rates. And trickier still with derivatives, because
x^y is complex if x < 0 and y is not an integer; so you
have to be careful about one-sided vs. two-sided limits.

- Bob

[Daniel Lee]

Nick,

I've spent some time trying to decompose your model, but it's a bit too complicated for me to make a lot of progress (By progress, I mean verifying that the model is what you intended it to be.)

With that said, here's what I've learned:

• There is no leak in the command line version. This is both under v1.3.0 and the current development branch.
  With that said, in some runs, the model is taking up to 12 GB, but it's not leaking. It typically runs at 141 MB on my machine for both v1.3.0 and the current development branch, when it thinks it's actually sampling.
• This doesn't really tell us what's going on with the rstan version, but I don't think there's an actual leak.
• I don't think the model is specified in a way that Stan can handle it. There's a lot of trickiness going on in the transformed parameters block and its showing up in the test gradient mode. From the command line v1.3.0, add the flag "--test_grad" and you can see that sometimes there are large differences between finite diff and the auto diff. Other times all gradients are nan. That's bad news. I was able to consistently get the nan's out by moving the transformed parameter tmp into a local variable, but that didn't really solve the underlying issue. The use of V1 and V2 in the transformed data block creeps me out a bit. It's just setting and resetting those variables too many times for me to be comfortable with it. I'm also having trouble keeping track since they span across different nested loops.
  
• I would really suggest simplifying the model and building up. If this was my model, I'd start by stripping it down to just one parameter and making sure all the data is read in properly. Then I'd add a transformed parameter and make sure that was sampling from the prior properly. I would continue to build up, running the test gradient mode at each step and verifying that the model is sampling properly.
• Is the likelihood right? If I squint, it looks like perhaps you're trying to implement an alternate parameterization of a multinomial distribution, but it's hard for me to see how that's implemented correctly. In Stan, you can drop additive normalizing constants (on the log scale) when they are functions of data, but not when they are functions of parameters or transformed parameters. Have you tried implementing that distribution in a standalone model? That's where I'd start.

At this point, I'm not sure how to help you model as it currently stands. If you can break it down into a simpler model, I might be able to make more progress.

I would approach it in this order: 1) put the likelihood in a model alone to make sure its sampling from the correct distribution, 2) make sure temporary variables are scoped as tight as possible to decouple the transformed parameters block, 3) build up the model starting with a parameter at a time.

If you want help breaking down the model, we can do that on the list. I'm willing to help, I just think I'm stuck with the current incarnation of the model.

Daniel

PS. Try out the Stan emacs mode. https://github.com/stan-dev/stan-mode. It handles indentation and code highlighting for Stan models.

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

[Nick Menzies]

Hi Daniel, Bob,

Thank you both for your detailed comments. I am glad to hear a memory leak is unlikely.

From my end:

• Thanks for the suggestions about areas where code could be written more efficiently, I will add them in.
  
• Similarly, thanks for the formatting tips, sorry it was hard to go through, I will try to iron these out.
  
• For the likelihood, your squinting is correct, this is related to a multinomial, but my data are not all integer due to survey adjustments, so I could not use the built-in function. I will look at it again.
  
• Regarding V1 and V2 being reset so many times: is the concern that this might throw out the autodiff, or make for high memory usage?  It might be useful for me to explain what I am trying to do with the model, as that may point to better ways to approach the task: the model follows cohorts of individuals over time, describing changes in smoking behavior and subsequent mortality. In the intended version of the model there would be 5200 cohorts (52 geographical regions x 100 age categories in the start year). Each cohort is followed for 70 quarters (i.e. 17.5 years), and in any given quarter V1 and V2 represent the distribution of a cohort across smoking behavior categories. The model parameters, once transformed in a couple of ways, determine the probabilities of individuals in one category moving to another category. Most transition parameters get updated each quarter, and the distribution across categories always does, hence the need to reset the values of V1 and V2 so often. At the moment I am only working with 1 state, but that still implies resetting V1 and V2 about 100*70*2 = 14,000 times. In the innermost loop I loop through  1:70 quarters for 1 cohort, the next loop out will go through 1:100 starting ages, and the final loop (once implemented) would go through 1:52 geographic regions.  The likelihood will compare the distributions across smoking behavior categories to empirical data on the same, hence the need to record them all.  Would it change anything to define V1 and V2 as local variables? At the moment I use V1 and V2 to perform all the intermediate calculations then save the results for each quarter/cohort in an array (e.g. out1). Would it be better to directly update the values of the array?  In the end the substantive question I am trying to answer involves relating cigarette taxes to smoking behavior and subsequent mortality, and the identification strategy relies on inter-state comparisons, so the model needs to get bigger still to fulfill its purpose (i.e. build in those other likelihoods and add more states to the outer loop). I stress this issue as it is useful for me to know if I am going to hit barriers that can be fixed by improving the way I code things. Given, from both of your gracious suggestions it sounds like the coding could improve a lot, so I have hope:) 
• Any other ideas why the memory spikes up to 12GB, if it is normally using a lot less memory?  With more testing on my end I was able to generate the same kinds of problems by increasing the number of states (i.e. the outer loop went through 1:5 rather than just 1:1). 

Thanks again - Nick

--

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/os-sGfUglBo/unsubscribe.

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

For more options, visit https://groups.google.com/groups/opt_out.

[Bob Carpenter]

On 9/26/13 10:26 AM, Nick Menzies wrote:
> Hi Daniel, Bob,
>
> Thank you both for your detailed comments. I am glad to hear a memory leak is unlikely.
>
> From my end:
>

> * Thanks for the suggestions about areas where code could be written more efficiently, I will add them in.
> * Similarly, thanks for the formatting tips, sorry it was hard to go through, I will try to iron these out.
> * For the likelihood, your squinting is correct, this is related to a multinomial, but my data are not all integer due

> to survey adjustments, so I could not use the built-in function. I will look at it again.

> * Regarding V1 and V2 being reset so many times: is the concern that this might throw out the autodiff, or make for
> high memory usage?

It won't have any effect on auto-diff or memory usage.

It just that it makes the code harder to understand.

> It might be useful for me to explain what I am trying to do with the model, as that may point to
> better ways to approach the task: the model follows cohorts of individuals over time, describing changes in smoking
> behavior and subsequent mortality. In the intended version of the model there would be 5200 cohorts (52 geographical
> regions x 100 age categories in the start year). Each cohort is followed for 70 quarters (i.e. 17.5 years), and in
> any given quarter V1 and V2 represent the distribution of a cohort across smoking behavior categories. The model
> parameters, once transformed in a couple of ways, determine the probabilities of individuals in one category moving
> to another category. Most transition parameters get updated each quarter, and the distribution across categories
> always does, hence the need to reset the values of V1 and V2 so often.


> At the moment I am only working with 1 state,
> but that still implies resetting V1 and V2 about 100*70*2 = 14,000 times. In the innermost loop I loop through 1:70
> quarters for 1 cohort, the next loop out will go through 1:100 starting ages, and the final loop (once implemented)
> would go through 1:52 geographic regions. The likelihood will compare the distributions across smoking behavior
> categories to empirical data on the same, hence the need to record them all. Would it change anything to define V1
> and V2 as local variables?

I think it'd make the code much clearer. You usually want
to declare something near where it's used so that it's scope
is clear and you can relate its use back to its definition.

> At the moment I use V1 and V2 to perform all the intermediate calculations then save the
> results for each quarter/cohort in an array (e.g. out1). Would it be better to directly update the values of the
> array? In the end the substantive question I am trying to answer involves relating cigarette taxes to smoking
> behavior and subsequent mortality, and the identification strategy relies on inter-state comparisons, so the model
> needs to get bigger still to fulfill its purpose (i.e. build in those other likelihoods and add more states to the
> outer loop). I stress this issue as it is useful for me to know if I am going to hit barriers that can be fixed by
> improving the way I code things. Given, from both of your gracious suggestions it sounds like the coding could
> improve a lot, so I have hope:)

> * Any other ideas why the memory spikes up to 12GB, if it is normally using a lot less memory?

What do you mean by normally using a lot less memory?

The only difference in memory between runs in Stan should
be how many leapfrog steps get taken. With a lot of parameters,
this can make a big difference and may be related to the
instability Daniel was seeing in numerical evaluations.

Some of these problems we can't really fix at a deep level
because Stan lets you write arbitrary code and arbitrary code
with lots of numerical operations (like you have) can be
badly behaved.

So I'd strongly second Daniel's suggestion to move up
to a full model gradually. He has a number of excellent
suggestions. (We practice what we preach here! Daniel and
I have been working on some complicated models over the past
couple of months and always start with the simplest model
we can, even if it doesn't make sense for the final problem.
It's easier to build up from things that work than it is to
diagnose something large that doesn't work.

It's usually good modeling practice to do this anyway because
it lets you compare the results of simpler models to more
complex ones. (And no, I don't want to get dragged into
a technical discussion of model comparison --- I'll leave
that to Andrew on his blog and in the 3rd Edition of
Bayesian Data Analysis.)

> With more testing on
> my end I was able to generate the same kinds of problems by increasing the number of states (i.e. the outer loop
> went through 1:5 rather than just 1:1).

What kinds of problems?

- Bob

> On Wed, Sep 25, 2013 at 10:59 PM, Daniel Lee <bea...@alum.mit.edu <mailto:bea...@alum.mit.edu>> wrote:
>
> Nick,
>
> I've spent some time trying to decompose your model, but it's a bit too complicated for me to make a lot of progress
> (By progress, I mean verifying that the model is what you intended it to be.)
>
> With that said, here's what I've learned:
>

> * There is no leak in the command line version. This is both under v1.3.0 and the current development branch.

> With that said, in some runs, the model is taking up to 12 GB, but it's not leaking. It typically runs at 141 MB
> on my machine for both v1.3.0 and the current development branch, when it thinks it's actually sampling.

> * This doesn't really tell us what's going on with the rstan version, but I don't think there's an actual leak.
> * I don't think the model is specified in a way that Stan can handle it. There's a lot of trickiness going on in

> the transformed parameters block and its showing up in the test gradient mode. From the command line v1.3.0, add
> the flag "--test_grad" and you can see that sometimes there are large differences between finite diff and the
> auto diff. Other times all gradients are nan. That's bad news. I was able to consistently get the nan's out by
> moving the transformed parameter tmp into a local variable, but that didn't really solve the underlying issue.
> The use of V1 and V2 in the transformed data block creeps me out a bit. It's just setting and resetting those
> variables too many times for me to be comfortable with it. I'm also having trouble keeping track since they span
> across different nested loops.

> * I would really suggest simplifying the model and building up. If this was my model, I'd start by stripping it

> down to just one parameter and making sure all the data is read in properly. Then I'd add a transformed
> parameter and make sure that was sampling from the prior properly. I would continue to build up, running the
> test gradient mode at each step and verifying that the model is sampling properly.

> * Is the likelihood right? If I squint, it looks like perhaps you're trying to implement an alternate

> http://www.math.hmc.edu/__funfacts/ffiles/10005.3-5.__shtml
> <http://www.math.hmc.edu/funfacts/ffiles/10005.3-5.shtml>
> http://www.askamathematician.__com/2010/12/q-what-does-00-__zero-raised-to-the-zeroth-__power-equal-why-do-__mathematicians-and-high-__school-teachers-disagree/
> <http://www.askamathematician.com/2010/12/q-what-does-00-zero-raised-to-the-zeroth-power-equal-why-do-mathematicians-and-high-school-teachers-disagree/>
>
> http://en.wikipedia.org/wiki/__Exponentiation#Zero_to_the___power_of_zero

> <http://en.wikipedia.org/wiki/Exponentiation#Zero_to_the_power_of_zero>
>
> In any case, if 0^0 = 1 then log(0^0) = 0 :-)
>
> Sergio
>
>
> --
> 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 <mailto:stan-users%2Bunsu...@googlegroups.com>.
>
> For more options, visit https://groups.google.com/__groups/opt_out <https://groups.google.com/groups/opt_out>.

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

> <mailto:stan-users%2Bunsu...@googlegroups.com>.

> For more options, visit https://groups.google.com/groups/opt_out.
>
>
>
>
> --
>
> Nick Menzies

> nick.m...@gmail.com <mailto:nick.m...@gmail.com>
> 404 217 1076 <tel:404%20217%201076>

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

[Ross Boylan]

Mostly I wanted to ask about the emacs mode, but I also have a few
comments on the memory issue.

As far as I could tell the code was not invoking the
categorical(softmax()) function that seemed to be related to the
memory leak in my case. However, it is doing a multinomial, and so
it's possible that it's hitting the same problematic code.

Since the original poster mentioned bailing out at 8G and the report
here is it worked with 12GB, maybe there isn't actually a leak in this
case. Alternatively, it suggests the problem is R specific.

On Wed, Sep 25, 2013 at 10:59:29PM -0400, Daniel Lee wrote:
> Nick,
>
> I've spent some time trying to decompose your model, but it's a bit too
> complicated for me to make a lot of progress (By progress, I mean verifying
> that the model is what you intended it to be.)
>
> With that said, here's what I've learned:
>

> - There is no leak in the command line version. This is both under

> v1.3.0 and the current development branch.
> With that said, in some runs, the model is taking up to 12 GB, but it's
> not leaking. It typically runs at 141 MB on my machine for both v1.3.0 and
> the current development branch, when it thinks it's actually sampling.

> - This doesn't really tell us what's going on with the rstan version,

> but I don't think there's an actual leak.

[snip]

>
> PS. Try out the Stan emacs mode. https://github.com/stan-dev/stan-mode. It
> handles indentation and code highlighting for Stan models.

Does that work inside of r-mode? In R typically one has
stanmodel <- '
stan code goes here
'

So the code is inside of R inside of a quotation.

I'd be suprised if it worked, but I thought I'd ask. There are other
mode-switching issues with R when it is combined with various
documentation systems; I'm not sure if those have been integrated or
not. If they have, it might suggest a way to get stan-mode working in
R.

Ross

[Nick Menzies]

Hi Daniel, Bob,

Thanks again for the long list of great suggestions. I have implemented them through-out the code and in general they are making things a lot faster, as well as hopefully making the code more understandable.

One aspect I wanted to check specifically: you had noted that I was generating a number of large arrays (e.g. out1) which might be using up a lot of memory. The reason I had out1 was to store a lot of values in the transformed parameters which would then be evaluated against the likelihood under the model block. I have tried to get around the need for out1 by evaluating the likelihood immediately each time I generate new simulation results (i.e. each iteration of the innermost loop). As lp__ cannot be updated outside the model block (correct?) I just use a new parameter to accumulate all the log-likelihoods, then add this to lp__ in the model block.  Acknowledging this does not follow the normal process, do you see any potential problems with this?

Overall the model seems to turn over about twice as fast as before (about 4 iterations per second). However, this does not seem to have fixed the problems mentioned earlier. I seem to encounter one of two general behaviors when I run the model with a low value for iter (e.g. 100) and no thinning :

a. The model adapts quickly and then continues to sample for a long time. This seems to involve many more iterations than 100 (I print lp__ to the console once per iteration, so just count how this cycles through), and lp_ seems to gradually increase. I assume this is just the model working as it should?  I (naively) thought setting iter to 100 might mean that the model only cycles trough 100 times, but it seems this is not the case, it just seem to keep going.

b. Other times the model compiles, everything looks good, it samples once or a couple of times, then memory usage goes up to 8GB (max of my computer) and things freeze. Sometimes this will be be with the same model that worked moments ago (i..e the situation a described above). When I try to expand the model in one of the ways I had planned to (i.e. add additional components to the likelihood, or run the model for multiple states) it seems to crash every time. 

Any thoughts welcome. My plan from here is to simplify the model down and see if that leads to revelations about what is wrong. Per Dan's advice, I have run a super simple model just to make sure the likelihood works, and it seems to be.  Also, Bob had an earlier suggestion to code it up in C++ and I might do this as well (with R and rcpp anyway) so I can run in with an optimization algorithm to generate inits that are closer to the posterior max, in case that helps.

Best - Nick

[Bob Carpenter]

On 9/26/13 10:40 PM, Nick Menzies wrote:
> Hi Daniel, Bob,
>

> Thanks again for the long list of great suggestions. I have implemented them through-out the code and in general they
> are making things a lot faster, as well as hopefully making the code more understandable.
>
> One aspect I wanted to check specifically: you had noted that I was generating a number of large arrays (e.g. out1)
> which might be using up a lot of memory.

The saved samples specifically. Not the arrays themselves
in the model execution.

> The reason I had out1 was to store a lot of values in the transformed
> parameters which would then be evaluated against the likelihood under the model block.

You can make them local variables in the model block if you
don't need to save them.

> I have tried to get around the
> need for out1 by evaluating the likelihood immediately each time I generate new simulation results (i.e. each iteration
> of the innermost loop). As lp__ cannot be updated outside the model block (correct?)

Also in transformed parameters, though I'd recommend only
including the Jacobian adjustments there.

> I just use a new parameter to
> accumulate all the log-likelihoods, then add this to lp__ in the model block. Acknowledging this does not follow the
> normal process, do you see any potential problems with this?

That'll work. But you can just move everything into the
model block as a local variable.

> Overall the model seems to turn over about twice as fast as before (about 4 iterations per second). However, this does
> not seem to have fixed the problems mentioned earlier. I seem to encounter one of two general behaviors when I run the
> model with a low value for iter (e.g. 100) and no thinning :
>
> a. The model adapts quickly and then continues to sample for a long time. This seems to involve many more iterations
> than 100 (I print lp__ to the console once per iteration, so just count how this cycles through), and lp_ seems to
> gradually increase. I assume this is just the model working as it should?

From random initialization, it should increase until convergence, then start
to bounce around. From an MLE initialization, you'd get the opposite, a
decrease until convergence, then bouncing.

> I (naively) thought setting iter to 100 might
> mean that the model only cycles trough 100 times, but it seems this is not the case, it just seem to keep going.

You can set the total iterations. NUTS controls the
number of leapfrog steps, though, unless you specify
the number of leapfrog steps explicitly, which we don't
recommend.

> b. Other times the model compiles, everything looks good, it samples once or a couple of times, then memory usage goes
> up to 8GB (max of my computer) and things freeze. Sometimes this will be be with the same model that worked moments ago
> (i..e the situation a described above). When I try to expand the model in one of the ways I had planned to (i.e. add
> additional components to the likelihood, or run the model for multiple states) it seems to crash every time.

Can you try setting the max tree depth to something lower
and see if you still get memory spikes? That'll control
the number of leapfrog steps per iteration.

Could I ask again which version of Stan or RStan you're using?

> Any thoughts welcome. My plan from here is to simplify the model down and see if that leads to revelations about what is
> wrong. Per Dan's advice, I have run a super simple model just to make sure the likelihood works, and it seems to be.
> Also, Bob had an earlier suggestion to code it up in C++ and I might do this as well (with R and rcpp anyway) so I can
> run in with an optimization algorithm to generate inits that are closer to the posterior max, in case that helps.

I didn't mean to suggest coding in C++! I may have suggested
running from the command line, because that doesn't store all
the samples --- it streams them out to disk as it goes. If you
can figure out good inits, that can help. You definitely don't
want to be starting too far out in the tail of the posterior.

- Bob

[Nick Menzies]

Hi Bob

 

Can you try setting the max tree depth to something lower
and see if you still get memory spikes?  That'll control
the number of leapfrog steps per iteration.

I will give it a go.   

 

Could I ask again which version of Stan or RStan you're using?

rstan_1.3.0 

Also:  

R version 3.0.1 (2013-05-16)

Platform: x86_64-w64-mingw32/x64 (64-bit)

Rtools: tools_3.0.1

Compiler: (GCC) 4.6.3

 

[Nick Menzies]

Regarding C++, the earlier suggestion had been to..

.. you might want to think about pulling out some of 

the calculations into C++.  Stan isn't really engineered for building 
large code bases --- no namespace support, no stepwise debuggers, very 
weak printing and exception control, etc. 

I am confused now.. is is possible to have some part of the model coded separately in C++ with the results passed back to Stan?  I assume not, but thought I should check.  If everything is running as intended, I see my possible problems as:

1. Large parameter space: if I ever get the full model working it will have ~250 parameters, which is a lot, however other users seem to have successfully dealt with problems with many more parameters, so I am assuming this okay for now.

2. Lots of data: in the full model I would be evaluating 4.9 x10^5 likelihoods. I assume this will just make iterations run more slowly, but otherwise things will still work (it seemed to when I just ran a toy model to estimate the overall mean).

3. Lots of action in the transformed parameters block: the relationship between model parameters and the outputs that get entered in the likelihood is very complicated, and takes a while to compute. If this is what is causing things to crash, and if this might be solved by somehow pulling these calculations out separately in C++, I think its worth giving it a try. 

My poorly-informed guess, as the model seems to run sometimes but not others, is that things are going awry when the sampler gets too far into the tail of the posterior. Does this sound possible?

N

[Bob Carpenter]

On 9/27/13 9:10 AM, Nick Menzies wrote:
> Regarding C++, the earlier suggestion had been to..
>

> /.. you might want to think about pulling out some of

> the calculations into C++. Stan isn't really engineered for building
> large code bases --- no namespace support, no stepwise debuggers, very

> weak printing and exception control, etc. /

My bad --- I lost the context there. Our e-mail volume is
starting to overwhelm my memory.

If the problem's really hard and requires involved code in submodules,
then you probably should do it in C++. It's what we're doing
right now for pharmacokinetic models requiring diff eq solvers.
We use C++ for the solvers. Of course, our plan's to get this all
into the Stan language eventually, but for now, the only option is
to write custom C++ functions based on solving diff-eqs then plug
them into Stan.

> I am confused now.. is is possible to have some part of the model coded separately in C++ with the results passed back
> to Stan?

Yes, you can do that, but for now, it requires recompiling
Stan (that is, we don't have a plugin architecture going yet,
but we're thinking about it).

Basically, you can define functions the same way the functions
in stan/agrad/rev are coded up. They need to support auto-diff
if they have inputs that are parameters. This can be as easy
as just supplying template inputs.

> I assume not, but thought I should check. If everything is running as intended, I see my possible problems as:

> *1. Large parameter space: *if I ever get the full model working it will have ~250 parameters, which is a lot, however

> other users seem to have successfully dealt with problems with many more parameters, so I am assuming this okay for now.

We regularly deal with thousands of parameters.

> *2. Lots of data:* in the full model I would be evaluating 4.9 x10^5 likelihoods. I assume this will just make

> iterations run more slowly, but otherwise things will still work (it seemed to when I just ran a toy model to estimate
> the overall mean).

Correct. You mean 50K data points contributing to the total likelihood?
If you need each likelihood to be saved (say for doing something like WAIC),
then it won't be a computational problem, just an I/O or memory size
problem if you're caching the results.

> *3. Lots of action in the transformed parameters block*: the relationship between model parameters and the outputs that

> get entered in the likelihood is very complicated, and takes a while to compute. If this is what is causing things to
> crash, and if this might be solved by somehow pulling these calculations out separately in C++, I think its worth giving
> it a try.

That may be.

> My poorly-informed guess, as the model seems to run sometimes but not others, is that things are going awry when the
> sampler gets too far into the tail of the posterior. Does this sound possible?

Yes. Have you tried initializing with something more
reasonable. We have an init=0 option that will set all
the init values to 0 in the unconstrained space (simplexes
are uniform, positive-constrained variables initialized to 1,
unbounded parameters initialized to 0, bounded parameters
initialized to middle of the bound).

- Bob

[Nick Menzies]

Sorry for sending fragmented updates, but I have now had a chance to run the model with max_treedepth set to a lower value. With max_treedepth set to 5 I am able to run the model with ~10 times as much data (with 10 states rather than just 1), and I was even able to get it working with the full set of 52 states for a while (i.e. it worked through the first chain nicely then choked while getting started on the second chain).

So what to do?  From earlier posts on the subject of max_treedepth I see that the default is a bit arbitrary, but small values mean that the benefits of HMC are not really being harnessed (if I am interpreting Mike's comment on this thread correctly). I can probably answer my research question with a smaller sample (10 states instead of 52), so I might have something that works, but I have a feeling (just from the earlier posts on the subject) that a max_treedepth of 5 may just be silly.

N

[Daniel Lee]

No worries about the fragmented updates. Just make sure to give us enough context so we can follow along.

Before getting too deep, would you mind sending along the simple model you had verifying your likelihood function works? And, when you run the gradient tests, are the gradients now well behaved?

Daniel

[Nick Menzies]

Hi Bob, thanks for the responses..

If the problem's really hard and requires involved code in submodules,
then you probably should do it in C++.  It's what we're doing
right now for pharmacokinetic models requiring diff eq solvers.
We use C++ for the solvers.  Of course, our plan's to get this all
into the Stan language eventually, but for now, the only option is
to write custom C++ functions based on solving diff-eqs then plug
them into Stan.

From this description, I think it is best for me to continue with the current approach!  

 

> *2. Lots of data:* in the full model I would be evaluating 4.9 x10^5 likelihoods. I assume this will just make 
> iterations run more slowly, but otherwise things will still work (it seemed to when I just ran a toy model to estimate
> the overall mean).

Correct.  You mean 50K data points contributing to the total likelihood?
If you need each likelihood to be saved (say for doing something like WAIC),
then it won't be a computational problem, just an I/O or memory size
problem if you're caching the results.

Sorry:  500,000 multinomial likelihoods, with 4 categories each, so ~2M individual data points in total. 

 

> My poorly-informed guess, as the model seems to run sometimes but not others, is that things are going awry when the 
> sampler gets too far into the tail of the posterior. Does this sound possible?

Yes.  Have you tried initializing with something more
reasonable.  We have an init=0 option that will set all
the init values to 0 in the unconstrained space (simplexes
are uniform, positive-constrained variables initialized to 1,
unbounded parameters initialized to 0, bounded parameters
initialized to middle of the bound).

Sounds like better inits is my next step, as well as deciding on an approach to max_treedepth, as mentioned on the previous email.

I hope it goes without saying, but is also worth writing out: thanks for your (and Dan's and Ben's, and others...)  continuing help on this. I really appreciate it.

Nick

 

[Bob Carpenter]

> On Fri, Sep 27, 2013 at 11:35 AM, Nick Menzies <nick.m...@gmail.com <mailto:nick.m...@gmail.com>> wrote:
>
> Sorry for sending fragmented updates, but I have now had a chance to run the model with max_treedepth set to a lower
> value. With max_treedepth set to 5 I am able to run the model with ~10 times as much data (with 10 states rather
> than just 1), and I was even able to get it working with the full set of 52 states for a while (i.e. it worked
> through the first chain nicely then choked while getting started on the second chain).

> So what to do? From earlier posts on the subject of max_treedepth I see that the default is a bit arbitrary, but
> small values mean that the benefits of HMC are not really being harnessed (if I am interpreting Mike's comment on
> this

> <https://groups.google.com/forum/#!searchin/stan-users/max_treedepth/stan-users/3-yyH1q4MGs/sZyjT-g5AKwJ> thread

> correctly). I can probably answer my research question with a smaller sample (10 states instead of 52), so I might
> have something that works, but I have a feeling (just from the earlier posts on the subject) that a max_treedepth of
> 5 may just be silly.

The question's really whether you get the same answer or
not. We haven't investigated this issue a lot ourselves
on NUTS.

Setting a max treedepth of 5 is still sound asymptotically (so is depth
0, which just reduces to Metropolis),
but it'll wind up being more of a random walk than with a deeper bound.

Generally, Matt and Michael (the ones who know the most about
this, along with Marcus) have been suggesting that you really
want to let NUTS go as far as it wants to, because that'll be
most efficient. Andrew, on the other hand, has argued for
lower upper limits. So depth 10 is something of a compromise.

The real question is how to verify that you're getting the
right answer. The danger's that you can get trapped in some
small region of the posterior and it can look like it's
converged (good R-hat and effective sample size estimates), but
it really hasn't. So if you can run with a deeper treedepth
and get the same answer, running with a shallower one should
be OK. Similarly, Michael's been suggesting all along to run
with a fixed lower step size than the one that adaptation finds
and make sure that provides the same answer.

- Bob

[Ross Boylan]

On Fri, Sep 27, 2013 at 01:33:45PM -0400, Bob Carpenter wrote:
...

> Setting a max treedepth of 5 is still sound asymptotically (so is depth
> 0, which just reduces to Metropolis),
> but it'll wind up being more of a random walk than with a deeper bound.
>
> Generally, Matt and Michael (the ones who know the most about
> this, along with Marcus) have been suggesting that you really
> want to let NUTS go as far as it wants to, because that'll be
> most efficient. Andrew, on the other hand, has argued for
> lower upper limits. So depth 10 is something of a compromise.

Another way of interpreting this is that if estimation hits a depth of
10 it is a clue to reconsider your model. My initial efforts hit the
bound. Although the usual advice is to simplify, adding a parameter
was what solved the problem in my case. I had fixed, diffuse
priors, and the fix was to make the dispersion of the priors a
hierarchical parameter.

Ross

[Nick Menzies]

Interesting thought. I have state and year fixed effects which I could change to random effects -- my only reason for the fixed effects was that the results of this analysis could end up playing to an econ audience (the research question is about the causal effect of cigarette excise taxes on various aspects of smoking behavior), and there appears to be a strong bias towards fixed effects to avoid any sniff of bias about the identification strategy. However, if I can get it to work with random effects that is a small battle to fight. 

N

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

For more options, visit https://groups.google.com/groups/opt_out.

[Nick Menzies]

Hi Daniel,

Sure thing -- simple model attached. It is a quick hack to make sure that the mean of the posterior conforms to the mean in the data, so forgive the unconventional prior. 

For the gradient tests -- sorry I didnt send those earlier, only just realized I could access these using rStan (not sure how to use the command line version, so had not moved on that yet). When I run stan with test_grad=TRUE, the gradients produced have no NaNs, though am unsure what represents good behavior. When I compare the finite difference result with the autodiff result the differences are pretty small except for some parameters right down the bottom (124-136). I am not sure how the parameters are indexed, but I assume these parameters at the bottom are the priors for the mortality rate ratios, which are gamma distributions with very strong priors. I am not sure why these would behave differently, but they are not the focus of inference, and prior should really equal posterior for these parameters. However, the inits for these parameters may be way off (given the strong prior) , so that may be something.

I have also attached this test-gradient output as a txt file as you probably have a better-calibrated understanding of whether these values are similar enough. 

N

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

For more options, visit https://groups.google.com/groups/opt_out.

[Ross Boylan]

On Fri, 2013-09-27 at 14:06 -0400, Nick Menzies wrote:
> Interesting thought. I have state and year fixed effects which I could
> change to random effects -- my only reason for the fixed effects was
> that the results of this analysis could end up playing to an econ
> audience (the research question is about the causal effect of
> cigarette excise taxes on various aspects of smoking behavior), and
> there appears to be a strong bias towards fixed effects to avoid any
> sniff of bias about the identification strategy. However, if I can get
> it to work with random effects that is a small battle to fight.
>
>
> N
>

Just to clarify, the change I made was not from fixed to random effects,
though it is perhaps analogous one level up. A simplified version is
the change from
sigma ~ cauchy(0, 5)
beta ~ normal(0, 5)
y ~ normal(beta*x, sigma)


to

sigma ~ cauchy(0, 5)
sigma2 ~ cauchy(0, 5)
beta ~ normal(0, sigma2)
y ~ normal(beta*x, sigma)

The estimated sigma2 was around 0.1


Ross

[Bob Carpenter]

I found exactly that same behavior with an IRT model I was fitting.

Trying to write a Bayesian model to approximate the usual
MLE model required some fiddling with hand-coded priors to identify
the model. Putting a hyperprior on the items and abilities (and
adding an intercept so that items and abilities could be centered
around 0) solved the problem right away.

I should write this up in the manual.

The simplification part is more just to understand the model
by building up from working pieces until you hit something
that doesn't work.

- Bob

[Bob Carpenter]

On 9/27/13 2:22 PM, Ross Boylan wrote:
> On Fri, 2013-09-27 at 14:06 -0400, Nick Menzies wrote:
>> Interesting thought. I have state and year fixed effects which I could
>> change to random effects -- my only reason for the fixed effects was
>> that the results of this analysis could end up playing to an econ
>> audience (the research question is about the causal effect of
>> cigarette excise taxes on various aspects of smoking behavior), and
>> there appears to be a strong bias towards fixed effects to avoid any
>> sniff of bias about the identification strategy. However, if I can get
>> it to work with random effects that is a small battle to fight.
>>
>>
>> N
>>
> Just to clarify, the change I made was not from fixed to random effects,
> though it is perhaps analogous one level up. A simplified version is
> the change from
> sigma ~ cauchy(0, 5)
> beta ~ normal(0, 5)
> y ~ normal(beta*x, sigma)
>
>
> to
>
> sigma ~ cauchy(0, 5)
> sigma2 ~ cauchy(0, 5)
> beta ~ normal(0, sigma2)
> y ~ normal(beta*x, sigma)
>
> The estimated sigma2 was around 0.1

Thanks for sharing the particular example. It's very much
like what I did for the ability parameters in my IRT model.

I'm very curious as to what happens to posterior inference
in these cases to the likelihood for y. My guess is that
with fat fixed priors (scale = 5), the inferences for y won't change
much compared to a tighter prior of scale=0.1. If they did,
I'd think beta would be more constrained by the data.

Although the posterior for beta will look very different in terms of
variance, I haven't found the mean changes much. But you have
to run a lot more iterations to get enough effective samples because
the beta is so unconstrained. I found with the IRT model
that if I used a very weak prior and ran for 100K iterations,
I'd get the same posterior mean for an item difficulty parameter
as running 1K iterations with a hierarchical prior.

I want a name for a model like this --- where if you run long enough
you get the right posterior mean, but where it's not so well identified
that this happens quickly. With a truly non-identified model, it
doesn't matter how long you run, you won't get convergence. And if
it's not identified as to scale, Stan will usually diagnose this and
halt with a warning that the posterior isn't proper because it'll try
to run the parameter off to infinity.

- Bob

[Michael Betancourt]

Right, this is what I think was Andrew's motivation behind having an upper bound. Theoretically you take full advantage of
HMC when you run until you turn around -- anything less drops your performance, especially in high dimensions. But if
it takes more than 2^{10} steps to cross the posterior your model is crazy and the folk theorem comes into play. Usually some
care in the choice of priors (and make sure all constraints are explicitly taken into account in your data types) makes a huge difference.

[Daniel Lee]

Nick,

Thanks for attaching the simple model. Unfortunately, I can't get it to work. I can't sample B1, B2, B3, B4 properly. I didn't want to spend too much time debugging what you have. Instead, I've created a model that works.

My guess is that you really want the Dirichlet distribution instead of a multinomial -- that's from seeing the data and realizing your data are continuous. I'll write out both models, in case you have integer data. I've renamed the variables to something more general.

Multinomial (discrete data)
------------------------------------

data {
  int <lower=0> N;
  int ns[N,4];
}
parameters {
  simplex[4] B;
}
model {
  for (n in 1:N) {
    ns[n] ~ multinomial(B);
  }
}

Dirichlet (continuous data)
------------------------------------
data {
  int <lower=0> N;
  vector[4] ns[N];
}
parameters {
  simplex[4] B;
}
model {
  for (n in 1:N) {
    B ~ dirichlet(ns[n]);
  }
}

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

First, notice the difference in data declaration. For multinomial, I've declared a 2D array of integers. For Dirichlet, it's a N-length array of 4-vectors.

Second, I've gotten rid of the transformed data block. Here, you want to declare the B's in the original model as a simplex. There's no need to try to do any trickiness -- Stan can handle it all under the hood for you.

If you wanted to try it out on more simple data first, you can replace the data block with a transformed data block like:

transformed data {
  int <lower=0> N;
  vector[4] ns[100];
  N <- 100;
  for (n in 1:100) {
    ns[n,1] <- 1;
    ns[n,2] <- 10;
    ns[n,3] <- 100;
    ns[n,4] <- 1000;
  }
}

Hopefully that helps and you can really simplify your model. If you want documentation on the distributions, you'll find them in the manual.

When we suggest breaking down your model, these are sort of the steps we take. Even with the Dirichlet example, I put in fake data (as shown), where I knew the expected parameter values. I had original swapped the parameters and the variables since I was working with the multinomial first (written as 'ns[n] ~ dirichlet(B)') and it was clearly not right. I wouldn't have been able to pick that out if I used the full dataset to begin with. I also wouldn't have been able to debug it as quickly if there were lots of things happening in transformation blocks or restrictive priors. I guess I don't have the ability to intuit parameters through a thick model.

Anyway, let us know how you get on with the rest of your model.

Daniel

[Nick Menzies]

Hi Daniel,

Thank you for laying this out so clearly. Sorry the simple model I sent did not work --> I forgot that you could specify simplex vectors, so instead that model had an awkward (and fatal, it seems) way of achieving the same result.

The suggestion of a Dirichlet is interesting. I raised a question about it earlier on this thread, and went with the current approach as there seemed to be some trickiness involved with the Dirichlet, in that (if I am interpreting the comments and the manual correctly), the Dirichlet should really be taking the counts plus 1 as its argument (i.e. if one had integer data,  the likelihood of counts ~ Multinomial(theta) would be proportional to theta ~ Dirichlet(counts+1)). I guess this extra step really is a trivial problem, but please confirm that this step (add one to the counts) is needed.   

Also, do you know if Dirichlet will be vectorised in Stan 2.0?  Presumably this would make Dirichlet the faster option as well.

For the transformed parameters section: I agree it is... busy. I need the detail to answer my research question, but there are many opportunities for bugs to come in. However, I have coded up a version in R and rcpp so will be able to use this to generate better inits (I think this has been a problem). This will also let me compare model results across R and Stan to check for bugs (I normally work in R so am more confident about spotting bugs there). 

Also -- great idea about fake data, will give it a go.  

Nick

[Bob Carpenter]

On 10/1/13 8:56 PM, Nick Menzies wrote:
> Hi Daniel,
>
> Thank you for laying this out so clearly. Sorry the simple model I sent did not work --> I forgot that you could specify
> simplex vectors, so instead that model had an awkward (and fatal, it seems) way of achieving the same result.
>
> The suggestion of a Dirichlet is interesting. I raised a question about it earlier on this

> <https://groups.google.com/forum/#!searchin/stan-users/bob$20multinomial$20dirichlet/stan-users/55-Dr6JW-yA/FQ2uwMtKd0cJ>thread,

> and went with the current approach as there seemed to be some trickiness involved with the Dirichlet, in that (if I am

> interpreting the comments and the manual correctly), the Dirichlet should really be taking the counts */plus 1/* as its

> argument (i.e. if one had integer data, the likelihood of counts ~ Multinomial(theta) would be proportional to theta ~
> Dirichlet(counts+1)). I guess this extra step really is a trivial problem, but please confirm that this step (add one to
> the counts) is needed.

The Stan parameterization of the Dirichlet is the usual one:

http://en.wikipedia.org/wiki/Dirichlet_distribution

You can see from the density that if the parameter vector is
all 1 values, it's a uniform prior. So if you've seen 4 outcomes
of category 1, 3 of category 2 and 1 of category 3, and you
use a Dirichlet({1,1,1}) prior, then the posterior is
Dirichlet({5,4,2}), which is why I said the parameter is 1 plus
the prior number of observations in the doc.

> Also, do you know if Dirichlet will be vectorised in Stan 2.0? Presumably this would make Dirichlet the faster option
> as well.

Not yet. We're still working through the vectorization of
multivariates. The metaprogramming's tricky, so we've been putting
it off. But it's a high priority for us, so it'll be done soon
(by that, I mean probably over the next six months).

> For the transformed parameters section: I agree it is... busy. I need the detail to answer my research question, but
> there are many opportunities for bugs to come in. However, I have coded up a version in R and rcpp so will be able to
> use this to generate better inits (I think this has been a problem).

Reasonable inits help with convergence.

- Bob

[Andrew Gelman]

Regarding the noninformative parameterization of the beta (and, by extension, the Dirichlet):
http://andrewgelman.com/2011/12/07/neutral-noninformative-and-informative-conjugate-beta-and-gamma-prior-distributions/
Jouni recommends Beta (1/3,1/3) as default. I'm not sure what he'd do with Dirichlet.
A