2016-01-11 4 views
0

Я хочу написать цикл в R, в котором моделируются образцы Пуассона, но я хочу отбросить образцы, которые не содержат нулей и «иметь другой ход». Как я могу это сделать?Отбраковка образцов в соответствии с условием

Например:

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

Есть ли способ, которым я мог бы сохранить образцы, для которых length(X[X==0])==0, а затем повторно образец, и продолжается до 100 средств из образцов, которые действительно содержат нули получаются?

+4

Использовать цикл 'while' с условием' if'? – Frank

+1

Вы имеете в виду 'X [X == 0]' или 'x [x == 0]'? – MichaelChirico

ответ

2

Как @Frank предложил, петля while - ваш лучший подход, хотя я не думаю, что if - лучший способ пойти.

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) 
} 

Затем закончить:

colMeans(draws) 

Альтернативой является использование replicate:

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

Мои быстрые тесты предполагают, что это последнее на самом деле быстрее.

Еще более сообразительным было бы просто устранить все плохие ничьи с самого начала (и в основном извлечь из усеченного распределения).

Мы знаем, что вероятность получения 0 на данном розыгрыше exp(-lambda), поэтому если мы инвертировать равномерная рисует на (exp(-lambda), 1], мы будем устанавливать:

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

Также конкурентоспособны с этим использовать data.table:

library(data.table) 
grps <- rep(1:NN, each = kk) 
data.table(qpois(runif(kk * NN, min = exp(-lam)), lam))[ , mean(V1), grps] 
+0

Спасибо, я действительно искал, чтобы сохранить образцы, которые содержат нули. Я на самом деле не получаю средства, но это метод получения образцов, которые меня интересуют. Как мне изменить вышеизложенное? – Pauljw11

0

Просто чтобы сказать, что я понял, что если я изменить код Мичил на:

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

Он будет делать то, что я желаю. Спасибо всем, кто ответил.

Смежные вопросы