Reputation: 115
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
Upvotes: 0
Views: 255
Reputation: 115
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)
Upvotes: 0