combination.r

data <- list(
  est1 = c(`0.025` = 0.84, `0.5` = 0.9, `0.975` = 0.96),
  est2 = c(`0.025` = 0.32, `0.5` = 0.53, `0.975` = 0.66)
)

fn_gamma <- function(par, x) {
  quantiles <- as.numeric(names(x))
  return(sum((qgamma(quantiles, shape = par[1], rate = par[2]) - x)**2))
}

fit_gamma <- function(quantiles, init) {
  res <- nloptr::sbplx(x0 = init, fn = fn_gamma, x = quantiles,
                       lower = c(shape = 0, rate = 0),
                       control = list(xtol_rel = 1.0e-6, ftol_rel = 1.0e-6))
  sol <- res$par
  names(sol) <- names(init)

  return(sol)
}

fits <- lapply(data, fit_gamma, init = c(shape = 1, rate = 1))

md <- mistr::mixdist(rep("gamma", length(fits)),
                     lapply(seq_along(fits),
                            function(i) c(shape = fits[[i]][["shape"]],
                                          rate = fits[[i]][["rate"]])),
                     weights = rep(1/length(fits), length(fits)))

round(mistr::q(md, c(0.025, 0.5, 0.975)), 2)
## [1] 0.37 0.81 0.95