Best practices...
68 views
Skip to first unread message
Aki Vehtari
Apr 7, 2015, 10:48:50 AM4/7/15
to stan...@googlegroups.com
Related to the discussion of best practices for translating BUGS models, I'll start asking questions
https://github.com/stan-dev/example-models/blob/master/bugs_examples/vol1/stacks/stacks_a_normal.stan
has
---clip---
// to standardize the x's
transformed data {
matrix[N,p] z;
row_vector[p] mean_x;
real sd_x[p];
for (j in 1:p) {
mean_x[j] <- mean(col(x,j));
sd_x[j] <- sd(col(x,j));
for (i in 1:N)
z[i,j] <- (x[i,j] - mean_x[j]) / sd_x[j];
}
}
---clip---
I think mean_x and sd_x should be same type. Which is better, and why row_vector and not vector? It's not affecting the performance, but it's just confusing now.
Aki
https://github.com/stan-dev/example-models/blob/master/bugs_examples/vol1/stacks/stacks_a_normal.stan
has
---clip---
// to standardize the x's
transformed data {
matrix[N,p] z;
row_vector[p] mean_x;
real sd_x[p];
for (j in 1:p) {
mean_x[j] <- mean(col(x,j));
sd_x[j] <- sd(col(x,j));
for (i in 1:N)
z[i,j] <- (x[i,j] - mean_x[j]) / sd_x[j];
}
}
---clip---
I think mean_x and sd_x should be same type. Which is better, and why row_vector and not vector? It's not affecting the performance, but it's just confusing now.
Aki
Bob Carpenter
Apr 7, 2015, 11:57:44 AM4/7/15
to stan...@googlegroups.com
OK, here's why I think this isn't going to be so easy.
First, keep in mind that what's best practice for efficiency
isn't the same as what's best practice for clarity. Which
are we aiming for? How direct should the translations of the
BUGS models be?
I was originally suggesting
* a line-by-line very direct translation
* a conversion of that to one using the priors, etc. we like,
possibly rearranging data structures for efficiency, etc.
* an efficient version of the good model
* optionally an efficient version of the direct translation
The best practices on arrays vs. matrices are in the
manual chapter on arrays vs. matrices, so I suppose the
duplication of effort can begin by shoveling that chapter
into the best practices Wiki.
Then there's the advice to always use matrix functions
and shorter functions like elementwise functions if
you can to shorten the code.
But then given that this is just the transformed data
block, the vectorization and so forth aren't that important
since they only get executed once. There's a section in the
manual on execution that's important to read to understand when
to put things in transformed data, transformed parameters, and
generated quantities.
Best practice is to always declare things as transformed data or
generated quantities if possible, always delcare things as local
variables if possible. Lots of users inadvertently make everything
a transformed parameter or local variable in the model that should
be transformed data or a generated quantity. So this all needs
to go in the best practices.
Another best practice is to use capital letters for
constants and then matching lower-case letters for loops.
This won't match BUGS, so I'm not sure where that goes
in the best practices vs. translation of BUGS models.
And follow Gelman et al. in variable naming, so no capital
letters like Y.
Are mean_x and sd_x used elsewhere? If not,
best practicses would dictate they be local variables.
This is so that their memory doesn't linger past when
they're used. I see they are, but this should be mentioned
as a best practice if we're making a list.
We also need best practices for formatting. I'm tired of
reading crazily formatted user models --- there's an appendix
to the manual for that.
Now onto the show.
On Apr 8, 2015, at 12:48 AM, Aki Vehtari <aki.v...@aalto.fi> wrote:
>
> Related to the discussion of best practices for translating BUGS models, I'll start asking questions
>
> https://github.com/stan-dev/example-models/blob/master/bugs_examples/vol1/stacks/stacks_a_normal.stan
> has
> ---clip---
> // to standardize the x's
> transformed data {
> matrix[N,p] z;
> row_vector[p] mean_x;
> real sd_x[p];
> for (j in 1:p) {
> mean_x[j] <- mean(col(x,j));
> sd_x[j] <- sd(col(x,j));
> for (i in 1:N)
> z[i,j] <- (x[i,j] - mean_x[j]) / sd_x[j];
> }
> }
So I'd suggest this:
matrix[N,M] x;
transformed data {
matrix[N,M] z;
vector[M] mean_x;
vector[M] sd_x;
for (m in 1:M) {
mean_x[m] <- mean(col(x,m));
sd_x[m] <- mean(col(x,m));
}
for (n in 1:N)
z[n] <- (x[n] - mean_x) ./ sd_x;
}
But maybe that ./ is too complicated to read?
Or maybe people don't like z[n] referring to a row
of a matrix?
This isn't ideal for speed, since the col(x,m)
operations have to copy. So you can make a local
col_x_m <- col(x,m)
to reuse it. A critical piece of best practices is
to define expressions only once and reuse them --- it's very
important for efficiency, but can make things harder to read.
And then you want to rearrange and the z[n] assignments are
not memory local, so aren't efficient for large matrices. But
then there's a tradeoff with wanting to be able to define mu
with a matrix operation. Also, do we really need to define
mu and sigma as transformed parameters? It's overhead in writing
and in print and whatnot. If not, use this:
y ~ normal(beta0 + z * beta, sqrt(sigma_sq));
I like the "_" to separate words in names, but your mileage may
vary.
So many pieces of advice!
For the rest of that model, we don't ever really want to
be using step functions --- relational operators are much clearer, so
I'd rewrite
outlier_3 <- step(fabs((Y[3] - mu[3]) / sigma) - 2.5);
as
outlier_3 <- fabs(Y[3] - mu[3]) / sigma - 2.5 > 0
or better still
outlier_3 <- fabs(Y[3] - mu[3]) / sigma > 2.5
or even better for efficiency (multiplication is faster than
division),
outlier_3 <- fabs(Y[3] - mu[3]) > 2.5 * sigma;
Best practice is not to include parens around things where the
language knows how to parameterize. One might be tempted to
write ... > (2.5 * sigma).
Of course, rather than writing the same thing three times
with cut-and-paste, best practice would dictate writing
a function and actually declaring the booleans properly as
constrained integers. First the function.
int is_outlier(double y, double mu, double sigma) {
return fabs(y - mu) > 2.5 * sigma;
}
Then declare outlier_3, etc. as int<lower=0,upper=1>.
Best practice is to always have bounds for error checking and
self-documenting code.
Also for Andrew, random variables are not capitalized, so our
best practice notation is to use lower-case y.
Of course, Andrew doesn't like those priors, so then there's
the step of actually making the priors better and stating
best practices for those. There's a bunch of advice extracted
from Andrew in the regression chapter of the manual. But it's
quite involved and varies based on whether you're doing MCMC or
optimization.
Also, this
for (j in 1:p)
b[j] <- beta[j] / sd_x[j];
should be rewritten as:
b <- beta ./ sd_x;
But then you see why you want row vectors vs. vectors
so that you can do all the matrix algebra with minimal
transpositions. With mean_x being a row vector and b being
a vector, (mean_x * b) is a dot product. Or you can use dot_product(a,b),
which is more efficient than making mean_x a vector and using
mean_x' * b (tranpose makes a copy).
Is BUGS (best practice spelling!) really taking
beta ~ normal(0, 1e5)?
If so, we want full precision, so the exact translation
should be:
beta ~ normal(0, sqrt(1e5));
or to eliminate the sqrt operation without losing precision,
either define a transformed data variable or use this
beta ~ normal(0, 316.227766016838);
I prefer 1e5 to 100000, so I think scientific notation should
be a best practice if we're drilling down to that level.
- Bob
> ---clip---
>
> I think mean_x and sd_x should be same type. Which is better, and why row_vector and not vector? It's not affecting the performance, but it's just confusing now.
>
> Aki
>
> --
> You received this message because you are subscribed to the Google Groups "stan development mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+u...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
First, keep in mind that what's best practice for efficiency
isn't the same as what's best practice for clarity. Which
are we aiming for? How direct should the translations of the
BUGS models be?
I was originally suggesting
* a line-by-line very direct translation
* a conversion of that to one using the priors, etc. we like,
possibly rearranging data structures for efficiency, etc.
* an efficient version of the good model
* optionally an efficient version of the direct translation
The best practices on arrays vs. matrices are in the
manual chapter on arrays vs. matrices, so I suppose the
duplication of effort can begin by shoveling that chapter
into the best practices Wiki.
Then there's the advice to always use matrix functions
and shorter functions like elementwise functions if
you can to shorten the code.
But then given that this is just the transformed data
block, the vectorization and so forth aren't that important
since they only get executed once. There's a section in the
manual on execution that's important to read to understand when
to put things in transformed data, transformed parameters, and
generated quantities.
Best practice is to always declare things as transformed data or
generated quantities if possible, always delcare things as local
variables if possible. Lots of users inadvertently make everything
a transformed parameter or local variable in the model that should
be transformed data or a generated quantity. So this all needs
to go in the best practices.
Another best practice is to use capital letters for
constants and then matching lower-case letters for loops.
This won't match BUGS, so I'm not sure where that goes
in the best practices vs. translation of BUGS models.
And follow Gelman et al. in variable naming, so no capital
letters like Y.
Are mean_x and sd_x used elsewhere? If not,
best practicses would dictate they be local variables.
This is so that their memory doesn't linger past when
they're used. I see they are, but this should be mentioned
as a best practice if we're making a list.
We also need best practices for formatting. I'm tired of
reading crazily formatted user models --- there's an appendix
to the manual for that.
Now onto the show.
On Apr 8, 2015, at 12:48 AM, Aki Vehtari <aki.v...@aalto.fi> wrote:
>
> Related to the discussion of best practices for translating BUGS models, I'll start asking questions
>
> https://github.com/stan-dev/example-models/blob/master/bugs_examples/vol1/stacks/stacks_a_normal.stan
> has
> ---clip---
> // to standardize the x's
> transformed data {
> matrix[N,p] z;
> row_vector[p] mean_x;
> real sd_x[p];
> for (j in 1:p) {
> mean_x[j] <- mean(col(x,j));
> sd_x[j] <- sd(col(x,j));
> for (i in 1:N)
> z[i,j] <- (x[i,j] - mean_x[j]) / sd_x[j];
> }
> }
matrix[N,M] x;
transformed data {
matrix[N,M] z;
vector[M] mean_x;
vector[M] sd_x;
for (m in 1:M) {
mean_x[m] <- mean(col(x,m));
sd_x[m] <- mean(col(x,m));
}
for (n in 1:N)
z[n] <- (x[n] - mean_x) ./ sd_x;
}
But maybe that ./ is too complicated to read?
Or maybe people don't like z[n] referring to a row
of a matrix?
This isn't ideal for speed, since the col(x,m)
operations have to copy. So you can make a local
col_x_m <- col(x,m)
to reuse it. A critical piece of best practices is
to define expressions only once and reuse them --- it's very
important for efficiency, but can make things harder to read.
And then you want to rearrange and the z[n] assignments are
not memory local, so aren't efficient for large matrices. But
then there's a tradeoff with wanting to be able to define mu
with a matrix operation. Also, do we really need to define
mu and sigma as transformed parameters? It's overhead in writing
and in print and whatnot. If not, use this:
y ~ normal(beta0 + z * beta, sqrt(sigma_sq));
I like the "_" to separate words in names, but your mileage may
vary.
So many pieces of advice!
For the rest of that model, we don't ever really want to
be using step functions --- relational operators are much clearer, so
I'd rewrite
outlier_3 <- step(fabs((Y[3] - mu[3]) / sigma) - 2.5);
as
outlier_3 <- fabs(Y[3] - mu[3]) / sigma - 2.5 > 0
or better still
outlier_3 <- fabs(Y[3] - mu[3]) / sigma > 2.5
or even better for efficiency (multiplication is faster than
division),
outlier_3 <- fabs(Y[3] - mu[3]) > 2.5 * sigma;
Best practice is not to include parens around things where the
language knows how to parameterize. One might be tempted to
write ... > (2.5 * sigma).
Of course, rather than writing the same thing three times
with cut-and-paste, best practice would dictate writing
a function and actually declaring the booleans properly as
constrained integers. First the function.
int is_outlier(double y, double mu, double sigma) {
return fabs(y - mu) > 2.5 * sigma;
}
Then declare outlier_3, etc. as int<lower=0,upper=1>.
Best practice is to always have bounds for error checking and
self-documenting code.
Also for Andrew, random variables are not capitalized, so our
best practice notation is to use lower-case y.
Of course, Andrew doesn't like those priors, so then there's
the step of actually making the priors better and stating
best practices for those. There's a bunch of advice extracted
from Andrew in the regression chapter of the manual. But it's
quite involved and varies based on whether you're doing MCMC or
optimization.
Also, this
for (j in 1:p)
should be rewritten as:
b <- beta ./ sd_x;
But then you see why you want row vectors vs. vectors
so that you can do all the matrix algebra with minimal
transpositions. With mean_x being a row vector and b being
a vector, (mean_x * b) is a dot product. Or you can use dot_product(a,b),
which is more efficient than making mean_x a vector and using
mean_x' * b (tranpose makes a copy).
Is BUGS (best practice spelling!) really taking
beta ~ normal(0, 1e5)?
If so, we want full precision, so the exact translation
should be:
beta ~ normal(0, sqrt(1e5));
or to eliminate the sqrt operation without losing precision,
either define a transformed data variable or use this
beta ~ normal(0, 316.227766016838);
I prefer 1e5 to 100000, so I think scientific notation should
be a best practice if we're drilling down to that level.
- Bob
> ---clip---
>
> I think mean_x and sd_x should be same type. Which is better, and why row_vector and not vector? It's not affecting the performance, but it's just confusing now.
>
> Aki
>
> You received this message because you are subscribed to the Google Groups "stan development mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+u...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
Andrew Gelman
Apr 7, 2015, 6:28:09 PM4/7/15
to stan...@googlegroups.com
Hi, there are 2 issues:
1. Translation of Bugs/Jags models to Stan. I don’t think we need a line by line translation. I think what is best is, for each model, for us to give the Stan equivalent, but changing anything that’s bad practice and explaining the changes in comments (for example, we don’t do inv-gamma(.001,.001) because of Gelman, 2006, and explaining how we set up the mixture model following the recommendations of the Stan manual, and explaining why we use normal(0,10) priors rather than normal with sd 1000. A line by line translation would be fine but I don’t think it’s necessary, and it does represent effort, it represents effort that could better be used elsewhere, I think.
Regarding this translation, I think we should use what you (Bob) consider best practice in every case. Including naming of variables (for example, n=1,…,N rather than i=1,…,n). I mean, we don’t have to kill ourselves. If the model’s already done we can leave it as is. But if we have the option, if we’re choosing what to do, I I’d say let’s do what Bob would do.
2. A wiki for best practices in Bayes and Stan. For this, ambiguity is ok. It’s fine that what’s best practice for efficiency isn’t the same as what’s best practice for clarity (and that clarity is in the eye of the beholder). We can explain these things in the wiki.
A
1. Translation of Bugs/Jags models to Stan. I don’t think we need a line by line translation. I think what is best is, for each model, for us to give the Stan equivalent, but changing anything that’s bad practice and explaining the changes in comments (for example, we don’t do inv-gamma(.001,.001) because of Gelman, 2006, and explaining how we set up the mixture model following the recommendations of the Stan manual, and explaining why we use normal(0,10) priors rather than normal with sd 1000. A line by line translation would be fine but I don’t think it’s necessary, and it does represent effort, it represents effort that could better be used elsewhere, I think.
Regarding this translation, I think we should use what you (Bob) consider best practice in every case. Including naming of variables (for example, n=1,…,N rather than i=1,…,n). I mean, we don’t have to kill ourselves. If the model’s already done we can leave it as is. But if we have the option, if we’re choosing what to do, I I’d say let’s do what Bob would do.
2. A wiki for best practices in Bayes and Stan. For this, ambiguity is ok. It’s fine that what’s best practice for efficiency isn’t the same as what’s best practice for clarity (and that clarity is in the eye of the beholder). We can explain these things in the wiki.
A
Bob Carpenter
Apr 7, 2015, 7:29:43 PM4/7/15
to stan...@googlegroups.com
The point was that Sebastian (and others) were asking for
us to translate
*** the exact same model ***
into Stan so that
*** they could compare the results are the same. ***
That would specifically entail
*** not changing the priors. ***
A secondary goal was just helping people learn how to code
in Stan by translating the BUGS models they already know and
love in the simplest possible way.
We know
* Michael has a bunch of issues with any such comparisons
from a statistical computation point of view (e.g., tail
probabilities and ergodicity),
* Andrew has a bunch of issues with doing this from a statistical
methodology and modeling point of view (e.g., priors), and
* I have a bunch of issues with doing this from a Stan programming
point of view (e.g., local variables).
- Bob
us to translate
*** the exact same model ***
into Stan so that
*** they could compare the results are the same. ***
That would specifically entail
*** not changing the priors. ***
A secondary goal was just helping people learn how to code
in Stan by translating the BUGS models they already know and
love in the simplest possible way.
We know
* Michael has a bunch of issues with any such comparisons
from a statistical computation point of view (e.g., tail
probabilities and ergodicity),
* Andrew has a bunch of issues with doing this from a statistical
methodology and modeling point of view (e.g., priors), and
* I have a bunch of issues with doing this from a Stan programming
point of view (e.g., local variables).
- Bob
Andrew Gelman
Apr 7, 2015, 7:55:22 PM4/7/15
to stan...@googlegroups.com
I can see that Sebastian and others might want to do this. But I don’t think such non-recommended models should be included in the official set of Stan models.
Bob Carpenter
Apr 7, 2015, 8:05:42 PM4/7/15
to stan...@googlegroups.com
> On Apr 8, 2015, at 9:55 AM, Andrew Gelman <gel...@stat.columbia.edu> wrote:
>
> I can see that Sebastian and others might want to do this. But I don’t think such non-recommended models should be included in the official set of Stan models.
to be the "official set of Stan models"?
I thought you just recently said you were OK having direct translations
of BUGS models, gamma(0.001,0.001), normal(0,316....), and all,
as long as we put big warning labels on them with directions on where
to find the recommended alternatives. If we do decide to distribute
direct translations ("official" or not), you and Sebastian will have
to negotiate how strong to make those warnings --- he made the compelling
point that we didn't want to offend BUGS users, but gradually bring lure
into the fold.
I see the advantages and disadvantages of both including them and
not including them and honestly don't care one way or the other. So I'm
going to bow out and let other people try to arrive at a final policy.
- Bob
Andrew Gelman
Apr 7, 2015, 8:11:32 PM4/7/15
to stan...@googlegroups.com
> On Apr 7, 2015, at 8:04 PM, Bob Carpenter <ca...@alias-i.com> wrote:
>
>> On Apr 8, 2015, at 9:55 AM, Andrew Gelman <gel...@stat.columbia.edu> wrote:
>>
>> I can see that Sebastian and others might want to do this. But I don’t think such non-recommended models should be included in the official set of Stan models.
>
> I must be misunderstanding some subtlety here. What do you take
> to be the "official set of Stan models"?
>
> I thought you just recently said you were OK having direct translations
> of BUGS models, gamma(0.001,0.001), normal(0,316....), and all,
> as long as we put big warning labels on them with directions on where
> to find the recommended alternatives.
Michael Betancourt
Apr 8, 2015, 2:54:14 AM4/8/15
to stan...@googlegroups.com
Andrew,
I agree that we don’t want to propagate bad models, but if we
were really serious about that then we wouldn’t be showing many
BUGS models at all! The problem is that people want to use
these translations to _learn_ the Stan modeling language and
so they need direct translations in order to transfer their knowledge
from BUGS into Stan. The _only_ way we’ll be able to have
any utility from this perspective is with two stages — first the direct
translation and only then the “proper” translation with modeling and
code improvements.
I agree that we don’t want to propagate bad models, but if we
were really serious about that then we wouldn’t be showing many
BUGS models at all! The problem is that people want to use
these translations to _learn_ the Stan modeling language and
so they need direct translations in order to transfer their knowledge
from BUGS into Stan. The _only_ way we’ll be able to have
any utility from this perspective is with two stages — first the direct
translation and only then the “proper” translation with modeling and
code improvements.
Aki Vehtari
Apr 8, 2015, 3:15:48 AM4/8/15
to stan...@googlegroups.com
I think the current BUGS example models are quite direct translations, and we don't need to repeat that. Just let them be there and add warning that they don't follow the best practices and meant only for people familiar with BUGS/JAGS.
I'm teaching Stan to people who have never used BUGS or gamma priors, and I don't want to spend time explaining the wrong models. That's why I'm willing to spend my time to add BUGS derived "data analysis examples" using the best practices in Stan. I wrote "data analysis examples", implying that we can make also other changes than just prior changes. For example, I could add log_lik computation for LOO for all models where it's applicable. Data analysis examples in BUGS example volumes are nice because they are small and have already the back story written, but also include many types of hierarchical models.
Aki
I'm teaching Stan to people who have never used BUGS or gamma priors, and I don't want to spend time explaining the wrong models. That's why I'm willing to spend my time to add BUGS derived "data analysis examples" using the best practices in Stan. I wrote "data analysis examples", implying that we can make also other changes than just prior changes. For example, I could add log_lik computation for LOO for all models where it's applicable. Data analysis examples in BUGS example volumes are nice because they are small and have already the back story written, but also include many types of hierarchical models.
Aki
Aki Vehtari
Apr 8, 2015, 6:28:01 AM4/8/15
to stan...@googlegroups.com
Thanks Bob for the detailed answer. I think quite many of the BUGS models could be translated just following this example. I know that most of these suggestions can be found in the manual, but they are scattered around.
> for (m in 1:M) {
> mean_x[m] <- mean(col(x,m));
> sd_x[m] <- mean(col(x,m));
> }
> for (n in 1:N)
> z[n] <- (x[n] - mean_x) ./ sd_x;
> }
>
> But maybe that ./ is too complicated to read?
> Or maybe people don't like z[n] referring to a row
> of a matrix?
I think I'll wait for the slicing operator and then write
for (m in 1:M)
z[:,m] <- (x[:,m] - mean(x[:,m])) / sd(x[:,m]);
> A critical piece of best practices is
> to define expressions only once and reuse them --- it's very
> important for efficiency, but can make things harder to read.
I think it's better to provide easy to read examples, and optimise for speed only when necessary. Naturally you could show some examples optimised for speed, but I don't think any of the BUGS examples are so slow that they would be good examples for that.
Aki
> for (m in 1:M) {
> mean_x[m] <- mean(col(x,m));
> sd_x[m] <- mean(col(x,m));
> }
> for (n in 1:N)
> z[n] <- (x[n] - mean_x) ./ sd_x;
> }
>
> But maybe that ./ is too complicated to read?
> Or maybe people don't like z[n] referring to a row
> of a matrix?
for (m in 1:M)
> A critical piece of best practices is
> to define expressions only once and reuse them --- it's very
> important for efficiency, but can make things harder to read.
Aki
Andrew Gelman
Apr 9, 2015, 12:46:28 AM4/9/15
to stan...@googlegroups.com
M
I guess we’d have to talk with some users. I actually think that a translation that keeps or improves functionality is fine, as long as we explain.
A
I guess we’d have to talk with some users. I actually think that a translation that keeps or improves functionality is fine, as long as we explain.
A
Reply all
Reply to author
Forward
0 new messages