##===================================================
## Proporção de canhotos
## Autor: Gligor
## ==================================================

estima_beta <- function(moda, q95, p = 0.95, kappa_max = 1e5){
    if (moda <= 0 || moda >= 1) stop("moda deve estar entre 0 e 1")
    if (q95 <= moda || q95 >= 1) stop("q95 deve ser maior que a moda e menor que 1")
    
    objetivo <- function(kappa){
        a <- 1 + kappa * moda
        b <- 1 + kappa * (1 - moda)
        valor <- suppressWarnings(qbeta(p, a, b))
    if (is.na(valor)) return(1e10)
        (valor - q95)^2
    }
    
    otimo <- optimize(objetivo, interval = c(0, kappa_max))
    kappa_est <- otimo$minimum
    c(a = 1 + kappa_est * moda, b = 1 + kappa_est * (1 - moda))
}


## ------------------------------------------------------------
## 2) Dados dos alunos
## ------------------------------------------------------------

alunos <- data.frame(
  aluno = c("Arthur", "Axel", "Camila", "Christian", "Cristiane", "João",
            "Julia", "Leonardo E", "Leo Bhering", "Maria", "Matheus", "Pietro",
            "Rodrigo", "Vinicius Santana", "Vinicius Oliveira", "Vitor"),
  moda  = c(20, 8, 10, 8, 20, 20, 18, 5, 11, 10, 1, 3, 10, 2, 6, 7)/100,
  q95   = c(21, 12, 20, 23, 30, 40, 35, 20, 15, 15, 20, 13, 15, 11, 15, 25)/100
)

params <- t(mapply(estima_beta, alunos$moda, alunos$q95))
alunos$a <- round(params[, 1], 1)
alunos$b <- round(params[, 2], 1)

alunos

# adiciona Wagner (priori nao informativa, uniforme) e Gligor (moda alta, 0.80)
alunos <- rbind(
  alunos,
  data.frame(aluno = "Wagner", moda = NA, q95 = NA, a = 1, b = 1),
  data.frame(aluno = "Gligor", moda = 0.80, q95 = 0.95,
             t(round(estima_beta(0.80, 0.95), 1)))
)

print(alunos)


# ------------------------------------------------------------
# 3) Grid de gráficos: priori de cada aluno
# ------------------------------------------------------------
x11()

n_alunos <- nrow(alunos)
nc <- ceiling(sqrt(n_alunos))
nr <- ceiling(n_alunos / nc)

par(mfrow = c(nr, nc), mar = c(3, 3, 1, 0.5), mgp = c(1.8, 0.8, 0))
for (i in 1:n_alunos){
  curve(dbeta(x, alunos$a[i], alunos$b[i]),
        from = 0, to = 1, n = 500,
        xlab = "p", ylab = "densidade",
        main = alunos$aluno[i], col = "steelblue", lwd = 2)
}
par(mfrow = c(1, 1))

### ==========================================================
# Mesma Escala:

n_alunos <- nrow(alunos)
nc <- ceiling(sqrt(n_alunos))
nr <- ceiling(n_alunos / nc)

par(mfrow = c(nr, nc), mar = c(3, 3, 1, 0.5), mgp = c(1.8, 0.8, 0))
for (i in 1:n_alunos){
  curve(dbeta(x, alunos$a[i], alunos$b[i]),
        from = 0, to = 1, n = 500,
        xlab = "p", ylab = "densidade",
        ylim = c(0,20),
        main = alunos$aluno[i], col = "steelblue", lwd = 2)
}
par(mfrow = c(1, 1))

# ------------------------------------------------------------
# 4) Todas as prioris no mesmo gráfico
# ------------------------------------------------------------
cores <- rainbow(n_alunos)
p_seq <- seq(0.001, 0.999, length.out = 500)

dens_mat <- sapply(1:n_alunos, function(i) dbeta(p_seq, alunos$a[i], alunos$b[i]))
ymax <- max(dens_mat)

plot(p_seq, dens_mat[, 1], type = "l", col = cores[1], lwd = 2,
     ylim = c(0, ymax), xlab = "p", ylab = "densidade",
     main = "Prioris de todos os alunos")
for (i in 2:n_alunos) lines(p_seq, dens_mat[, i], col = cores[i], lwd = 2)
legend("topright", legend = alunos$aluno, col = cores, lwd = 2, cex = 0.6)


# ======================================================================
# Escala de 0 a 20:


cores <- rainbow(n_alunos)
p_seq <- seq(0.001, 0.999, length.out = 500)

dens_mat <- sapply(1:n_alunos, function(i) dbeta(p_seq, alunos$a[i], alunos$b[i]))
ymax <- max(dens_mat)

plot(p_seq, dens_mat[, 1], type = "l", col = cores[1], lwd = 2,
     ylim = c(0, 20), xlab = "p", ylab = "densidade",
     main = "Prioris de todos os alunos")
for (i in 2:n_alunos) lines(p_seq, dens_mat[, i], col = cores[i], lwd = 2)
legend("topright", legend = alunos$aluno, col = cores, lwd = 2, cex = 0.6)


# ------------------------------------------------------------
# 5) Verossimilhança + priori + posteriori
# ------------------------------------------------------------
y <- 2
n <- 18

# posteriori = Beta(a + y, b + n - y)
alunos$a_post <- alunos$a + y
alunos$b_post <- alunos$b + n - y

# verossimilhança representada como Beta(y+1, n-y+1) apenas para plotagem na mesma escala
a_verossim <- y + 1
b_verossim <- n - y + 1
d_verossim <- dbeta(p_seq, a_verossim, b_verossim)

selecionados <- c("Arthur", "Julia", "Matheus", "Wagner", "Gligor")

# --- 5a) grid: priori, verossimilhança e posteriori de cada selecionado ---

nc_sel <- ceiling(sqrt(length(selecionados)))
nr_sel <- ceiling(length(selecionados) / nc_sel)

#par(mfrow = c(nr_sel, nc_sel), mar = c(4, 4, 3, 1))
for (nome in selecionados){
  idx <- which(alunos$aluno == nome)
  d_prior <- dbeta(p_seq, alunos$a[idx], alunos$b[idx])
  d_post  <- dbeta(p_seq, alunos$a_post[idx], alunos$b_post[idx])
  ymax_i  <- max(d_prior, d_verossim, d_post)

  x11()
  plot(p_seq, d_prior, type = "l", col = "steelblue", lwd = 2,
       ylim = c(0, ymax_i), xlab = "p", ylab = "densidade",
       main = nome)
  lines(p_seq, d_verossim, col = "darkorange", lwd = 2, lty = 2)
  lines(p_seq, d_post, col = "forestgreen", lwd = 2)
  legend("topright", legend = c("Priori", "Verossimilhança", "Posteriori"),
         col = c("steelblue", "darkorange", "forestgreen"),
         lwd = 2, lty = c(1, 2, 1), cex = 0.7, bty = "n")
}
par(mfrow = c(1, 1))