R help - Feb 2011 - Gamma mixture models with flexmix

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}}