[R] Three-component Negative Binomial Mixture: R code
Bert Gunter
bgunter.4567 at gmail.com
Mon Nov 7 19:04:12 CET 2016
Opinion only (and therefore ignorable):
Unless you have a lot (?) of data, fitting such a model successfully
is likely to be very challenging. I suggest that you consult a local
statistical expert to see if you might take a different approach.
On Mon, Nov 7, 2016 at 4:42 AM, <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.
>
> 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.
Isn't that your job?
-- Bert
>
> 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.
More information about the R-help
mailing list