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.
