exposing stan functions to R requires at least two parameters???

385 views
Skip to first unread message

WallSt

unread,
Oct 19, 2015, 4:32:14 PM10/19/15
to Stan users mailing list
I am getting an odd error when trying to expose a Stan function to R.  The function works if I have two input parameters in the Stan function, the second being a dummy, but returns an error if only one parameter is required.  The following R code with "mc" works, but "mc2" returns an error "In Rcpp::sourceCpp(code = paste(lines, collapse = "\n")) :
  No Rcpp::export attributes or RCPP_MODULE declarations found in source"

The only difference appears to be that the first function has two inputs, the second has one. (There also seems to be a bug when you write a function with no input parameters)

mc <- '
functions {
    vector fib(int n, real x) {
      vector[n] fib;
      fib[1] <- 1;
      fib[2] <- 1;
      for (i in 3:n) {
        fib[i] <- fib[i-1] + fib[i-2];
      }
      return fib;
    }
}

model {

}  
'

mc2 <- 
  '
functions {
    vector fib2(int n) {
      vector[n] fib2;
      fib2[1] <- 1;
      fib2[2] <- 1;
      for (i in 3:n) {
        fib2[i] <- fib2[i-1] + fib2[i-2];
      }
      return fib2;
    }
}

model {

}  
'
cppcode <- stanc(model_code = mc, model_name = "Demonstration2")
expose_stan_functions(cppcode)
fib(200, x=1)

cppcode2 <- stanc(model_code = mc2, model_name = "Demonstration3")
expose_stan_functions(cppcode2)
fib2(200)

Ben Goodrich

unread,
Oct 19, 2015, 4:48:28 PM10/19/15
to Stan users mailing list
On Monday, October 19, 2015 at 4:32:14 PM UTC-4, WallSt wrote:
I am getting an odd error when trying to expose a Stan function to R.  The function works if I have two input parameters in the Stan function, the second being a dummy, but returns an error if only one parameter is required.  The following R code with "mc" works, but "mc2" returns an error "In Rcpp::sourceCpp(code = paste(lines, collapse = "\n")) :
  No Rcpp::export attributes or RCPP_MODULE declarations found in source"

I believe this was fixed after rstan 2.8.0 was released in response to bug #214. Try executing


before calling expose_stan_functions().

Ben

WallSt

unread,
Oct 19, 2015, 5:23:15 PM10/19/15
to Stan users mailing list
Thanks Ben.  This fixes the one vs. two parameter problem.

However, please note  that if the function has no parameters, then the model does not compile.  This is not fixed by your workaround link.  For example, try the following fibonacci algorithm which is fixed at 20 steps.  "fib3" is not returned. (It was not clear to me from the Stan manual whether functions could be declared without any input parameters.)

mc3 <- 
  '
functions {
    vector fib3() {
      vector[20] fib3;
      fib3[1] <- 1;
      fib3[2] <- 1;
      for (i in 3:20) {
        fib3[i] <- fib3[i-1] + fib3[i-2];
      }
      return fib3;
    }
}

model {}  
'
cppcode3 <- stanc(model_code = mc3, model_name = "Demonstration3")
expose_stan_functions(cppcode3)
fib3()

Ben Goodrich

unread,
Oct 19, 2015, 5:45:55 PM10/19/15
to Stan users mailing list
On Monday, October 19, 2015 at 5:23:15 PM UTC-4, WallSt wrote:
However, please note  that if the function has no parameters, then the model does not compile.  This is not fixed by your workaround link.  For example, try the following fibonacci algorithm which is fixed at 20 steps.  "fib3" is not returned. (It was not clear to me from the Stan manual whether functions could be declared without any input parameters.)

It is legal but it hadn't occurred to me that anyone would do so. I'll try to fix it before the next rstan release.

Thanks,
Ben
 

Jan Tünnermann

unread,
Oct 19, 2015, 7:09:33 PM10/19/15
to stan-...@googlegroups.com
Dear all,

I'm still having trouble with non-efficient sampling in my hierarchical
model (I already tried some things, see earlier e-mails). Intriguingly,
often all four chains behave nicely except one. Or suddenly, out of
nowhere, a chain get trapped for a long period.

I prepared a simplified* version of this model and attached it. Maybe
someone can have a quick look if there are some obvious mistakes. We
haven't been working with Bayesian statistics for too long and
especially I switched only recently to Stan (in the hope HMC can better
deal with this model).

The general structure of the model is as follows (in bottom-up order):
- on the stimulus level layer: positive outcomes are binomially
distributed at each level
- on the subject layer: a custom psychometric function describes the
success rates at different levels
- on the group layer: parameters of the mentioned function are drawn
from normal group distributions

Many thanks and best regards
Jan

* I made the following simplifications:
- deleted model structure of a second identical condition
- deleted a huge generated quantities block, mostly containing
differences between the two conditions
- included dummy data within the script. Only 10 data sets are
included, whereas the model is intended to run on 20+ data sets
eventually. With half the data it is even more stubborn ...

P.S.: I jut obtained a 3000 iterations test sample with the file that
doesn't look too bad. However, I've seen this before and faced the
mention problems then after 20000 iterations or so ...
evaluate-toj-full-simplified.R

Andrew Gelman

unread,
Oct 19, 2015, 7:17:20 PM10/19/15
to stan-...@googlegroups.com
Hi, here are some quick thoughts:
1. We recommend putting the Stan program in a separate .stan file rather than encoding it as a character string. That’ll be easier when you want to fit the model more than once.

2. Hard constraints on parameters are typically a good idea. Inherent constraints (such as sigma>0) are fine, but constraints such as <-.2,.2> are typically a badd idea. Better to use a prior, for example replacing a hard <-.2,.2> constraint by a prior such as N(0,.1) or N(0,.2).
3. Similarly for the parameter you bound between 0.01 and 0.99. If it can’t be 0.995, it probably is highly unlikely to be 0.985. Etc.
4. No need to write 0 and 8 as 0.0 and 8.0.
5. The Cauchy priors could be a bit weak. I suggest replacing by normal or something like t_3 or t_7 if you want some long tails.
I suspect that none of these tips will solve your problem, but who knows? Maybe! I’m thinking it can’t hurt to put in some priors.
A
> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
> <evaluate-toj-full-simplified.R>


Jan Tünnermann

unread,
Oct 19, 2015, 7:41:07 PM10/19/15
to stan-...@googlegroups.com
Dear Andrew,

thanks for the quick reply.
On 10/20/2015 01:17 AM, Andrew Gelman wrote:
> Hi, here are some quick thoughts:
> 1. We recommend putting the Stan program in a separate .stan file rather than encoding it as a character string. That’ll be easier when you want to fit the model more than once.
Great, I was already wondering how to make use of that feature (because
I read people talking about fitting models without recompiling).
> 2. Hard constraints on parameters are typically a good idea. Inherent constraints (such as sigma>0) are fine, but constraints such as <-.2,.2> are typically a badd idea. Better to use a prior, for example replacing a hard <-.2,.2> constraint by a prior such as N(0,.1) or N(0,.2).
OK, I'll try that. I must admit I was using the uniform ranges becaue I
can imagine the range from which samples are drawn more intuitively.
> 3. Similarly for the parameter you bound between 0.01 and 0.99. If it can’t be 0.995, it probably is highly unlikely to be 0.985. Etc.
That's true. The idea was to not put in any information (which I
probalby still do) but avoid extreme values that may lead to problems in
the functions ...
These parameters are acutally biases between 0 and 1 that are neutral at
0.5. Is there any prior distribution favorable for such cases?
> 4. No need to write 0 and 8 as 0.0 and 8.0.
I always do so if a language allows it. Can avoid hard-to-find erros is
statically typed languages. Also helps me to keep track whether there is
a continous range or discrete integer values.
> 5. The Cauchy priors could be a bit weak. I suggest replacing by normal or something like t_3 or t_7 if you want some long tails.
I think you may be right here. I have been following the advice from
your (2003?) paper on priors for variance parameters.Hence the Cauchys.
But I observed in the posteriors that the sigmas of many variables look
rather normally distributed.

The problem is, I have little feeling for how informative certain priors
(or contraints) may be. For simpler models in JAGS I just had it
generate prior predictive distirbutions by running the model without
data. I'm not sure if that is an option here as well.

> I suspect that none of these tips will solve your problem, but who knows? Maybe! I’m thinking it can’t hurt to put in some priors.

> A
Thanks again,
Jan
>> On Oct 19, 2015, at 7:09 PM, Jan Tünnermann <je...@mail.upb.de> wrote:
>>
>> Dear all,
>>
>> I'm still having trouble with non-efficient sampling in my hierarchical model (I already tried some things, see earlier e-mails). Intriguingly, often all four chains behave nicely except one. Or suddenly, out of nowhere, a chain get trapped for a long period.
>>
>> I prepared a simplified* version of this model and attached it. Maybe someone can have a quick look if there are some obvious mistakes. We haven't been working with Bayesian statistics for too long and especially I switched only recently to Stan (in the hope HMC can better deal with this model).
>>
>> The general structure of the model is as follows (in bottom-up order):
>> - on the stimulus level layer: positive outcomes are binomially distributed at each level
>> - on the subject layer: a custom psychometric function describes the success rates at different levels
>> - on the group layer: parameters of the mentioned function are drawn from normal group distributions
>>
>> Many thanks and best regards
>> Jan
>>
>> * I made the following simplifications:
>> - deleted model structure of a second identical condition
>> - deleted a huge generated quantities block, mostly containing differences between the two conditions
>> - included dummy data within the script. Only 10 data sets are included, whereas the model is intended to run on 20+ data sets eventually. With half the data it is even more stubborn ...
>>
>> P.S.: I jut obtained a 3000 iterations test sample with the file that doesn't look too bad. However, I've seen this before and faced the mention problems then after 20000 iterations or so ...
>>
>> --
>> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
>> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
>> To post to this group, send email to stan-...@googlegroups.com.
>> For more options, visit https://groups.google.com/d/optout.
>> <evaluate-toj-full-simplified.R>
>


--
M.Sc. Jan T U E N N E R M A N N _________________________________

University of Paderborn
Faculty of Computer Science
Electrical Engineering and Mathematics
GET Lab

phone: +49/ (0)5251/ 60 - 2219 Pohlweg 47-49
fax: +49/ (0)5251/ 60 - 3238 D-33098 Paderborn
email: tuenn...@get.uni-paderborn.de

Faculty of Cultural Sciences
Department of Psychology
PsyLab
phone: +49/ (0)5251/ 60 - 3211 Warburger Straße 100
fax: +49/ (0)5251/ 60 - 3528 D-33098 Paderborn
email: jan.tue...@uni-paderborn.de

Bob Carpenter

unread,
Oct 20, 2015, 10:12:36 AM10/20/15
to stan-...@googlegroups.com
Yes, you can define constants or PRNGs as functions with no arguments.

- Bob

Bob Carpenter

unread,
Oct 20, 2015, 10:17:41 AM10/20/15
to stan-...@googlegroups.com
The cauchy's just a weakly informative prior --- if you have more info
about your distribution's tails, by all means use it!

The only place integers should matter in Stan is (1.0 / 3.0) vs. (1 / 3).

Chains typically get stuck in hierarchical models with wildly varying posterior
scales. You can try reducing initial step size and increasing the
target acceptance rate, which will slow down the model but make it more
robust to getting stuck.

- Bob

Jan Tünnermann

unread,
Oct 20, 2015, 10:28:07 AM10/20/15
to stan-...@googlegroups.com
Hi Bob,

On 10/20/2015 04:17 PM, Bob Carpenter wrote:
> The cauchy's just a weakly informative prior --- if you have more info
> about your distribution's tails, by all means use it!
Do observations of serious trial and error runs count as information in
this case (say, I try a lot of different priors, but the posterior
always looks like a normal with roughly the same location and scale)? Or
do you mean theoretically justified knowledge?
>
> The only place integers should matter in Stan is (1.0 / 3.0) vs. (1 / 3).
>
> Chains typically get stuck in hierarchical models with wildly varying
> posterior
> scales.
I'm just trying something that may improve the model in this aspect.
Most parameters are in the roughly range between 0 and 1 in that model.
Two however are rather in the of 0 and 100. I try to bring the ranges
closer together by adding constant factors. I'm not sure though if that
is a good (or possible way) ...

> You can try reducing initial step size and increasing the
> target acceptance rate, which will slow down the model but make it more
> robust to getting stuck.
I have not came across yet how to control these sampling parameters.
I'll look it up in the manual.

Jan
phone: +49/ (0)5251/ 60 - 2907 Warburger Straße 100

Bob Carpenter

unread,
Oct 20, 2015, 6:14:13 PM10/20/15
to stan-...@googlegroups.com
Yes, fitting other models counts. It doesn't match the
prescriptive theory of Bayes (choose a prior, observe some
data, update to posterior), but it matches the practice of
model criticism and updating.

If you have one parameter in (0,1) and another in (0,100)
range, you probably want to bring them both down into the (-1/2, 1/2)
range or (-1,1) scale. (Not a hard constraint, just most of
the probability mass.) But that'd require multiplying and
shifting.

The way to control those parameters varies in name by interface,
but the doc of what they do is in the Stan language manual in the
part on algorithms.

- Bob

Andrew Gelman

unread,
Oct 20, 2015, 6:24:41 PM10/20/15
to stan-...@googlegroups.com

> On Oct 20, 2015, at 10:28 AM, Jan Tünnermann <je...@mail.upb.de> wrote:
>
>
> On 10/20/2015 04:17 PM, Bob Carpenter wrote:
>>
>>
>> Chains typically get stuck in hierarchical models with wildly varying posterior
>> scales.
> I'm just trying something that may improve the model in this aspect. Most parameters are in the roughly range between 0 and 1 in that model. Two however are rather in the of 0 and 100. I try to bring the ranges closer together by adding constant factors. I’m not sure though if that is a good (or possible way) ...
>

I think it might work better to multiply by a constant factor, rather than adding.

Jan Tünnermann

unread,
Oct 20, 2015, 6:29:17 PM10/20/15
to stan-...@googlegroups.com
Dear Bob,

thanks again for this hints.

On 10/21/2015 12:14 AM, Bob Carpenter wrote:
> Yes, fitting other models counts. It doesn't match the
> prescriptive theory of Bayes (choose a prior, observe some
> data, update to posterior), but it matches the practice of
> model criticism and updating.
There is another thing I have been thinking recently about: I would have
no problem obtaining a nice sample if I just would throw out the 1/4
chains that behave ugly, or so. E.g., fit 8 and keep the nicest 4. This
of course is fishy somehow in terms of representative results. It would
somehow solve the problem without really identifying what was going
wrong. But on the other hand I'm currently wildly trying stuff and may
also not know in the end what was going on (of course I really would
like to understand it). Either way, the final result is most likely the
same ...

> If you have one parameter in (0,1) and another in (0,100)
> range, you probably want to bring them both down into the (-1/2, 1/2)
> range or (-1,1) scale. (Not a hard constraint, just most of
> the probability mass.) But that'd require multiplying and
> shifting.
OK, I'll try that.
> The way to control those parameters varies in name by interface,
> but the doc of what they do is in the Stan language manual in the
> part on algorithms.
OK, thanks

Bob Carpenter

unread,
Oct 20, 2015, 6:50:13 PM10/20/15
to stan-...@googlegroups.com
The problem is that those non-converging chains may be indicating
that there's something biased in the sample you do get. When they
get stuck, it's usually because the sampler fell into a tightly
constrained neck of a funnel somewhere --- if you don't explore
the necks of the funnels, you get biased estimates. But if you
try to explore them with too-large step sizes, the sampler virtually
hangs, never being able to find the escape direction.

- Bob

Jan Tünnermann

unread,
Oct 20, 2015, 6:59:22 PM10/20/15
to stan-...@googlegroups.com
I'm using rstan and have passed "control=list(stepsize=0.01)" to the
stan command (after seeing this somewhere; hope it's right).
Also I'm not sure if 0.01 is a small step size ...

The problem is that I always need several runs to see if sometihng
changed because often the sample looks OK but the I get the bad one
after running it the third time or something.

Mabye running a lot (8 or 10?) of independent chains in one run would be
good at this point ...

Jan

Bob Carpenter

unread,
Oct 20, 2015, 7:02:19 PM10/20/15
to stan-...@googlegroups.com
RStan's set up to parallelize by default as of 2.8.0.

I believe 1 or 0.1 is the default, so 0.01 is a small-ish initial
step size. It's all problem dependent. We set things up assuming
all parameters were roughly normal(0, 1) distributed.

You can also look at where you're having problems in ShinyStan
(or even in some of the RStan output, I believe). If you're getting
stuck way out in the tail, then you might be able to initialize in
a more reasonable place and avoid the problem. But if you're getting
stuck in a volume of high posterior mass, then you might need
to reparameterize into a more favorable posterior geometry.

- Bob

Jan Tünnermann

unread,
Oct 20, 2015, 7:11:11 PM10/20/15
to stan-...@googlegroups.com
On 10/21/2015 01:02 AM, Bob Carpenter wrote:
> RStan's set up to parallelize by default as of 2.8.0.
>
> I believe 1 or 0.1 is the default, so 0.01 is a small-ish initial
> step size. It's all problem dependent. We set things up assuming
> all parameters were roughly normal(0, 1) distributed.
Apparently, 0.01 increases warmup but reduced sampling time (which I
guess is a good sign).
If I go further down it sampling takes longer again. So maybe it's a
good value. I think I will try the shifting and scaling into the -1/2
--- 1/2 range business first before further optimizing the step size.
> You can also look at where you're having problems in ShinyStan
> (or even in some of the RStan output, I believe). If you're getting
> stuck way out in the tail, then you might be able to initialize in
> a more reasonable place and avoid the problem. But if you're getting
> stuck in a volume of high posterior mass, then you might need
> to reparameterize into a more favorable posterior geometry.
Hmm, if I look at the trace plots, it is rather the second case (atleast
the chain gets stuck somewhere within the range where it normally mixes).

I already initialize the group level parameters close to the values I
expect. I have not explicitly initialized the subject level parameters,
though.

Thanks again

Jan Tünnermann

unread,
Oct 21, 2015, 7:01:01 PM10/21/15
to stan-...@googlegroups.com
Hi,

I'm currently busy scaling and shifting things as suggested below.

On 10/21/2015 12:14 AM, Bob Carpenter wrote:
> [...]
>
> If you have one parameter in (0,1) and another in (0,100)
> range, you probably want to bring them both down into the (-1/2, 1/2)
> range or (-1,1) scale. (Not a hard constraint, just most of
> the probability mass.) But that'd require multiplying and
> shifting.
>
> [...]
I have done it for some parameters now.

Apparently it slows down the sampling a lot. But the results look better
on the first sight (less auto-correlation; chain agree better).
Is this expected? Due to scaling and shifting, I could remove the
constraints for these parameters. Maybe this leads to the slowing but
better exploration?

Best regards
Jan

Bob Carpenter

unread,
Oct 22, 2015, 1:41:07 PM10/22/15
to stan-...@googlegroups.com

> On Oct 21, 2015, at 7:00 PM, Jan Tünnermann <je...@mail.upb.de> wrote:
>
> Hi,
>
> I'm currently busy scaling and shifting things as suggested below.
>
> On 10/21/2015 12:14 AM, Bob Carpenter wrote:
>> [...]
>>
>> If you have one parameter in (0,1) and another in (0,100)
>> range, you probably want to bring them both down into the (-1/2, 1/2)
>> range or (-1,1) scale. (Not a hard constraint, just most of
>> the probability mass.) But that'd require multiplying and
>> shifting.
>>
>> [...]
> I have done it for some parameters now.
>
> Apparently it slows down the sampling a lot. But the results look better on the first sight (less auto-correlation; chain agree better).
> Is this expected?

Yes, sometimes the stable parameters are slower than the
ones that don't work. Failing can be fast. Or it can be
slow --- depends on whether the step size vanishes to zero
(slow failure) or you get arithmetic errors (fast failure).

> Due to scaling and shifting, I could remove the constraints for these parameters. Maybe this leads to the slowing but better exploration?

You should have the constraints on parameters defined so
that any values meeting the constraint have nonzero density.

- Bob

Jan Tünnermann

unread,
Nov 3, 2015, 6:28:30 PM11/3/15
to stan-...@googlegroups.com
Dear all,

thanks to the help provided in this group, my model is now finally
converging properly.
In the end I re-parametrized 6 of the 7 group-level parameters with the
non-central form and roughly scaled them into similar ranges. For one
parameter this resulted in problems, so I left this one with in the
original version. Furthermore, I added some slightly(?) informative
priors for the group level parameters.

Here are diagnostics plots of ALL the parameters on all levels.

http://homepages.uni-paderborn.de/jeti/diag/

I think everything looks fine, doesn't it? Somewhere there is a bi-modal
posterior on the subject level, but I don't think this is any problem in
principle.

One thing remains where I am a bit unsure. I find it difficult to judge
how informative the priors are. I would like to keep it pretty
conservative, so I tried to only avoid extreme values. Furthermore, I
have an experimental and one neutral control condition. I used priors
(at least for the means) that reflect the neutral state (e.g., no shift
effect, equally distributed weights, etc.) and are quite diffuse for
both conditions.

I plotted the posteriors for parameters from each condition together
with a visualization of their common prior:

http://homepages.uni-paderborn.de/jeti/priors/

I'm very interested in your opinions regarding the informativeness of
these priors? From what I see the data clearly dominate the distributions.

And do you think using the "neutral" prior as described above for both
conditions is a conservative choice? As I'm interested in the
differences neutral-experimental, it could only reduce the detected effects.

Many thanks & best regards
Jan

Andrew Gelman

unread,
Nov 3, 2015, 9:06:00 PM11/3/15
to stan-...@googlegroups.com
I have an intuition that a little bit of prior information goes a long way—that is, using soft constraints to knock out truly unreasonable values will have payoff throughout the analysis because of near-collinearity and similar issues. The key is to knock out the unreasonable things. So if you have a parameter that really can’t be more than 1 or less than -1, then N(0,1) is conservative but can do the job. But N(0,10) is too weak to really regularize.
A

Jan Tünnermann

unread,
Nov 4, 2015, 2:51:52 AM11/4/15
to stan-...@googlegroups.com
Dear Andrew,

thanks for the reply. I'm not sure I understand. Is the problem that my
priors are essentially flat over the relvant range?
What would the mentioned payoff be (for example)?

Many thanks and best regards
Jan

Michael Betancourt

unread,
Nov 4, 2015, 8:06:53 AM11/4/15
to stan-...@googlegroups.com
“Noninformative” and “conservative” are loaded terms with little formal
meaning, so you’re probably not going to get a particularly satisfactory
answer.  When you model you have to make assumptions, whether
it be in the form of the likelihood and the prior (the hierarchical structure)
or the hyperpriors.  The best you can really judge is whether or not those
assumptions are qualitatively consistent with the given data which requires
doing posterior predictive checks.

Andrew is just saying that O(10) is probably still extreme if you have
standardized all of your data so that your random effects should be O(1).
Tightening the prior will improve your posterior inferences provided that
you don’t tighten it so far that you start introducing tension with the data.

Jan Tünnermann

unread,
Nov 4, 2015, 3:42:19 PM11/4/15
to stan-...@googlegroups.com
Hi Michael,

OK, I should be able to select tighter normal priors for the means. But how about the gamma priors on the standard deviations? Any hints on how to set them?  Up to now I have been adjusting shape and rate until I thought it looks visually OK. But maybe you or someone has hints how to better set them if I have some rough expectancy for them?

Thanks and best regards
Jan





Sent from Samsung tablet.

Jan Tünnermann

unread,
Nov 5, 2015, 5:31:23 PM11/5/15
to stan-...@googlegroups.com
Does anyone have a suggestion for how to find an intuition for the parameters of the gamma distribution? What I mean is, when I expect the parameter (a group standard deviation) around a certain value, how can I select a reasonable gamma prior for it?

Thanks and best regards
Jan

Krzysztof Sakrejda

unread,
Nov 5, 2015, 5:40:41 PM11/5/15
to Stan users mailing list
On Thursday, November 5, 2015 at 5:31:23 PM UTC-5, dqsIJan Tünnermann wrote:
Does anyone have a suggestion for how to find an intuition for the parameters of the gamma distribution? What I mean is, when I expect the parameter (a group standard deviation) around a certain value, how can I select a reasonable gamma prior for it?


curve(from=0, to=5, expr=dgamma(shape=2, rate=1, log=FALSE))

also, the wikipedia article on the gamma has expressions for the mean/median/variance which help.

Krzysztof

 

Jan Tünnermann

unread,
Nov 5, 2015, 6:18:42 PM11/5/15
to stan-...@googlegroups.com
Hi Krzysztof,

thanks for the reply. That is pretty much what I have been doing so far and it felt a bit unsystematic and I was hoping for some rule of thumb (how to set shape and rate in relation to the expected values and outliers, etc.). I thought there might be some reason why stan uses the shape / rate parametrization?

Best regards
Jan

Andrew Gelman

unread,
Nov 5, 2015, 6:23:20 PM11/5/15
to stan-...@googlegroups.com
Maybe no good reason to use a gamma anyway.

Jan Tünnermann

unread,
Nov 5, 2015, 6:30:08 PM11/5/15
to stan-...@googlegroups.com
I was thinking about switching back again to half-cauchy again, because I can control there shape more intuitively.
Any other suggestions are welcome!

Krzysztof Sakrejda

unread,
Nov 5, 2015, 7:23:17 PM11/5/15
to Stan users mailing list
Jan, rules of thumb:

1) rate is just 1/scale
2) mean ~ 1/rate when shape > 1, stick to shape = 2 unless you need a tighter tail on the high end
3) the mass for rates ~ 2 is mostly between zero and 2 or 3 * scale
4) the slope at zero is positive for shape > 1 (?) so essentially it weakly avoids that edge... not sure how this compares to half-cauchy with zero-bounded parameter

It would have been nice if all Stan distributions went either with scale or rate parameterization but....

Krzysztof

Ben Goodrich

unread,
Nov 5, 2015, 7:50:03 PM11/5/15
to Stan users mailing list
On Thursday, November 5, 2015 at 7:23:17 PM UTC-5, Krzysztof Sakrejda wrote:
It would have been nice if all Stan distributions went either with scale or rate parameterization but....

Well, we mostly went with BDA parameterizations which often followed Bayesian conventions and Bayesian conventions were usually designed to facilitate Gibbs sampling. And as we went along, we came to realize how important it is to establish a reasonable scale for parameters when using the default sampler in Stan. So, it probably would have been better to use scale parameterizations from the outset. But for priors on unknown quantities, it is typically a good idea to define a primitive in the parameters block whose prior entails scale = 1 = rate and than multiply that primitive by your prior scale inside the transformed parameters block to get the quantity that corresponds to the algebraic representation of the log-likelihood.

Ben

 

Andrew Gelman

unread,
Nov 5, 2015, 8:11:24 PM11/5/15
to stan-...@googlegroups.com
No, Bayesian conventions predate Gibbs sampling.  I think the conventions come from various sources; different stories for different distributions.  One awkward example we’ve discussed on this list is that the exponential distribution is parameterized in terms of rate, and the Laplace double-exponential distribution is parameterized in terms of scale = 1/rate.

Ben Goodrich

unread,
Nov 5, 2015, 10:34:03 PM11/5/15
to Stan users mailing list, gel...@stat.columbia.edu
On Thursday, November 5, 2015 at 8:11:24 PM UTC-5, Andrew Gelman wrote:
No, Bayesian conventions predate Gibbs sampling.  I think the conventions come from various sources; different stories for different distributions.  One awkward example we’ve discussed on this list is that the exponential distribution is parameterized in terms of rate, and the Laplace double-exponential distribution is parameterized in terms of scale = 1/rate.

What is Bayesian justification for using the rate parameterization for the exponential if not for the fact that it is conjugate with the rate-parameterized gamma distribution?


Andrew Gelman

unread,
Nov 5, 2015, 10:45:39 PM11/5/15
to Ben Goodrich, Stan users mailing list
I don’t know, it’s just the way they always wrote it in the textbooks.  When I learned probability theory in 1981, my textbook had the exponential distribution and it was parameterized according to the rate parameter.  Nothing Bayesian about it, it’s just the convention.

Ben Goodrich

unread,
Nov 6, 2015, 12:02:44 AM11/6/15
to Stan users mailing list, gel...@stat.columbia.edu
On Thursday, November 5, 2015 at 10:45:39 PM UTC-5, Andrew Gelman wrote:
I don’t know, it’s just the way they always wrote it in the textbooks.  When I learned probability theory in 1981, my textbook had the exponential distribution and it was parameterized according to the rate parameter.  Nothing Bayesian about it, it’s just the convention.

I would agree that the frequentists had conventions about parameterizations prior to Gibbs sampling but, at least for maximum likelihood, the frequentists didn't have a frequentist justification for one parameterization over another due to the invariance property. Substantively, in the context of a duration model, it could either make sense to think of events occurring at some rate or to think about the mean of an exponential process.

But I don't think substantive reasons should be that relevant to how Stan parameterizes distributions because a user could just declare a rate parameter in the parameters block, specify an informative prior on the rate in the model block, and flip it over when putting into a likelihood that is parameterized in terms of scale (or vice versa). For technical reasons, it seems that HMC is more suited to scale (re)parameterizations, but those parameterizations tend to run counter to the parameterizations Bayesians were using for technical reasons when they were developing Gibbs samplers.

Ben
 

Krzysztof Sakrejda

unread,
Nov 6, 2015, 6:11:31 AM11/6/15
to Stan users mailing list, gel...@stat.columbia.edu

I'm pretty sure this just comes from the failure rate literature. Krzysztof
 
Reply all
Reply to author
Forward
0 new messages