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, <[email protected]> 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 > > ______________________________________________ > [email protected] 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. ______________________________________________ [email protected] 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.
