EM algorithm with Beta Distribution in R

Question

I would like to write R code to build the dirichlet mixture model. The loglikelihood I used for the beta distribution is as below:

(,)=(−1)ln¯+(−1)ln(1−)¯+lnΓ(+)−lnΓ()−lnΓ()

and I need help for initialising parameters (alpha, beta and setting extra parameters for the gamma distribution used in the code).

###E step: compute the likelihood of each cell over each component
log_like <- function(theta, X) {
  N = nrow(X)
  alpha <- theta[1]
  beta <-  theta[2]
  log_lik <- N*log(pgamma(alpha + beta) - N*log(pgamma(alpha) - N*log(pgamma(beta)) + (alpha -1)*N*log(X) + (beta - 1) *N* log(1-X)
 
 return(-log_lik)
}

###M-step: 
MLE_estimates <- optim(fn = log_like, 
                       par = c(1,1),
                       lower = c(-Inf, -Inf),
                       upper = c(Inf, Inf),
                       hessian = TRUE,
                       method = "L-BFGS-B",
                       # custom Inputs
                       x = vaf_trim$vaf
                       )

And here's the data set I want to fit. (randomly generated from the negative binomial distribution.

0.25 2 0.23 3 0.22 4 0.21 5 0.21 6 0.21 7 0.21 8 0.21 9 0.20 10 0.20 11 0.20 12 0.20 13 0.19 14 0.19 15 0.19 16 0.19

Answer 1

I think I figured out. This might work. Feel free to comment if you find any errors.

log_like <- function(theta, x) {
  log_lik <- dbeta(x, theta[1], theta[2])
  return(-sum((log(log_lik))))
}

vaf_trim <- as.numeric(unlist(vaf_trim))

# optimise likelihood above over parameters
# given parameters, update weights (normalise per row, column means)
#Find MLE estimates
MLE_estimates <- optim(fn = log_like, #Likelihood function
                       par = c(2,5), #Initialization
                       lower = c(-Inf, -Inf), #Lower bound on parameters
                       upper = c(Inf, Inf), #Upper bound on parameters
                       hessian = TRUE, #Return Hessian for SEs
                       method = "L-BFGS-B",
                       # custom Inputs
                       x = vaf_trim
)

#Examine estimates
MLE_par <- MLE_estimates$par
MLE_SE <- sqrt(diag(solve(MLE_estimates$hessian)))
MLE <- data.table(param = c("alpha", "beta"), estimates = MLE_par, sd = MLE_SE)
knitr::kable(MLE)
MLE_mean <- MLE_par[1] / (MLE_par[1] + MLE_par[2])
print(MLE_mean)