0

我写了下面的代码，我需要重复 100 次，我知道我需要使用另一个 for 循环，但我不知道该怎么做。这是代码：

mean <- c(5,5,10,10,5,5,5)
x <- NULL
u <- NULL
delta1 <- NULL
w1 <- NULL

for (i in 1:7 ) {
x[i] <- rexp(1, rate = mean[i])
u[i] <- (1/1.2)*runif(1, min=0, max=1)
y1 <- min(x,u) 


if (y1 == min(x)) {
delta1 <- 1
}
else {
delta1 <- 0
}

if (delta1 == 0)
{
w1 <- NULL
}

else {

if(y1== x[[1]])
{
w1 <- "x1"
} 
}

}
output <- cbind(delta1,w1)
output

我希望最终输出为 100 行 * 3 列矩阵，表示运行编号、delta1 和 w1。

任何想法都会受到真正的赞赏。

4

2 回答 2

2

Here's what I gather you're trying to achieve from your code:

  1. Given two vectors drawn from different distributions (Exponential and Uniform)
  2. Find out which distribution the smallest number comes from
  3. Repeat this 100 times.

Theres a couple of problems with your code if you want to achieve this, so here's a cleaned up example:

rates <- c(5, 5, 10, 10, 5, 5, 5)  # 'mean' is an inbuilt function
# Initialise the output data frame:
output <- data.frame(number=rep(0, 100), delta1=rep(1, 100), w1=rep("x1", 100))
for (i in 1:100) {
    # Generating u doesn't require a for loop. Additionally, can bring in
    # the (1/1.2) out the front.
    u <- runif(7, min=0, max=5/6)
    # Generating x doesn't need a loop either. It's better to use apply functions
    # when you can!
    x <- sapply(rates, function(x) { rexp(1, rate=x) })
    y1 <- min(x, u)

    # Now we can store the output
    output[i, "number"] <- y1
    # Two things here: 
    #   1) use all.equal instead of == to compare floating point numbers
    #   2) We initialised the data frame to assume they always came from x.
    #      So we only need to overwrite it where it comes from u.
    if (isTRUE(all.equal(y1, min(u)))) {
        output[i, "delta1"] <- 0
        output[i, "w1"] <- NA # Can't use NULL in a character vector.
    }
}
output
于 2013-06-05T06:44:29.063 回答
1

这是另一种更有效的方法replicate：

Mean <- c(5, 5, 10, 10, 5, 5, 5)

n <- 100 # number of runs

res <- t(replicate(n, {
           x <- rexp(n = length(Mean), rate = Mean)
           u <- runif(n = length(Mean), min = 0, max = 1/1.2)
           mx <- min(x)
           delta1 <- mx <= min(u)
           w1 <- delta1 & mx == x[1]
           c(delta1, w1)
           }))

output <- data.frame(run = seq.int(n), delta1 = as.integer(res[ , 1]), 
                     w1 = c(NA, "x1")[res[ , 2] + 1])

结果：

head(output)

#     run delta1   w1
#   1   1      1 <NA>
#   2   2      1 <NA>
#   3   3      1 <NA>
#   4   4      1   x1
#   5   5      1 <NA>
#   6   6      0 <NA>
于 2013-06-05T08:03:05.190 回答