R中患者健康的编程概率

Question

首先，我是一个完全的 R 新手。
我有以下问题：

考虑一个简单的进步模型 year-by-year。在第 i 年，令 W_i = 患者健康，I_i = 患者生病，D_i = 患者死亡。转换可以建模为一组条件概率。

令 L = 患者健康的年数。

我得出 L 的概率质量函数为 P(L)=(1-p)(p)^{L-1}。
给定的信息是患者在第 1 年身体状况良好，并考虑到他们的年龄和风险因素，对于所有 i

，P(W_{i+1}|W_{i})=0.2

问题是在 R 中编写一个函数来模拟单个患者的轨迹和 returns 患者健康的年数。

我认为可以使用 rbinom 函数在 R 中将其编程为二项分布。对于单个患者，

rbinom(1, 1, 0.2)

但我认为这不会 return 患者康复的年数。我认为 rbinom 函数应该是开始，它需要搭配一种方法来计算患者健康的年数，但我不知道该怎么做。

下一个问题是使用 R 模拟 1000 个患者轨迹并找到健康年数的样本均值。我假设这将是前一部分的扩展，只需用 1000 名患者替换 1 名患者。但是我不太清楚在哪里用 1000 名替换 1：n 或 size

rbinom(n, size, prob)

这是假设首先使用 rbinom 是正确的做法... 如果我要用另一种编程语言（比如 Python）执行此操作，我将使用一个以 patient_status=W 为条件的 while 循环，并从 L=0 开始迭代循环并在每次成功的迭代中加 1。我不确定 R 是否以同样的方式工作。

Answer 1

让我们从 rbinom(1, 1, 0.2) 的作用开始：它将 returns 1 1 个独立伯努利（即 0-1）个随机变量的实例加在一起有 0.2 等于 1 的概率。因此，该行只会输出 0（80% 的时间执行）或 1（另外 20% 的时间执行） .如您所述，这不是您想要的。

这里的问题是随机变量的选择。二项式变量非常适合 "I roll ten dice. How many land on 6?"，因为它具有以下基本组成部分：

结果分为成功/失败
固定次数（十次）的试验
一致的成功概率 (1/6)
独立试验（骰子互不影响）

您描述的情况不具备这些特征。那么，怎么办？

选项 1：根据您的直觉进行 while() 循环。作为序言，由于各种原因（主要是效率低下），R 中不鼓励使用 while() 循环。但是，既然您已经理解了这个概念，那么让我们运行来理解它。

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 参数更改为更乐观的参数（例如 0.99）来对此进行试验，以模拟长寿命。

选项 2：重新考虑随机变量。我在上面提到变量不是二项式的；事实上，它是几何。像 "I roll a die until it lands on 6. How many rolls did it take?" 这样的情况是几何的，因为它具有以下基本组成部分：

结果分为成功/失败
一致的成功概率
在第一次成功时终止的重复试验
独立试验

就像二项式变量在 R 中如何具有有用的函数，例如 rbinom(), pbinom(), qbinom(), dbinom()，几何变量有一个相应的集合：rgeom(), pgeom(), qgeom(), dgeom().

要使用rgeom()，我们需要注意一个细节：这里，"success"被表征为病人生病，因为那是当实验结束时。（上文中，通过将患者健康编码为 1，我们隐含地使用了反向视角。）这意味着 "success" 概率为 0.8。 rgeom(1, 0.8) 将 return 抽取次数 严格在 第一次成功之前，这相当于患者从一口井到另一口井的年数，如上所述。请注意，参数 1 指的是我们想要运行这个实验的次数，而不是其他的。因此：

rgeom(1, 0.8)

将完成与我们上面定义的 one_patient() 函数相同的任务。（也就是说，每个输出的分布都是相同的。）

对于多个患者，可以将one_patient()函数包裹在replicate()中，也可以直接调整rgeom(1, 0.8)的第一个参数。第二个选项快很多，尽管如果只模拟 1000 名患者，两者都很快。

附录

两者作用相同的证明：

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

R中患者健康的编程概率

Programming probability of patients being well in R

r

probability

附录