Bayesian inferenceExample 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 |
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"))
)
))