Давайте начнем с того, что делает rbinom(1, 1, 0.2): он возвращает 1 экземпляр 1 независимых бернуллиевских (то есть, 0-1) случайных величин, сложенных вместе, которые имеют вероятность 0.2 быть равной 1. Таким образом, эта строка будет выдавать только 0 (что она будет делать в 80% случаев) или 1 (что она будет делать в остальные 20% времени). Как вы заметили, это не то, что вам нужно.
Проблема здесь заключается в выборе случайной величины. Биноминальная переменная отлично подходит для чего-то вроде: «Я бросаю десять кубиков. Сколько земли на 6?» потому что он имеет следующие существенные компоненты:
- результаты, дихотомические в успех / неудачу
- фиксированное число (десять) испытаний
- постоянная вероятность успеха (1 / 6)
- независимые испытания (игральные кости не влияют друг на друга)
Ситуация, которую вы описываете, не имеет этих особенностей. Итак, что делать?
Вариант 1: Go с вашим инстинктом для while() l oop. Я предскажу это, сказав, что циклы while() не рекомендуются в R по разным причинам (главным образом, из-за неэффективности). Но, поскольку вы уже понимаете эту концепцию, давайте поработаем с ней.
one_patient <- function(){
status <- 1 # 1 = healthy, 0 = ill
years <- (-1) # count how many years completed while healthy
while(status == 1){
years <- years + 1 # this line will run at least one time
status <- rbinom(1, 1, 0.2) # your rbinom(1, 1, 0.2) line makes an appearance!
}
return(years)
}
Теперь выполнение one_patient() приведет к числу лет, в течение которых пациент успешно перешел из скважины в скважину. Это будет как минимум 0, так как years начинается с -1 и увеличивается как минимум один раз. Если пациенту повезет, он может быть очень высоким, хотя, скорее всего, нет. Вы можете поэкспериментировать с этим, изменив параметр 0.2 на что-то более оптимистичное c, например 0,99, для имитации длинных периодов жизни.
Опция 2: переосмыслить случайную переменную. Я упоминал выше, что переменная не была биномиальной; на самом деле это геометрия c. Ситуация типа: «Я бросаю объявление ie, пока оно не приземлится на 6. Сколько роллов потребовалось?» является геометрией c, потому что он имеет следующие существенные компоненты:
- результаты, дихотомические в успех / неудачу
- постоянная вероятность успеха
- повторные испытания, которые прекращаются при первый успех достигнут
- независимых испытаний
Так же, как биномиальные переменные имеют полезные функции в R, такие как rbinom(), pbinom(), qbinom(), dbinom(), существует соответствующая коллекция для геометрических переменных c: rgeom(), pgeom(), qgeom(), dgeom().
Чтобы использовать rgeom(), мы должны быть осторожны с одной деталью: здесь «успех» характеризуется как пациент, заболевший , потому что именно тогда эксперимент концы. (Выше, кодируя пациента как 1, мы неявно используем обратную перспективу.) Это означает, что вероятность «успеха» равна 0,8. rgeom(1, 0.8) вернет число др aws строго до первого успеха, что эквивалентно количеству лет, в течение которых пациент прошел путь от нуля к здоровью, как указано выше. Обратите внимание, что параметр 1 относится к числу раз, когда мы хотим запустить этот эксперимент, а не к чему-то еще. Следовательно:
rgeom(1, 0.8)
выполнит sh ту же задачу, что и функция one_patient(), которую мы определили выше. (То есть распределение выходов для каждого будет одинаковым.)
Для нескольких пациентов вы можете либо обернуть функцию one_patient() внутри replicate(), либо просто настроить первый параметр rgeom(1, 0.8). Второй вариант - намного быстрее, хотя оба быстры, если просто имитировать 1000 пациентов.
Addendum
Доказательство того, что оба имеют одинаковый эффект:
sims1 <- replicate(10000, one_patient())
hist(sims1, breaks = seq(-0.5, max(sims1) + 0.5, by = 1))
sims2 <- rgeom(10000, 0.8)
hist(sims2, breaks = seq(-0.5, max(sims2) + 0.5, by = 1))
Доказательство того, что rgeom() быстрее:
library(microbenchmark)
microbenchmark(
replicate(10000, one_patient()),
rgeom(10000, 0.8)
)
#Unit: milliseconds
# expr min lq mean median uq max neval
# replicate(10000, one_patient()) 35.4520 38.77585 44.135562 43.82195 46.05920 73.5090 100
# rgeom(10000, 0.8) 1.1978 1.22540 1.273766 1.23640 1.27485 1.9734 100