Question: Fit Mixture Of Gamma And Negative Binomial Distributions In R
gravatar for FGV
7.5 years ago by
FGV110 wrote:


I'm trying to fit a mixture of a gamma and negative binomial distribution to a data set in R. In the end the idea is to get the parameters of the two fitted distributions as well as a likelihood; something like this but, as far as I can tell, this package only does mixtures of the same type of distributions (two normal, two poisson, etc..)

I've looked arround and the R package that looks more promising is flexmix but I cannot get it to work. I've been trying with simpler examples (gaussian and poisson) but no good. So far this is what I've got:

ex1 <- stepFlexmix(y~x,
                   data=data, k=2:5,
                   model=list(FLXMRglm(.~., family="gaussian"),
                              FLXMRglm(.~., family="poisson")),
                   control=list(verb=5, iter=100))

Anyone has any idea on how to get it?

Thanks in advance,


R • 5.2k views
ADD COMMENTlink modified 5.4 years ago by muhammad.daniyal0 • written 7.5 years ago by FGV110

(may be irrelevant, but out of curiosity) why do you fit with gamma and negative binomial? Negative binomial could already be obtained by modeling $\lambda$ parameter of poisson distribution as a random variable that is gamma distributed.

ADD REPLYlink written 7.5 years ago by Arun2.3k

I'm analyzing some data that is supposed to be a mixture of two different processes: one has been described as following a negative binomial but the other is a bit more obscure (so I'am trying the gamma). My objective is to try to disentangle both effects.

I know that a neg binomial can be obtained by a poisson and a gamma but am not entirely sure how that can help...

ADD REPLYlink written 7.5 years ago by FGV110

Is this read depth data?

ADD REPLYlink written 7.5 years ago by Zev.Kronenberg11k

Yes, the data are NGS read depth

ADD REPLYlink written 7.5 years ago by FGV110

Yes indeed, you're right. It has nothing to do with your question :). I was thinking more in the line of Zev. Reads from RNA-seq data is usually modeled as a NB distribution as well. I thought you had it gamma and NB confused. Sorry about that.

ADD REPLYlink written 7.5 years ago by Arun2.3k

respected sir


i am working on the mixture of burr XII and weibull dstributions. i have to estimate the parameters through R language. can you kindly guide me how i will do it in R. thanx

ADD REPLYlink written 5.4 years ago by muhammad.daniyal0
gravatar for Arun
7.5 years ago by
Arun2.3k wrote:

How about this? you'll need the BB package

# install.packages("BB") # if you already don't have it.

# generate test data to check model efficiency
gamma.nb.mix <- function(n, prob=0.5, shape=2, scale=15, mu=1, sig=0.1) {
    u <- runif(n) 
    out <- apply( as.matrix(u), 1, function(x) ifelse(x<=0.5, rgamma(1, shape=shape, scale=scale), rnorm(1, mu, sig) ) )
out <- gamma.nb.mix(1000)

# maximizing function
f = function(shape, scale, mu, sig, prob) { 
     -sum(log(prob*dgamma(out, shape=shape, scale=scale) + 
       (1-prob)*dnorm(out, mean=mu, sd=sig))) 

# arbitrary start (probability parameter is important to remain as close as possible)
start0 <- list("shape"=1.2, "scale"=4, "mu"=0.35, "sig"=0.6, "prob"=0.5) 

fcn <- function(x) f(x[1], x[2], x[3], x[4], x[5])   
res <- spg(par=as.numeric(start0), fn=fcn, lower=c(0, 0, -Inf, 0, 0), upper=c(Inf, Inf, Inf, Inf, 1), control=list( maxit=10000 )  )

After 10k iterations:

Initial model parameters: shape=2, scale=15, mu=1, sig=0.1
Final model parameters (non-converged, but after limit reached):  shape=1.92 scale=15.399  mu=0.999  sigma=0.102  prob=0.5149593

Any good?

ADD COMMENTlink modified 7.5 years ago • written 7.5 years ago by Arun2.3k

In general probability/statistics questions are more appropriate on stats.stackexchange

ADD REPLYlink written 7.5 years ago by Arun2.3k

I just realized that I implemented for gamma and normal distribution. I'll try once again with gamma and NB and report back.

ADD REPLYlink written 7.5 years ago by Arun2.3k
Please log in to add an answer.


Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1777 users visited in the last hour