Hot answers tagged

19

Firstly, feel free to ask questions like this on our users' list (http://mc-stan.org/mailing-lists.html) where we discuss not only issues related to Stan implementations/optimizations/etc but also practical statistical and modeling questions. As to your question, it's absolutely a fine approach. There are many ways to justify it more rigorously (for ...


16

Quick thoughts: 1) The key issue is what applied question you are trying to answer for your audience, because that determines what information you want from your statistical analysis. In this case, it seems to me that you want to estimate the magnitude of differences between groups (or perhaps the magnitude of ratios of the groups if that is the measure ...


14

A divergent transition in Stan tells you that the region of the posterior distribution around that divergent transition is geometrically difficult to explore. For example here is a quote from the manual: The primary cause of divergent transitions in Euclidean HMC (other than bugs in the code) is highly varying posterior curvature, for which small step ...


13

This is of course a diverse set of people with a range of opinions getting together and writing a wiki. I summarize I know/understand with some commentary: Choosing your prior based on computational convenience is an insufficient rationale. E.g. using a Beta(1/2, 1/2) solely because it allows conjugate updating is not a good idea. Of course, once you ...


12

Because all of the parameters of this distribution are known, and we merely want to draw samples from this distribution, coding the model in rstan is straightforward. Note that this is, by far, one of the least efficient paths to sampling from this particular model, in terms of time that I spent coding it (15 minutes). The author of the original post is ...


11

I would distinguish between the prior distribution and the parameters for the prior distribution. When I started with Stan, I would set the parameters to the prior distributions just as some values. So in the model step, I would have something like model { mu ~ normal(0, 1) y ~ normal(mu, s) } for a normal prior on the mean coefficient for ...


9

You can define a proper or improper prior in the Stan language using the increment_log_prob() function, which will add its input to the accumulated log-posterior value that is used in the Metropolis step to decide whether to accept or reject a proposal for the parameters. In your example, the model block would need to include the new line ...


9

Hamiltonian Monte Carlo is often quoted as being rotation invariant, but what does that actually mean? In theory Hamiltonian Monte Carlo (HMC) is independent of the coordinates chosen to represent your distribution, which actually means it is both rotation and scale invariant. The problem is that, in practice, we can't run HMC exactly and must instead use ...


9

If you have a disjoint set of possible events (say three of them) that exhaust all possiblites, the probabilities of each of these events must sum to one: $$p_1 + p_2 + p_3 = 1$$ Being probabilities, each of these is also bounded between zero and one: $$ 0 \leq p_i \leq 1 $$ If consider the three probabilities as a point in euclidean space $(p_1, p_2, ...


9

Objects declared in the transformed parameters block of a Stan program are: Unknown but are known given the values of the objects in the parameters block Saved in the output and hence should be of interest to the researcher Are usually the arguments to the log-likelihood function that is evaluated in the model block, although in hierarchical models the line ...


9

Fitting random slopes with the population-level slope fixed to zero is not out of the question - it's not mathematically or statistically ill-posed - but it's a rather weird model that would require some extra justification. Why would you expect that the average slope across cities would be exactly zero (which is what is implied by the model that omits the ...


8

Keep in mind that in Bayesian inference, you are looking at the entire range of values in the posterior, rather than point-estimate summaries. After convergence, MCMC samples are samples drawn from that posterior density. Using just the latest sample is the same as picking a random point from the posterior and using that as your sample estimate. For an ...


8

Likelihood For a mixture of two Gaussians, the likelihood can be written as: $$ y_i \sim \pi N(y_i|\alpha_0 + x_i\beta, \sigma_0) + (1-\pi) N(y_i|\alpha_1 + x_i\beta, \sigma_1) $$ where $\pi \in [0, 1]$. This is fine, but having two components in the likelihood makes sampling more difficult. A trick when dealing with mixture models is to augment the ...


8

They do not provide any scientific/mathematical justification for doing so. Most of the developers do not work on this kind of priors, and they prefer to use more pragmatic/heuristic priors, such as normal priors with large variances (which may be informative in some cases). However, it is a bit strange that they are happy to use PC priors, which are based ...


7

$\hat{R}$ and "potential scale reduction factor" refer to the same thing. See Chapter 6 of the Handbook of Markov Chain Monte Carlo, "Inference from Simulations and Monitoring Convergence" by Andrew Gelman and Kenneth Shirley. In Stan, the number reported is actually split $\hat{R}$; the calculation of $\hat{R}$ is computed with each of the chains split in ...


7

The primary reason that your code does not yield the expected answer is that you are using the multi_normal_prec likelihood rather than the multi_normal likelihood. The former expects a precision matrix (the inverse of a covariance matrix) as its second argument, while the latter expects a covariance matrix. For what it is worth, you should be able to ...


7

As my previous answer was deleted, here is a more explicit one, with an example using a simple sampling from the prior: library(rstan) model = " parameters { real p; } model { p ~ normal(1,3); } " fit = stan(model_code = model, pars = c('p'), control=list(adapt_delta=0.99, max_treedepth=10), iter = 5000, chains = 1, ...


7

As mentioned by previous answers, Stan, JAGS, and WinBUGS require that priors be specified as mathematical functions. If you've already got an MCMC-represented posterior from a previous analysis, and you want to use that MCMC posterior as a prior for subsequent data, you must approximate the MCMC posterior in a mathematical form. Unless you have a simple ...


7

This is a famous problem known as Gull's Lighthouse, from an example by Gull in 1988. It has deep implications when taken one additional step in both the social sciences and in physics. You actually have enough information to solve this problem through acceptance-rejection testing, but if you want to use MCMC feel free. Let's look at your problem another ...


7

In No-U-Turn-Sampler a maximum tree depth of 10 is a sensible default, but occasionally you have to increase it. In my experience not usually by much. I might try 12 next and I have never had to go beyond 15, so far. As it increases computation time, it is not a good trade-off to make it higher by default - unless you always run into this problem in a ...


7

I'll leave this as an "answer" as I don't have enough reputation to "comment" on this post. This webpage might be of interest to you. The development team describses here, although quite shortly, the main implications and solutions to warnings in STAN. Reaching the Maximum treedetph is, according to their explanation, far less an issue than divergent ...


7

Each value of y_miss can either be 0 or 1, so you need to marginalize over them with a statement such as for (i in 1:n_miss) { real eta = a + beta * x_miss[i]; target += log_mix(inv_logit(eta), bernoulli_logit_lpmf(1 | eta), bernoulli_logit_lpmf(0 | eta)); } where x_miss is a vector of observations on the predictor ...


6

The posterior predictive distribution of a logit model involves four conceptual steps: Draw from the posterior distribution of the parameters Create linear predictors using those parameter draws Transform those linear predictors into probabilities via the inverse link function Draw from a Bernoulli (or Binomial with size = 1) distribution with those ...


6

The main question was the result of a bug in rstanarm that has since been fixed on GitHub. However, in general, we do not recommend rstanarm models that exclude the intercept. A better alternative is to place a tight prior with mean zero on the intercept. In this case, including the intercept yields a better fit and similar results between lme4::lmer and ...


6

First stan_lm only accepts one type of prior: Must be a call to R2 with its location argument specified or NULL, which would indicate a standard uniform prior for the R^2. If you want to specify any of the other types of priors, you'll need to switch to stan_glm. With the default family (gaussian), it will be equivalent to stan_lm. Second, you've tried to ...


6

Let's work it out from first principles, beginning with the hard work of computing a convolution. As an auxiliary calculation, consider the distribution of $W=X+Y$ where $Y$ has an Exponential distribution with pdf $$f_Y(y) = e^{-y}\,\mathcal{I}(y\gt 0)$$ and $X$ has a Normal$(\mu,\sigma^2)$ distribution with pdf $f_X(x;\mu,\sigma) = \phi((x-\mu)/\sigma)/\...


6

One of the more interesting choices in R is rstan, where you could code this up yourself in the Stan modeling language (which tends to be amazing in that it can produce inference for models that we used to be unable to do for a long time). However, getting started can be a little challenging and it sounds like you'd like a higher level interface. That could ...


5

From the manual, 2.8.0, page 101, just swap the normal function for Weibull, "One way to model censored data is to treat the censored data as missing data that is constrained to fall in the censored range of values. Since Stan does not allow unknown values in its arrays or matrices, the censored values must be represented explicitly, as in the ...


5

In BUGS/JAGS, the order in which statements are written does not matter. In Stan statements execute in the order in which they are written (see Stan 2.2.0 Reference Manual, pg. 405). Thus your last statement is in the wrong place: lambda is sampled from a gamma distribution, but that happens after the previous statements, so it's sampled 'in the air'. ...


5

$y \mid p,\lambda$ is Poisson! Your marginalization, or at least the end result, is correct. The form you have obtained for the distribution is the probability mass function of a Poisson distribution -- just write $p^y\lambda^y$ as $(\lambda\,p)^y$ and behold. That is, \begin{equation} y \mid p, \lambda \sim \mathrm{Poisson}(\lambda\,p). \end{equation} ...


Only top voted, non community-wiki answers of a minimum length are eligible