Bernoulli model, data cloning

{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 = 1, max = 50, value = 10, step = 5),
  sliderInput("a", label = "Beta prior shape parameter a",
              min = 0, max = 2, value = 1, step = 0.1),
  sliderInput("b", label = "Beta prior shape parameter b",
              min = 0, max = 2, value = 1, step = 0.1),
  sliderInput("K", label = "Number of clones",
              min = 1, max = 100, value = 1, step = 10),
  radioButtons("scale", label="Scale",
             c("Probability (0, 1)" = "prob",
               "Logit (-Inf, Inf)" = "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 = input$n, size = 1, p = input$p)
    yk <- rep(y, input$K)
    BY <- 0.0005
    pval <- seq(0.001, 0.999, by = BY)
    fLik <- function(p, y)
        sum(dbinom(y, size = 1, prob = p, log=TRUE))
    Lik <- exp(sapply(pval, fLik, y=yk))
    if (all(Lik <= 0)) {
        est <- optimize(fLik, c(0.001, 0.999), y=yk, maximum=TRUE)$maximum
        Lik[which.min(abs(pval - est))] <- 1
    }
    if (input$scale == "prob") {
        p <- input$p
        fPri <- function(p, shape1=0.5, shape2=0.5)
            dbeta(p, shape1, shape2)
        Pri <- sapply(pval, fPri, input$a, input$b)
    } else {
        p <- qlogis(input$p)
        N <- 10^5
        x <- rbeta(N, input$a, input$b)
        br <- c(0.001, seq(0.001+BY/2, 0.999-BY/2, by = BY), 0.999)
        d <- as.numeric(table(cut(x, breaks=br))) / N
        pval <- qlogis(pval)
        g <- diff(qlogis(br))
        gy <-  d / g
        Pri <- smooth.spline(pval, gy)$y
    }
    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="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"))
})