Permalink
Switch branches/tags
Nothing to show
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
118 lines (97 sloc) 3.86 KB
BetaBinoMLE = function (success, trials, start = NULL, optim.method = "default",
lower = 0, upper = Inf) {
#################################################################################
# MLE estimate of Beta-Binomial parameters
#
# Args:
# success: vector of #success; trials:=vector of #trials;
# start: initial parameters (must be a list with name shape1, shape2)
# optim.method: optimization methods in optim(){stats}
# lower(upper): lower(upper) bound for parameters
#
# Returns:
# $estimate: MLE estimate for beta parameters
# $convergence: convergence code from optim(). 0 means good.
# $loglik: Loglikelihood with estimated parameters
# $starting: initial parameters from the method of moments
#
# Dependent package: VGAM
#
# Note: The structure of the function heavily relies on
# mledist(){fitdistrplus} by Marie Laure Delignette-Muller.
#
# Summer2013 @Supplyframe
#################################################################################
if(!is.element("VGAM", installed.packages()[,1])){
stop("Please install and load package VGAM before using this function.")
}
require(VGAM)
distname <- "betabinom.ab"
ddistname <- paste("d", distname, sep = "")
if (is.null(start)) {
if (distname == "betabinom.ab") {
if (any(success/trials < 0) | any(success/trials > 1)) {
stop("Proportion must be in [0-1] to fit a betabinom distribution")
}
start.mu <- mean(success/trials)
start.var <- var(success/trials)
start.a <- ((1 - start.mu) / start.var - 1 / start.mu) * start.mu ^ 2
start.b <- start.a * (1 / start.mu - 1)
start <- list(shape1 = start.a, shape2 = start.b)
}
if (!is.list(start))
stop("'start' must be defined as a named list for this distribution")
}
vstart <- unlist(start)
argddistname <- names(formals(ddistname))
m <- match(names(start), argddistname)
if (any(is.na(m)) || length(m) == 0)
stop("'start' must specify names which are arguments to 'distr'")
fnobj <- function(par, x, n, ddistnam) {
-sum(do.call(ddistnam, c(list(x), list(n), par, log = TRUE)))
}
if (optim.method == "default") {
if (is.infinite(lower) && is.infinite(upper)) {
if (length(vstart) > 1)
meth <- "Nelder-Mead"
else meth <- "BFGS"
}
else meth <- "L-BFGS-B"
}
else meth <- optim.method
opttryerror <- try(opt <- optim(par = start, fn = fnobj,
x = success, n = trials, ddistnam = ddistname,
hessian = TRUE, method = meth, lower = lower,
upper = upper), silent = TRUE)
if (inherits(opttryerror, "try-error")) {
warnings("The function optim encountered an error and stopped")
print(opttryerror)
return(list(estimate = rep(NA, length(vstart)), convergence = 100,
loglik = NA, hessian = NA))
}
if (opt$convergence > 0) {
warnings("The function optim failed to converge, with the error code ",
opt$convergence)
return(list(estimate = rep(NA, length(vstart)), convergence = opt$convergence,
loglik = NA, hessian = NA, message = opt$message))
}
res <- list(estimate = opt$par, convergence = opt$convergence,
loglik = -opt$value, initial = vstart)
return(res)
}
EbPropEstor = function(success, trials){
# Empirical Bayes estimator for binomial proportion
#
# Args: success: a vector of #success
# trials: a vector of #trials
#
# Returns: a vector of estimated proportion.
#
# Dependent Function: BetaBinoMLE()
# Summer2013 @Supplyframe
est = BetaBinoMLE(success, trials)$estimate
a = est[1]
b = est[2]
proportion = (a + success)/(a + b + trials)
return(proportion)
}