prior_predictive_from_expert.R

set.seed(42) # The answer to life, the universe and everything

## This is what the expert said
p10 <- 8
p50 <- 12
p90 <- 25

# I can approximate it with Myerson distribution defined by quantiles
prior_predictive_elicited <- tidyear::rMyerson(1e5, p10, p50, p90, tl=0.2)

summary(prior_predictive_elicited)
#  Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
# 6.332   9.320  12.016  15.090  17.014 315.607 
quantile(prior_predictive_elicited, c(0.1, 0.5, 0.9))
#      10%       50%       90% 
# 7.991898 12.016108 25.106989

# Trying to pick parameters
mu <- rlnorm(1e5, log(12), log(1.5))
sigma <- rlnorm(1e5, log(1.5), log(1.1))

# That would produce the same prior predictive as the expert said
prior_predictive_fitted <- rlnorm(1e5, log(mu), log(sigma))

summary(prior_predictive_fitted)
#   Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
# 0.8165   8.1428  11.9824  14.1638  17.5939 183.0402 
quantile(prior_predictive_fitted, c(0.1, 0.5, 0.9), na.rm = TRUE)
#      10%       50%       90% 
# 5.699318 11.982383 25.013896 

library(hrbrthemes)

ggplot()+
  geom_density(aes(prior_predictive_elicited), col="blue")+
  geom_density(aes(prior_predictive_fitted), col="red")+
  scale_x_continuous(limits = c(0,75))+
  theme_ipsum_rc(grid_col = "gray95")+
  labs(x="prior predictive")

ggsave("prior_predictive.png")

Probably lognormal mu and lognormal sigma are not suitable distributions to emulate this prior predictive. Perhaps inverse gamma would help? Reality is that all of these "advaced" distributions are really difficult to relate to - especially when you need to go back to the expert and explain how his input has "moved the needle".