A function to estimate the beta model (w/ standard errors)

beta-model.R

# log-likelihood for beta
ll.fn.beta <- function(theta, y) {
  alpha <- theta[1]  # optim() requires a single parameter vector
  beta <- theta[2]
  ll <- alpha*sum(log(y)) + beta*sum(log(1 - y)) - 
    length(y)*log(beta(alpha, beta))
  return(ll)
}

# function to estimate beta model
est.beta <- function(y) {
  # optimize log-likelihood
  mle <- optim(par = c(1, 1), fn = ll.fn.beta, y = y,
               control = list(fnscale = -1),
               hessian = TRUE,
               method = "Nelder-Mead") # for >1d problems
  # check convergence
  if (mle$convergence != 0) print("Model did not converge!")
  # pull out and calculate important quantities
  est <- mle$par
  V <- solve(-mle$hessian)
  se <- sqrt(diag(V))
  # put important quantities in list to return
  res <- list(est = est,
              V = V,
              se = se)
  return(res)
}

# set alpha, beta, and n
alpha <- 2
beta <- 2
n <- 1000

# simulate the data
set.seed(1234)
y <- rbeta(n, alpha, beta)

# estimate shape parameters
m1 <- est.beta(y)
m1