[R] Gamma mixture models with flexmix

[Wittner, Ben, Ph.D.]

I've been trying with no success to model mixtures of Gamma distributions using
the package flexmix (see examples below). Can anyone help me get it to model
better? Thanks very much.

-Ben

##
##  Please help me get flexmix to correctly model mixtures of
##  Gamma distributions. See examples below.
##
library('flexmix')

##
##  Plot a histogram of dat and the Gamma mixture model given by
##  shapes, rates and pis that is intended as a model of the
##  distribution from which dat was drawn.
##
plotGammaMixture <- function(dat, shapes, rates, pis) {

  KK <- length(pis)
  stopifnot(KK == length(shapes))
  stopifnot(KK == length(rates))

  ho <- hist(dat, plot=FALSE)
  x <- seq(ho$breaks[1], ho$breaks[length(ho$breaks)],
           length.out=1000)
  y <- list()
  for (ii in 1:KK) {
    y[[ii]] <- pis[ii]*dgamma(x, shape=shapes[ii],
                              rate=rates[ii])
  }

  uy <- unlist(y)
  ylim <- if (any(uy == Inf)) {
    c(0, 2*max(ho$intensities))
  } else {
    c(0, max(c(ho$intensities, uy)))
  }
  plot(ho, col='lightgray', freq=FALSE, ylim=ylim,
       main=paste(KK, 'component(s) in model\nmodel prior(s) =',
         paste(round(pis, 2), collapse=', ')))
  cols <- rainbow(KK)
  for (ii in 1:KK) {
    lines(x, y[[ii]], col=cols[ii])
  }
}

##
##  Model dat as a mixture of Gammas then plot.
##
modelGammas <- function(dat, which='BIC') {
  set.seed(939458)
  fmo <- stepFlexmix(dat ~ 1, k=1:3,
                     model=FLXMRglm(family='Gamma'))
  mdl <- getModel(fmo, which=which)
  print(smry <- summary(mdl))
  print(prm <- parameters(mdl))
  plotGammaMixture(dat, prm['shape', ],
                   prm['shape', ]*prm['coef.(Intercept)', ],
                   smry at comptab[, 'prior'])
}

##
##  Works well for a single Gamma distribution.
##
set.seed(78483)
dat1 <- rgamma(6000, shape=2, rate=0.5)
modelGammas(dat1)

set.seed(78483)
dat2 <- rgamma(6000, shape=5, rate=0.3)
modelGammas(dat2)

##
##  Please help me get it to work for mixtures of
##  two Gamma distributions.
##
set.seed(78483)
dat3 <- c(rgamma(6000, shape=3, rate=.5),
          rgamma(4000, shape=5, rate=.1))
modelGammas(dat3)

## Even telling it that there are two components
## does not help. I get two nearly identical distributions,
modelGammas(dat3, which=2)

## whereas what I want to see is something like this.
plotGammaMixture(dat3, shapes=c(3, 5), rates=c(.5, .1),
                 pis=c(.6, .4))

###############################
Here's the output of sessionInfo()
R version 2.12.1 (2010-12-16)
Platform: x86_64-unknown-linux-gnu (64-bit)

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=C              LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     



The information in this e-mail is intended only for the ...{{dropped:11}}