1

感谢任何帮助提高这个循环代码的效率。对不起！

我正在模拟“可变寿命调整显示”。更多细节在这里。

当应用于手术时，这是一个图表，作为预期死亡人数减去实际死亡人数 (e) 的连续统计。假设外科手术在 p=0.1 的人群中具有预期的死亡概率。一个手术单元的性能可以被认为是 0s 和 1s 的向量 (r_(1:n))，其中 0 表示成功，1 表示死亡。在开始 (r_0) 时，e_0 = 0。如果患者存活，则将 p 添加到运行分数中。如果患者死亡，则从运行分数中减去 1-p（例如，其中 r_1=0，e_1 = e_0 + p；其中 r_1=1，e_1 = e_0 - (1-p)。

可以（可能！）看到，如果给定的一组结果的总体死亡概率等于总体的死亡概率，则运行分数 (e) 应该在 0 附近波动。如果死亡率高于预期（更多1s)，e的趋势为负。如果结果好于预期（更多的 0），则趋势为正。

# Simulate VLAD
# m=number of simulations, n=number of procedures, p1=expected mortality
#    p2=actual mortality

vlad_sim <- function(m,n,p1,p2){
  e<-matrix(nrow=n, ncol=m)
  e[1,]<-0
  r<-vector()
  for (j in 1:m){
    r<-rbinom(n,1, p2)
    for (i in 2:n){
      e[i, j] <- ifelse(r[i]==0, e[i-1,j] + p1, e[i-1,j] - (1-p1))
    }
  }
return(e)
}

# Test example using m=100, n=100, p1=0.1, p2=0.2
e <- vlad_sim(100, 100, 0.1, 0.2)

这段代码可以工作并且可以做我想要的。我可以用 ggplot2 制作可爱的情节。我想更改这两个 for 循环以应用函数，但不知道如何。首先，制作尺寸为 nxm 的结果矩阵可能更容易：

r<-matrix(rep(rbinom(n,1, p2),m), nrow=n, ncol=m)

然后我如何将我的函数应用于这个矩阵？谢谢！

4

1 回答 1

3

以下解决方案按行应用计算。您仍然需要一个 for 循环，因为一个阶段的计算是基于前一阶段的结果。此外，现在的代码效率更高。

vlad_sim <- function(m, n, p1, p2){
  e <- matrix(0, nrow = n, ncol = m)
  r <- matrix(sample(c(FALSE, TRUE), size = (n - 1) * m, 
                     replace = TRUE, prob = c(1 - p2, p2)), ncol = m)
  for (i in seq(2L, n)) {
    e[i, ] <- e[i - 1, ] + p1 - r[i - 1, ]
  }   
  return(e)
}
于 2013-07-27T15:59:46.170 回答