bayes tails examples

bayes_tails.R

## original by Richard McElreath at https://twitter.com/rlmcelreath/status/1701165075493470644, 
##   https://gist.github.com/rmcelreath/39dd410fc6bb758e54d79249b11eeb2f
## originally based on https://doi.org/10.1006/jmva.1999.1820
## with improvements from https://gist.github.com/murphyk/94205bcf335837108ff0e6d51331785a

post_prior <- function(
                       m_pr = 10, m_lik = 0,        ## mean values for prior and likelihood
                       sd_pr = 1, sd_lik = 1,       ## standard deviations
                       df_pr = 1000, df_lik = 1000, ## dfs (df = 1000 is effectively Gaussian)
                       xlim = c(-5, 15), n = 1001,   ## range and delta for x-vector
                       ymin = NA,
                       lab_offset = 0.1,
                       ...)           ## additional args to matplot()
{
    xvec <- seq(xlim[1], xlim[2], length.out = n)
    ## scaled/shifted t distribution (dt() only does standardized t)
    tfun <- function(x,m,s,df, ...) dt((x-m)/s, df)/s
    ## compute prior & likelihood curves
    prvec <- tfun(xvec, m_pr, sd_pr, df_pr)
    likvec <- tfun(xvec, m_lik, sd_lik, df_lik)
    ## posterior (Bayes rule!)
    postvec <- prvec*likvec/(sum(prvec*likvec)*diff(xvec[1:2]))
    ylim <- c(max(ymin, min(pmin(prvec, likvec)), na.rm = TRUE),
              max(prvec, likvec, postvec)*(1+lab_offset))
    ## plot it all
    matplot(xvec, cbind(prvec, likvec, postvec),
            type = "l",
            lty = c(2,1,1),
            col = c(1,1,2),
            lwd = c(1,1,2),
            ylab = "density",
            xlab = "",
            ylim = ylim, ...)
   maxval <- max(prvec, likvec)*(1+lab_offset)
   text(m_lik, maxval, "likelihood")
   text(m_pr, maxval, "prior")
}

par(mfrow=c(2,2), las = 1, bty ="l", mar = c(3,4,2,2)+0.1)
post_prior(main="prior:N, lik:N") # log = "y")
post_prior(df_lik = 2, df_pr =2, main="prior:T, lik:T")
post_prior(df_pr = 2, df_lik =1000, main="prior:T, lik:N")
post_prior(df_pr = 1000, df_lik = 2, main="prior:N, lik:T")

par(mfrow=c(2,2), las = 1, bty ="l", mar = c(3,4,2,2)+0.1)
post_prior(main="prior:N, lik:N", log = "y", ymin = 1e-10)
post_prior(df_lik = 2, df_pr =2, main="prior:T, lik:T", log = "y", ymin = 1e-10)
post_prior(df_pr = 2, df_lik =1000, main="prior:T, lik:N", log = "y", ymin = 1e-10)
post_prior(df_pr = 1000, df_lik = 2, main="prior:N, lik:T", log = "y", ymin = 1e-10)