R help - Aug 2011 - Possible Error in generic function rzigp in ZIGP Package

Dear All, 
I would like to report a possible mistake in the generic function  rzigp in ZIGP
package, and seek for help.

Yesterday I tried the "rzigp" function, in this way:
1. generate poisson means from a pareto distribuion,
2. generate dispersion for each poisson from gamma distribution
3. generate ratios between poisson means from uniform
4. generate dispersed poisson counts by rzigp.

However, it went fine, when the number (" nmiu" in the codes I ran) of
poisson means was 100, but  when it is 500, or any number larger than 500, I got
this error :

Error in while (s[i] < upper) { : missing value where TRUE/FALSE needed 
In addition: Warning messages: 
1: In log((1 - omega) * mu * (mu + (phi - 1) * (i - 1))^(i - 2)/exp(lgamma(i - 
:
  NaNs produced 
2: In log((mu + (i - 1) * (phi - 1))/(phi * i) * (1 + (phi - 1)/(mu +  : 
  NaNs produced 


I then took a look at the defintion of rgzip at
http://rss.acs.unt.edu/Rdoc/library/ZIGP/html/rzigp.html, and this message seems
to say that the upper bound which is max(runif(n,0,1)) is in comparable to a
NsNs, possibly produced by the two functions displayed in the error message.

Any help on this? Since I will have to use rzigp to do large scale simulation. 


Hereunder are the codes I ran: 

rm(list=ls(all=TRUE)) 
set.seed(123); 
library(MCMCpack) 

# You can use rzigp() in the ZIGP package. rzigp(n,mu,phi,omega=0) for
generalized poisson.

library(ZIGP) # for generalized poisson 
library(VGAM) # for pareto ditribution 


# parameters for pareto to generate mius 
nmiu <- 500;   miupareto <- 3  ; 
disppareto <- 7; 
gmshape <- 1; 
gminvscale <- 0.5; 
gmscale <- 1/gminvscale; 

#generate means for control from Pareto distribution 
miuveccontrol <-rpareto(nmiu,miupareto,disppareto); 

# generate disperssions with gamma distribution 
dispvectpoiss <- rgamma(nmiu, gmshape, gminvscale, gmscale) 

foldchange <- runif(nmiu,min=0.5,max=3); 

miuvectreat <- foldchange*miuveccontrol; 

countvectreat <- matrix(0,nrow=nmiu,ncol=1); 
countveccontrol <- matrix(0,nrow=nmiu,ncol=1); 

  #### rzigp function gives error 
    lgh3 <- length(miuvectreat); 
    for (i in 1:lgh3) 
    { 
    countvectreat[i] <- rzigp(1,miuvectreat[i],dispvectpoiss[i],0); 
    countveccontrol[i] <- rzigp(1,miuveccontrol[i],dispvectpoiss[i],0); 
    }

	[[alternative HTML version deleted]]

Chee Chen <chee.chen <at> yahoo.com> writes:
> 
> Dear All, 
> I would like to report a possible mistake in the generic function
>  rzigp in ZIGP package, and seek for help.
  Unless you get very lucky and somewhere here volunteers to dig into
the package and find out what's wrong, your best bet is to contact
the package maintainer [maintainer("ZIGP")].  However, since the
package
is no longer maintained on CRAN (I had to Google for it and dig into the
archives --  the last version is from 18 months ago), you may be out
of luck ...

  Ben Bolker