[R] Three-component Negative Binomial Mixture: R code

[Achim Zeileis]

On Mon, 7 Nov 2016, danilo.carita at 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 at 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.