Your problems is that your code is generating NaN and 'err' is set to
NaN and therefore the 'while (err > eps)' fails since NaN is not a
logical variable.
You should learn to use debug or the browser. I put the following
statement in your code:
new.parms=c(k1,mu1,k2,mu2,prob)
err=max(abs((old.parms-new.parms)/new.parms))
if (is.nan(err)) browser() # added statement
iter=iter+1
}
that caught the error and the printed out the values of some of your objects:
Called from: mixnbinom(y, 10, 7, 20, 12, 0.8)
Browse[1]> old.parms
[1] 8.233552e+00 5.225557e+00 -2.487524e+05 1.512247e+01 9.668136e-01
Warning message:
In dnbinom(x, size, prob, log) : NaNs produced
Browse[1]> new.parms
[1] NaN NaN NaN NaN NaN
Browse[1]>
This shows that 'new.parms' at that point was all NaN. You need to
check your logic since is does say that "In dnbinom(x, size, prob,
log) : NaNs produced".
On 10/23/07, francogrex <[EMAIL PROTECTED]> wrote:
>
> 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.
>
> ______________________________________________
> R-help@r-project.org mailing list
> 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.
>
--
Jim Holtman
Cincinnati, OH
+1 513 646 9390
What is the problem you are trying to solve?
______________________________________________
R-help@r-project.org mailing list
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.
