Bernoulli model, Normal prior

```{r, echo=FALSE}
inputPanel(
  sliderInput(“p”, label = “Probability (true)”,
              min = 0, max = 1, value = 0.3, step = 0.05),
  sliderInput(“n”, label = “Sample size”,
              min = 0, max = 1000, value = 10, step = 10),
  sliderInput(“mu_1”, label = “Normal prior mean”,
              min = -10, max = 10, value = -2, step = 1),
  sliderInput(“sig2_1”, label = “Normal prior variance”,
              min = 0.001, max = 10, value = 1, step = 1),
  sliderInput(“mu_2”, label = “Normal prior mean”,
              min = -10, max = 10, value = 2, step = 1),
  sliderInput(“sig2_2”, label = “Normal prior variance”,
              min = 0.001, max = 10, value = 2, step = 1),
  radioButtons(“scale”, label=”Scale”,
             c(“Probability (0, 1)” = “prob”,
               “Logit (-Inf, Inf)” = “logit”),
             selected = “logit”),
  sliderInput(“seed”, label = “Random seed”,
              min = 0, max = 100, value = 0, step = 10)
)
renderPlot({
    par(las = 1)
    set.seed(input$seed)
    y <- rbinom(n = 1000, size = 1, p = input$p)
    BY <- 0.0005
    pval <- seq(0.001, 0.999, by = BY)
    fLik <- function(p, y)
        prod(dbinom(y, size = 1, prob = p))
    Lik <- sapply(pval, fLik, y=y[1:input$n])

fPri <- function(p, mu_1, sig2_1, mu_2, sig2_2)
    0.5 * (dnorm(qlogis(p), mu_1, sqrt(sig2_1)) + 
        dnorm(qlogis(p), mu_2, sqrt(sig2_2)))
Pri <- sapply(pval, fPri, input$mu_1, input$sig2_1, 
    input$mu_2, input$sig2_2)
if (input$scale == "prob") {
    p <- input$p
    br <- qlogis(c(0.001, seq(0.001+BY/2, 0.999-BY/2, by = BY), 0.999))
    dx <- diff(qlogis(pval))
    dx <- c(dx[1], dx)
    d <- Pri * dx / diff(plogis(br))
    Pri <- smooth.spline(pval, d)$y
} else {
    pval <- qlogis(pval)
    p <- qlogis(input$p)
}
    
Pos <- Lik * Pri
M <- cbind(Pri=Pri/max(Pri),
    Lik=Lik/max(Lik),
    Pos=Pos/max(Pos))
Col <- c("#cccccc", "#3498db", "#f39c12")
matplot(pval, M, type = "l", 
    col=Col, lwd=2, lty=1,
    ylab = "Density", 
    xlab=ifelse(input$scale == "logit", "logit(p)","p"),
    sub=paste0("Mean = ", round(mean(y[1:input$n]), 2), " (", 
        sum(1-y[1:input$n]), " 0s & ", sum(y[1:input$n]), " 1s)"),
    main = paste0("True value = ", round(p, 2), 
        ", Posterior mode = ", round(pval[which.max(Pos)], 2)))
abline(v = p, lwd = 2, col = "#c7254e")
abline(v = pval[which.max(Pos)], lwd = 2, col = "#18bc9c")
legend("topleft",lty=1, lwd=2, col=Col, bty="n",
       legend=c("Prior","Likelihood","Posterior")) }) ```