Question: Fit Mixture Of Gamma And Negative Binomial Distributions In R
0
FGV110 wrote:

Hi,

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?

FGV

R • 5.2k views
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. http://en.wikipedia.org/wiki/Negative_binomial_distribution#Gamma.E2.80.93Poisson_mixture

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

Yes, the data are NGS read depth

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.

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

2
Arun2.3k wrote:

How about this? you'll need the BB package

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

# generate test data to check model efficiency
set.seed(100)
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
}
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, x, x, x, x)
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?

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