On Mon, 7 Nov 2016, danilo.car...@uniparthenope.it wrote:
I need a function for R software which computes a mixture of Negative Binomial distributions with at least three components.
The package "countreg" on R-Forge provides a driver FLXMRnegbin() that can be combined with the "flexmix" package (i.e., functions flexmix() and stepFlexmix()). The manual page provides some worked illustrations in example("FLXMRnegbin", package = "countreg").
Note that the driver is mainly designed for negative binomial _regression_ models. But if you just regress on a constant (y ~ 1) you can also get negative binomial mixture distributions without covariates.
I found on another site the following function "mixnbinom". It works very well, but it computes a mixture of only two components:mixnbinom=function(y,k1,mu1,k2,mu2,prob,eps=1/100000) { new.parms=c(k1,mu1,k2,mu2,prob) err=1 iter=1 maxiter=100 hist(y,probability=T,nclass=30,col="lightgrey",main="The EM algorithm") xvals=seq(min(y),max(y),1) lines(xvals,prob*dnbinom(xvals,size=k1,mu=mu1)+ (1-prob)*dnbinom(xvals,size=k2,mu=mu2),col="green") while(err>eps){ if(iter<=maxiter){ lines(xvals,prob*dnbinom(xvals,size=k1,mu=mu1)+ (1-prob)*dnbinom(xvals,size=k2,mu=mu2),lty=3) } bayes=(prob*dnbinom(y,size=k1,mu=mu1))/((prob* dnbinom(y,size=k1,mu=mu1))+((1-prob)*dnbinom(y,size=k2,mu=mu2))) mu1=sum(bayes*y)/sum(bayes) mu2=sum((1-bayes)*y)/sum((1-bayes)) var1=sum(bayes*(y-mu1)^2)/sum(bayes) var2=sum((1-bayes)*(y-mu2)^2)/sum((1-bayes)) k1=abs(mu1/((var1/mu1)-1)) k2=abs(mu2/((var2/mu2)-1)) prob=mean(bayes) old.parms=new.parms new.parms=c(k1,mu1,k2,mu2,prob) err=max(abs((old.parms-new.parms)/new.parms)) iter=iter+1 } lines(xvals,prob*dnbinom(xvals,size=k1,mu=mu1)+ (1-prob)*dnbinom(xvals,size=k2,mu=mu2),col="red") print(list(k1=k1,mu1=mu1,k2=k2,mu2=mu2,p=prob,iter=iter,err=err)) }I would be grateful if someone can modify the previous function to model a three-component mixture instead of a two-component one. I also tried to look for a package which does the same job: I have used the package "mixdist", but I am not able to set up a suitable set of starting parameters for the function mix (they always converge to zero). Hereafter, I found the package "DEXUS", but the related function does not provide good estimates for the model, even in the event that I already know what results I have to expect. Any help is highly appreciated. Danilo Carità ------------------------------------------------------------- Danilo Carità PhD Candidate University of Naples "Parthenope" Dipartimento di Studi Aziendali e Quantitativi via G. Parisi, 13, 80132 Napoli - Italy ______________________________________________ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
______________________________________________ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
