0

I wish to write a loop in R within which Poisson samples are simulated, but I wish to discard samples that do not contain any zeros and "have another go". How may I do this?.

For example:

X<-rep(999,100)
for(j in 1:100){
x<-rpois(100,4)
X[j]<-mean(x)
}

Is there any way I could keep samples for which length(X[X==0])==0, and then reselect a sample, and continue until 100 means from samples which do contain zeros are obtained?

MichaelChirico
  • 33,841
  • 14
  • 113
  • 198
Pauljw11
  • 1
  • 1

2 Answers2

2

As @Frank suggested, a while loop is your best approach, though I don't think if is the best way to go.

NN <- 100
kk <- 100
lam <- 4

draws <- matrix(rpois(kk * NN, lam), ncol = NN)

while (!all(idx <- apply(draws, 2, all))){
  draws[ , nidx] <- matrix(rpois(sum(nidx <- !idx) * NN, lam), ncol = NN)
}

Then to finish:

colMeans(draws)

An alternative is to use replicate:

colMeans(replicate(NN, {draws <- rpois(kk, lam)
while (!all(draws)) draws <- rpois(kk, lam)
draws}))

My quick benchmarks suggest this latter is actually faster.

Even more savvy would be to simply eliminate all bad draws from the start (and essentially draw from the truncated distribution).

We know that the probability of getting 0 on a given draw is exp(-lambda), so if we invert uniform draws on (exp(-lambda), 1], we'll be set:

colMeans(matrix(qpois(runif(kk * NN, min = exp(-lam)), lam), ncol = NN))

Also competitive with this is to use data.table:

library(data.table)
grps <- rep(1:NN, each = kk)
data.table(qpois(runif(kk * NN, min = exp(-lam)), lam))[ , mean(V1), grps]
MichaelChirico
  • 33,841
  • 14
  • 113
  • 198
  • Thanks, I was actually looking to retain the samples that DO contain zeroes. I am not actually getting the means, but it is the method of obtaining the samples I am interested in. How do I modify the above? – Pauljw11 Jan 12 '16 at 17:15
0

Just to say that I have realised that if I edit Micheal's code to:

replicate(NN, {draws <- rpois(kk, lam)
while (all(draws)) draws <- rpois(kk, lam)
draws})

It will do what I wish. Thanks to all who answered.

Pauljw11
  • 1
  • 1