0

我需要生成审查百分比不能为 0 或 1 的模拟数据。这就是我使用 while 循环的原因。问题是如果我将计数增加到 10,000（而不是 5），程序会非常慢。我必须用 400 种不同的场景重复这个，所以它非常慢。我试图找出可以将我的代码逐个矢量化的地方。我怎样才能避免while循环并且仍然能够保持条件？

另一种方法是保留 while 循环并生成符合我的条件的 10,000 个数据集的列表，然后将该函数应用于列表。这里我以汇总函数为例，但我的实际函数同时使用 X_after 和 delta（即 mle(X_after,delta)）。如果我必须使用 while 循环，这是一个更好的选择吗？

我担心的另一个问题是内存问题。在进行如此大的模拟时如何避免耗尽内存？

mu=1 ; sigma=3 ; n=10 ; p=0.10
dset <- function (mu,sigma, n, p) {              
   Mean <- array()
   Median <- array()
   Pct_cens_array <- array()
   count = 0
   while(count < 5) { 

     lod <- quantile(rlnorm(100000, log(mu), log(sigma)), p = p)
     X_before <- rlnorm(n, log(mu), log(sigma))
     X_after <-  ifelse(X_before <= lod, lod,  X_before)
     delta <- ifelse(X_before <= lod, 1,  0) 
     pct_cens <- sum(delta)/length(delta)
     # print(pct_cens)
     if (pct_cens == 0 | pct_cens == 1 ) next
     else {
        count <-  count +1
        if (pct_cens > 0 & pct_cens < 1) {
             sumStats <- summary(X_after)
             Median[count] <- sumStats[3]
             Mean [count]<- sumStats[4]
             Pct_cens_array [count] <- pct_cens 
             print(list(pct_cens=pct_cens,X_after=X_after, delta=delta, Median=Median,Mean=Mean,Pct_cens_array=Pct_cens_array))
          }
       }
    }

          return(data.frame(Pct_cens_array=Pct_cens_array, Mean=Mean, Median=Median)) 
 }
4

2 回答 2

2

我从 C 编程中学到的第一条规则：分而治之！我的意思是你应该首先创建多个函数并将它们调用到你的循环中，因为这个循环做了太多不同的事情。我担心你的算法：

if (pct_cens == 0 | pct_cens == 1 ) next
            else {count <-  count +1

你有什么理由使用while而不是for？while 和 for 循环之间有区别：使用 while，你总是有第一个循环，而不是 for。

最后，关于您的问题：使用更多内存和数组来提高速度。例子：

lod <- quantile(rlnorm(100000, log(mu), log(sigma)), p = p)
            X_before <- rlnorm(n, log(mu), log(sigma))

log(mu) 和 log(sigma) 计算两次：使用变量来存储结果，您会节省时间，但当然会花费更多内存。

于 2012-04-10T06:34:07.590 回答
2

我对您的代码做了一些小调整，但没有改变它的整体风格。听从 Yoong Kim 的建议并尝试将代码分解成更小的部分，以使其更具可读性和可维护性会很好。

  • 您的函数现在有两个“n”参数，用于表示每行中有多少个样本，以及您想要多少次迭代（列）。

  • 您正在增长数组Median并Mean在循环中，这需要大量重新分配内存和复制内容，这会减慢一切。我已经预定义X_after并在循环之后移动了平均值和中值计算以避免这种情况。（作为奖励，mean只median被调用一次而不是n_iteration多次。）

  • ifelse真的不需要打电话。

  • 调用rlnorm一次，为 x 和 lod 生成足够的值，比调用两次要快一些。

这是更新的功能。

dset2 <- function (mu, sigma, n_samples, n_iterations, p) {    
  X_after <- matrix(NA_real_, nrow = n_iterations, ncol = n_samples)
  pct_cens <- numeric(n_iterations)
  count <- 1
  while(count <= n_iterations) {     
    random_values <- rlnorm(2L * n_samples, log(mu), log(sigma))
    lod <- quantile(random_values[1:n_samples], p = p)
    X_before <- random_values[(n_samples + 1L):(2L * n_samples)]
    X_after[count, ] <- pmax(X_before, lod)
    delta <- X_before <= lod
    pct_cens[count] <- mean(delta)
    if (pct_cens > 0 && pct_cens < 1 ) count <- count + 1
  }

  Median <- apply(X_after, 1, median)
  Mean <- rowMeans(X_after)
  data.frame(Pct_cens_array=pct_cens, Mean=Mean, Median=Median) 
}

比较时间，例如，

mu=1
sigma=3
n_samples=10L
n_iterations = 1000L
p=0.10
system.time(dset(mu,sigma, n_samples, n_iterations, p))
system.time(dset2(mu,sigma, n_samples, n_iterations, p))

在我的机器上，有 3 倍的加速。

于 2012-04-10T09:50:03.637 回答