Weird output for a very simple model with sample_prior = "only"

Hi,

I’m fitting a model to generate prior predictive distributions, but they are not looking as expected. (I know that the priors are terrible, but that’s the whole point.)

library(dplyr)
library(brms)
df_spacebar_empty <- tibble(rt = runif(361, 0, 10000))
fit_prior_press <- brm(rt ~ 1,
  data = df_spacebar_empty,
  family = gaussian(),
  prior = c(
    prior(uniform(0, 60000), class = Intercept),
    prior(uniform(0, 2000), class = sigma)
  ), chains = 1,
  sample_prior = "only")

First of all I would expect the intercept to be 60000/2, but it’s 2604.49. (I got even smaller numbers when I set warmup =0, which shouldn’t be a problem, or am I wrong?)

fit_prior_press
(...)
Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept  2604.49   1511.28   622.08  6442.92 1.53        2       11

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma   997.40    559.00    51.91  1929.68 1.00       93      130
(...)
Warning messages:
... Ton of warnings here

I checked the prior predictive distributions, and I have the same problem:

posterior_predict(fit_prior_press) %>% 
  rowMeans() %>% 
  mean()
2606

Am I misunderstanding what brms is supposed to do? I’ve checked the stan code and it’s what I expect, but I’m not sure how brms is generating the predictive distribution.

  • Operating System: Ubuntu 20.04
  • brms Version: 2.13.5

The model does not like as if it has converged if you check the ESS and Rhat.

ok, that’s what I don’t get. If it’s not fitting to any data, how can it not converge?

Because it still samples the parameters using hmc if sample_prior = “only”

oh, is there any way to workaround this? I mean to just sample from the priors, as if the priors are also in the generated quantities block in Stan. (I know brms doesn’t use it, but something similar).

Well, if sample_prior = "yes" priors are sampled from the in generated quantities block but only additionally to the posterior samples and not for the intercept (unless you use bf(..., center = FALSE)).

sampled from the in generated quantities block but only additionally to the posterior samples
I didn’t understand that.

If you have target += normal_lpdf(beta,0,1);
then the sampler will give values of beta, how can I also have them from the generated quantities? As in
beta = normal_rng(0,1)

And by the way, is bf(..., center = FALSE), the same as doing 0 + intercept +....? (If so, is the use of bf, the preffered one)?

You can only have one of the two beta options at the same time. Practically 0 + Intercept is the same as center = FALSE

Yes, that’s what I don’t understand.

Say here:

fit_prior_press <- brm(rt ~ 1,
  data = df_spacebar_empty,
  family = gaussian(),
  prior = c(
    prior(uniform(0, 60000), class = Intercept),
    prior(uniform(0, 2000), class = sigma)
  ), chains = 1,
  sample_prior = "only")

Is, for example, here sigma taken from the samples that the HMC finds, or is it generated? (I guess not in the generated quantities, because it’s empty, right? If it’s generated it has to be in R)

Correct, it is sampled via HMC in this case.

Sorry to keep going with this, I’m trying to understand exactly what’s going on

I’ve tried this

fit_prior_press <- brm(bf(rt ~ 1, center = FALSE),
  data = df_spacebar_empty,
  family = gaussian(),
  prior = c(
    prior(uniform(0, 60000), class = b, lb = 0, ub =60000),
    prior(uniform(0, 2000), class = sigma)
  ), chains = 1,
  sample_prior = "only")

and now the values make sense. But I still get divergent transitions, I guess it’s because of this

  target += uniform_lpdf(b | 0, 60000)
    - 1 * log_diff_exp(uniform_lcdf(60000 | 0, 60000), uniform_lcdf(0 | 0, 60000));
  target += uniform_lpdf(sigma | 0, 2000)
    - 1 * uniform_lccdf(0 | 0, 2000);

I think we are talking past to each other.
Let me verify if this correct:
Even if sample_prior = “yes”, brms takes samples from the “posterior” distribution of b and sigma based on HMC and generates (using some RNG) the prior predictions based on that, right?

brms never uses a RNG for the parameters, I mean something like the following is never done

b = rnunif(ndraws, 0,2000)
sigma = rnunif(ndraws, 0,60000)

That means that I can never ignore the warnings, is that right?

Yes, to all of your questions.

1 Like

ok, now I’m satisfied :)

I think it would be nice if everything could be generated with RNGs when we sample from the priors. But I’m aware that this would be a ton of work.

1 Like

I understand, but as you say, this is a lot of work, and it is unlikely that I will be able to invest the time into this feature.

Oh, and in addition, such a feature, even if implemented would not fully support all priors are not all priors can be sample from directly via _rng.

ok, yes I know. And also for mixture models, and custom likelihoods, etc…
Anyways, thanks for the help, now I understand what’s going on.

In my defense, I wasn’t the only one confused by this:

I fully understand your confusion. The shortness of my answers is just related to my time budget for discourse right now so that I don’t spend too much time answering questions every day :-D

1 Like