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]]

