R/bayes_util.R

In statsr: Companion Software for the Coursera Statistics with R Specialization

BF_plot = function(den_H2, res, parameter)
{
  d_H2 = data.frame(x = den_H2$x, 
                    y = den_H2$y * res$post_H2 / max(den_H2$y),
                    Hypothesis = "H2") 

  li = min(which(d_H2$x >= res$ci_H2[1]))  
  ui = max(which(d_H2$x <  res$ci_H2[2]))

  d_H2_poly = data.frame(x = c(d_H2$x[c(li,li:ui,ui)]), 
                         y = c(0, d_H2$y[li:ui], 0),
                         Hypothesis = "H2")

  d_H1 = data.frame(x = c(res$null, res$null), 
                    y = c(0, res$post_H1),
                    Hypothesis = "H1")

  d = rbind(data.frame(x=NA, y=NA, Hypothesis="H1"),
            d_H2,
            data.frame(x=NA, y=NA, Hypothesis="Overall"))

  # H2 Features
  p = ggplot(d, aes_string(x="x", y="y", color="Hypothesis", fill="Hypothesis")) + 
      geom_line(alpha=0.8) +
      geom_polygon(data = d_H2_poly, linetype="blank",alpha=0.8) + 
      ylab("Density") +  
      xlab(parameter)

  # H2 Features
  p = p + geom_line(data = d_H1, size=1.5, alpha=0.8)
      

  # Marginal plot features 
  y_min = ggplot_build(p)$panel$ranges[[1]]$y.range[1]

  d_Marg = data.frame(x = rep(res$ci_Marg, c(2,2)),
                      y = c(y_min*1/2, y_min, y_min, y_min*1/2),
                      Hypothesis = "Overall")

  p = p + geom_line(data = d_Marg, size=0.75, alpha=0.8)

  print(p)
}

coda_density = function(x, from, to)
{
    bwf = 1.06 * min(sd(x), IQR(x)/1.34) * length(x)^-0.2
    
    return(density(x, from=from, to=to, bw=bwf))
}


check_beta_prior = function(beta_prior, group="")
{
    arg_name = paste(substitute(beta_prior))
    if (arg_name == "") arg_name = "beta_prior"

    param = ifelse(group == "", "p", paste0("p_",group))    
    
    if (is.null(beta_prior))
    {
        warning("No beta prior for ",param," was specified, assuming a uniform prior (p ~ Beta(a=1,b=1)).\n",
                "  This beta prior is specified using the argument ",arg_name,"=c(a,b),\n",
                "  where a and b are your desired hyperparameters.")
        beta_prior = c(a=1,b=1)
    }
    
    stopifnot(length(beta_prior) == 2)
    
    if (is.null(names(beta_prior)))
        names(beta_prior) = c("a","b")
    stopifnot(all(sort(names(beta_prior)) == c("a","b")))
    beta_prior = beta_prior[c("a","b")]
    
    return(beta_prior)
}

check_hypothesis_prior = function(prior)
{
    if (is.null(prior))
    {
        warning("No prior set for H1 and H2, assuming a uniform prior of P(H1) = 0.5 and P(H2) = 0.5. The hypothesis prior is assigned using the argument  prior=c(H1=a,H2=b). ")
        prior = c(H1=0.5,H2=0.5)
    }

    if (length(prior) == 1)
    {   
        if (names(prior) %in% c("H1","H2"))
            prior[ setdiff(c("H1","H2"), names(prior)) ] = 1 - prior
    }

    stopifnot(length(prior) == 2)
    stopifnot(all(prior >= 0))
    stopifnot(sum(prior) == 1)

    if (is.null(names(prior)))
        names(prior) = c("H1","H2")

    stopifnot(all(sort(names(prior)) == c("H1","H2")))
    
    return(prior[c("H1","H2")])
}