## 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

``````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[], prior[]))
maxy2 <- max(dbeta(xvals, posterior[], posterior[]))
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"))
)
))``````