[R] How to avoid the NaN errors in dnbinom?

[francogrex]

Hi, The code below is giving me this error message:
Error in while (err > eps) { : missing value where TRUE/FALSE needed
In addition: Warning messages:
1: In dnbinom(x, size, prob, log) : NaNs produced
2: In dnbinom(x, size, prob, log) : NaNs produced

I know from the help files that for dnbinom "Invalid size or prob will
result in return value NaN, with a warning", but I am not able to find a
workaround for this in my code to avoid it. I appreciate the help. Thanks.
Below is the reproducible code:

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=mu1/((var1/mu1)-1)
k2=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))
}

y=c(rnbinom(900,size=10,mu=5),rnbinom(100,size=15,mu=10))

vals=mixnbinom(y,10,7,20,12,0.8)

-- 
View this message in context: http://www.nabble.com/How-to-avoid-the-NaN-errors-in-dnbinom--tf4680211.html#a13373226
Sent from the R help mailing list archive at Nabble.com.