rm(list=ls())

require(rpanel)

rpBeta <- function() {
    plotbeta <- function(panel) {
        with(panel,
             plot(1, type = "n", xlim = c(0, 1), ylim = c(0, 7.5),
                  xlab = "Parâmetro", ylab = "Densidade", 
                  main = paste("Priori: Beta(", alpha, ",", beta, ")", sep="")
                  ))
        x <- seq(0, 1, length.out = 501)
        y <- with(panel, dbeta(x, shape1 = alpha, shape2 = beta))
        lines(x, y, col = "navyblue", lwd = 2)
        media <- with(panel, alpha / (alpha + beta))
        moda <- with(panel, (alpha - 1) / (alpha + beta - 2))
        var <- with(panel, (alpha * beta) / (((alpha + beta)^2) * (alpha + beta + 1)))
        q95 <- with(panel, qbeta(0.95, shape1 = alpha, shape2 = beta))
        text(0.75, 7.3, paste("Moda:", round(moda, 3)), adj = 0)
        text(0.75, 7.0, paste("Q95:", round(q95, 3)), adj = 0)
        text(0.75, 6.7, paste("Media:", round(media, 3)), adj = 0)
        text(0.75, 6.4, paste("Variancia:", round(var, 4)), adj = 0)
        panel
    }
    
    panel <- rp.control(title = "Parametros", alpha = 2, beta = 2)
    rp.slider(panel = panel, variable = alpha, from = 0, to = 50, 
              action = plotbeta, title = "Alpha (α)",
              showvalue = TRUE, resolution = 0.5)
    rp.slider(panel = panel, variable = beta, from = 0, to = 50, 
              action = plotbeta, title = "Beta (β)",
              showvalue = TRUE, resolution = 0.5)
    rp.do(panel, action = plotbeta)
}
rpBeta()

##
##
##
achaBeta <- function(Moda, q95, plot = FALSE){
    f <- function(logalpha, Moda, q95){
        alpha <- exp(logalpha)
        beta <- (alpha-1)/Moda - alpha + 2
        (qbeta(0.95, shape1 = alpha, shape2 = beta) - q95)^2
    }
    options(warn = -1)
    res <- optimize(f, interval = c(-3,6), Moda = Moda, q95 = q95)$min
    alpha <- exp(res)
    beta <- (alpha-1)/Moda - alpha + 2
    if(plot){
        x <- seq(0, 1, length.out = 501)
        y <- dbeta(x, shape1 = alpha, shape2 = beta)
        plot(x, y, type = "l", col = "navyblue", lwd = 2,
             xlab = "Parâmetro", ylab = "Densidade", 
             main = paste("Priori: Beta(", round(alpha, 3), " , ", round(beta, 3), ")", sep="")
             )
        media <- alpha / (alpha + beta)
        moda <- (alpha - 1) / (alpha + beta - 2)
        var <- (alpha * beta) / (((alpha + beta)^2) * (alpha + beta + 1))
        q95 <- qbeta(0.95, shape1 = alpha, shape2 = beta)
        mtext(paste("Moda:", round(moda, 3), " | Q95:", round(q95, 3), " | Media:", round(media, 3), " | Variancia:", round(var, 4)))
    }
    return(c(alpha, beta))
}
achaBeta(Moda = 0.1, q95 = 0.40)
achaBeta(Moda = 0.1, q95 = 0.40, plot = TRUE)
achaBeta(Moda = 0.05, q95 = 0.10)
achaBeta(Moda = 0.5, q95 = 0.80)

rpBeta2 <- function() {    
    plotbeta <- function(panel) {
        ab <- with(panel, achaBeta(Mo = Mo, q95 = q95))
        panel.ab <- list(alpha = ab[1], beta = ab[2])
        with(panel.ab,
             plot(1, type = "n", xlim = c(0, 1), ylim = c(0, 7.5),
                  xlab = "Parâmetro", ylab = "Densidade", 
                  main = paste("Priori: Beta(", round(alpha, 3), " , ", round(beta, 3), ")", sep="")
                  ))
        x <- seq(0, 1, length.out = 501)
        y <- with(panel.ab, dbeta(x, shape1 = alpha, shape2 = beta))
        lines(x, y, col = "navyblue", lwd = 2)
        media <- with(panel.ab, alpha / (alpha + beta))
        moda <- with(panel.ab, (alpha - 1) / (alpha + beta - 2))
        var <- with(panel.ab, (alpha * beta) / (((alpha + beta)^2) * (alpha + beta + 1)))
        q95 <- with(panel.ab, qbeta(0.95, shape1 = alpha, shape2 = beta))
        mtext(paste("Moda:", round(moda, 3), " | Q95:", round(q95, 3), " | Media:", round(media, 3), " | Variancia:", round(var, 4)))
        panel
    }
    
    panel <- rp.control(title = "Opiniao", Mo = 0.5, q95 = 0.80)
    rp.slider(panel = panel, variable = Mo, from = 0, to = 1, 
              action = plotbeta, title = "Moda", showvalue = TRUE,
              resolution = 0.01)
    rp.slider(panel = panel, variable = q95, from = 0, to = 1, 
              action = plotbeta, title = "Quantil 0.95", showvalue = TRUE,
              resolution = 0.01)
    rp.do(panel, action = plotbeta)
}
rpBeta2()