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