#Utility functions used in sis.seal, written in Splus 6.1 rbinom.0 <- function(n,N,p) { if(length(p) != n) p <- rep(p,n) p[(N==0)|is.na(N)|(p<0)] <- 0 p[p>1] <- 1 N[N==0|is.na(N)] <- 1 out <- rbinom(n,N,p) return(out) } rnbinom.0 <- function(n,k,p) { if(length(p) !=n) p <- rep(p,n) i<-((k==0)|(p==0)|(is.na(k))|(is.na(p))) k[i]<-1 p[i]<-0.5 out <- rnbinom(n,k,p) out[i]<-0 return(out) } dpois.0 <- function(y,mu) { temp <- mu temp[mu==0] <- 1 out <- dpois(y,temp) out[mu==0] <- 0 return(out) } rmultinom <- function(nn, n, pvec) { #generates random variables from a multinomial distribution #nn is the number of multinomial random variates to output # (number of particles in our application) #n is a vector of number to distribute amoung bins. If vector length is 1, and nn is # greater than 1, it assumes you have the same n to distribute for each particle. Then # it is fast. Otherwise, length of n must be nn. #pvec is a matrix containing the ps. #rows=nn; #columns=#bins. if(length(n) == 1) { k <- length(pvec) temp <- matrix(sample(1:k, n * nn, replace = T, prob = pvec),nn, n) out <- apply(temp, 1, tabulate, k) } else { k <- dim(pvec)[[2]] out <- matrix(0,nn,k) if(k <= 5) { out[, 1] <- rbinom.0(nn, n, pvec[, 1]) temp.p <- pvec[, 1] temp.y <- out[, 1] for(i in 2:(k - 1)) { #relies on rbinom.0 being able to cope with p's <0 and >1 out[, i] <- rbinom.0(nn, n - temp.y, pvec[,i]/(1-temp.p)) temp.p <- temp.p + pvec[, i] temp.y <- temp.y + out[, i] } out[, k] <- n - temp.y } else { Q <- matrix(cumsum(t(pvec)) - rep(0:(nn - 1), rep(k, nn)), ncol = k, byrow = T) R <- 1 + (Q[rep(1:nn, n), ] < runif(sum(n))) %*% rep(1, k) temp <- split(R, rep(1:nn, n)) temp <- lapply(temp, tabulate, k) out <- matrix(unlist(temp), nn, k, byrow = T) } } return(out) }