Edward Wijaya
2008-May-23 06:40 UTC
[R] Est. Component Size with AIC/BIC under Gamma Distribution
Dear all, I am trying to model number of samples from a given series. The series are modelled according Gamma function. In order to estimate the # samples, I use BIC/AIC with MLE (computed from dgamma function). Here is the code I have. __BEGIN__ mlogl <- function( x_func, theta_func, samp) { # computing log_likelihood return( - sum(dgamma(samp, shape = x_func, scale=theta_func, log = TRUE))) } find_bic <- function(mll,smpl,k) { bic <- (-2 * mll) + (k * log(length(smpl))) bic } find_aic <- function(mll,smpl,k) { aic <- (-2 * mll) + (k * 2) aic } mlogl_process <- function(smpl,error,start ) { # EM algorithm to estimate max loglikelihood thetalast <- 0 thetacurrent <- start theta <- start difference <- start thetaprocess <- start best <- 0 while (difference > error) { mlogl_out <- nlm(mlogl, mean(smpl), theta_func=thetacurrent, samp=smpl) theta <- mlogl_out$estimate thetalast <- thetacurrent thetacurrent <- theta difference <- abs(thetalast - thetacurrent) # E-STEP new_maxlogl <- nlm(mlogl, mean(smpl), theta_func=theta, samp=smpl) # M-STEP if (new_maxlogl$minimum > best) { best <- new_maxlogl$minimum } } best } # main program # my samples vsamples<- c(14.7, 18.8, 14, 15.9, 9.7, 12.8) # initialize start_ <- 300 error_ <- 0.001 thetastart <- 1 k <- 5 # compute AIC/BIC for k component for (nofc in 1:k { cat("k = ", nofc, "\n") maxlogl <- -(mlogl_process(vsamples,error_,thetastart)) bic_k <- find_bic(maxlogl,vsamples,nofc) aic_k <- find_aic(maxlogl,vsamples,nofc) cat("BIC = ", bic_k, " - AIC = ", aic_k, "- MLL= ", maxlogl,"\n") } __ END__ We finally expect to choose K with the smallest AIC/BIC value. However, the problem I have is that the AIC/BIC value is always on the increase. Can anybody advice what's wrong with my above approach? Regards, Edward