Bayesian inference

Example E of section 3.5 of Rice 3rd edition, page 94-95

beta prior (blue), likelihood is binomial, posterior is beta (green), plus likelihood (brown)

See Figure 3.16 for results from a beta(1, 1) prior and 13 successes out of 20 attempts

Prior is blue, posterior is green

show with app
options(bitmapType='cairo')   # only needed temporarily

shinyServer(function(input, output) {
  output$distPrior <- renderPlot({
 
    require(mosaic)
    require(latticeExtra)
    trellis.par.set(theme=theme.mosaic())
    prior <- list(input$alpha, input$beta)
    posterior <- list((input$alpha + input$y), (input$beta + input$n - input$y))

    # Watch out for Michael!
    if (input$n < input$y) {
       stop("Number of successes must be at least as large as the number of trials!")
    }

    # pesky details to tweak y axis limits
    xvals <- seq(0, 1, length=101)
    eps <- 0.10
    maxy1 <- max(dbeta(xvals, prior[[1]], prior[[2]]))
    maxy2 <- max(dbeta(xvals, posterior[[1]], posterior[[2]]))
    maxy <- min(8, max(maxy1, maxy2)+eps)

    # pesky detail to rescale the likelihood
    maxbinom <- dbinom(input$y, input$n, input$y/input$n)
    mle <- makeFun(0.99*maxy/maxbinom*choose(n, y)*p^y*(1-p)^(n-y) ~ p, n=input$n, y=input$y)

    # display the results
    plot1 <- plotDist('beta', params=prior, 
                     ylim=c(0, maxy), xlim=c(0, 1), lwd=2)
    plot2 <- plotDist('beta', params=posterior, lwd=2, col="green")
    plot3 <- plotFun(mle(p) ~ p, lwd=2, col="brown")
    tmpplot <- plot1 + as.layer(plot2)
    tmpplot + as.layer(plot3)
    })
  }
)
library(shiny)

shinyUI(fluidPage(
  # Application title
  titlePanel("Bayesian inference"),
  p("Example E of section 3.5 of Rice 3rd edition, page 94-95"),
  p("beta prior (blue), likelihood is binomial, posterior is beta (green), plus likelihood (brown)"),
  p("See Figure 3.16 for results from a beta(1, 1) prior and 13 successes out of 20 attempts"),

  sidebarLayout(
    sidebarPanel(
      p("Prior is blue, posterior is green"),
      numericInput("alpha", "alpha parameter for Beta prior:", 1,
                   min = 1, max = 100),
      numericInput("beta", "beta parameter for Beta prior:", 1,
                   min = 1, max = 100),
      numericInput("y", "number of successes observed:", 13,
                   min = 1, max = 100),
      numericInput("n", "out of this many trials (must be larger than # successes):", 20,
                   min = 1, max = 100)),
    # Show a plot of the generated distribution
    mainPanel(plotOutput("distPrior"))
  )
))