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

[danilo.carita at uniparthenope.it]

I need a function for R software which computes a mixture of Negative
Binomial distributions with at least three components.

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