Intervalos de credibilidad para la beta: una alternativa

Estadística, R

A partir de los comentarios de Olivier Núñez a mi entrada anterior casi homónima, se nos ha ocurrido a ambos de forma independiente y simultánea una manera alternativa de calcular el intervalo: minimizando su longitud.

a <- 3
b <- 5
alfa <- 0.05

# versión de la entrada anterior:
f <- function(x){
  (dbeta(x[2], a, b) - dbeta(x[1], a, b))^2 +
    (pbeta(x[2], a, b) - pbeta(x[1], a, b) -1 +  alfa)^2
}

res <- optim(c(a/(a+b), a/(a+b)), f)
res$par
#[1] 0.08052535 0.68463436

# nueva versión
f.alt <- function(x){
  qbeta(x+0.95, a, b) - qbeta(x, a, b)
}

res.alt <- optim(0.025, f.alt)
qbeta(c(res.alt$par, res.alt$par + 0.95), a, b)
#[1] 0.08054388 0.68464900