When trying to find the suitable priors for my models, I’d like to quickly check how the prior would look like by plotting it e.g.
hist(rskew_normal(n = n, mu = 4, sigma = 0.5, alpha = -10)) but I found that the distribution that I create like this and the distribution that I get when running a
brms model with
sample_prior = "only" are not the same when I want a skew normal distribution. Am I setting the paramters wrong? Are they different? Do they priors influence each other some how?
Here is an example:
and the reproducible code:
# Library library(brms) library(ggplot2) # Seed set.seed(20230307) # Generate data with rskew_normal n <- 4000 dists <- data.frame(Value = rskew_normal(n = n, mu = 4, sigma = 0.5, alpha = -10), Method = rep("rskew_normal", n)) # Create a random data set to fit a model fake_data <- data.frame(x = rnorm(100), y = rnorm(100)) # Set brms priors test_prior <- c(prior(skew_normal(4, 0.5, -10), class = "b")) # Fit brms model currentModel <- brm(y ~ x, data = fake_data, prior = test_prior, sample_prior = "only", seed = 1234, backend = "cmdstanr", silent = 2, refresh = 0) # Extract the prior from the fit dists <- rbind(dists, data.frame(Value = as_draws_df(currentModel)$b_x, Method = rep('brms_fit', ndraws(currentModel)))) # Visualise the difference ggplot(dists, aes(x = Value, fill = Method)) + geom_density(alpha = 0.5) + labs(title = "Difference between rskew_normal & skew_normal")
For comparison, I also did the same for the normal distribution, which looks like expected.