R 仿真编程

Question

赌徒废墟。在这种情况下，赌徒以 6 美元开始。想象一下游戏是抛硬币，它有 1/2 的赢/输概率。现在，每场胜利给你 1 美元，每场失败 -1 美元。下面的代码将多次模拟这种情况，如果他达到 0$ 或某个金额，例如 10$，它将停止。但问题是我不知道如何存储他的股份轨迹，例如在试验 1 中将显示 6 5 4 3 2 1 0。我该怎么做？

gamble <- function(k,n,p) {                                             
   stake <- k                                   
   while (stake > 0 & stake < n) {
         bet <- sample(c(-1,1),1,prob=c(1-p,p))
         stake <- stake + bet }                                                     
    if (stake == 0) return(1) else return(stake)}  
         storage <- vector("list", 100)                                     
         k <- 6       
         n <-  10  
         p <- 1/2  
         trials <- 100
    simlist <- replicate(trials, gamble(k, n, p))              
    print(simlist)

Answer 1

我修改了gamble，这样stake 不是每次都更新一个stake 值，而是一个向量，我们用i 跟踪我们在其中的位置。一种可怕的方法是在每次迭代时将一个新值附加到stake - 一次使向量更长一个项目是非常低效的。相反，我们用一个慷慨的 10k NA 值初始化 stake。如果我们用完了，我们会在最后再粘上 10k。

否则我会尽可能多地保留您的代码。

gamble <- function(k, n, p) {
  stake <- rep(NA_real_, 1e4)
  i <- 1
  stake[1] <- k
  while (stake[i] > 0 & stake[i] < n) {
    bet <- sample(c(-1, 1), 1, prob = c(1 - p, p))
    stake[i + 1] <- stake[i] + bet
    i <- i + 1
    if (length(stake) == i) stake <- c(stake, rep(NA_real_, 1e4))
  }
  return(stake[!is.na(stake)])
}

k <- 6
n <-  10
p <- 1 / 2
trials <- 100
simlist <- replicate(trials, gamble(k, n, p))
head(simlist)
# [[1]]
# [1] 6 5 4 3 4 3 2 1 0
# 
# [[2]]
#  [1]  6  7  6  5  6  5  4  3  2  3  4  3  4  5  4  5  4  5  6  7  8  7  8  7  6  7  8  7
# [29]  8  9  8  7  6  7  6  7  8  9 10
# 
# [[3]]
# [1] 6 5 4 3 2 1 0
# 
# [[4]]
#  [1] 6 7 8 9 8 7 6 5 6 5 4 5 6 7 6 7 6 5 4 3 2 3 2 3 4 3 2 1 2 1 0
# 
# [[5]]
#  [1] 6 5 6 5 4 3 4 3 2 3 4 3 4 3 4 3 4 3 4 5 4 5 6 5 6 7 6 5 4 5 4 5 4 3 4 3 2 1 2 1 2 3
# [43] 2 3 2 3 2 1 0
# 
# [[6]]
#  [1]  6  7  6  7  8  7  6  7  8  9 10

Answer 2

这是 gamble 函数的修改版本：在 while 循环之前初始化的空 track 数组将跟踪赌注的不同值，直至达到最小值或最大值

gamble <- function(s, mi, ma, p){
  stake <- s
  track <- array()
  counter <- 1
  while(stake > mi & stake < ma) {
    bet <- sample(c(-1,1),1,prob=c(1-p,p))
    stake <- stake + bet
    track[counter] <- stake
    counter = counter + 1
    if (counter > 20) break
  }
  return(track)
}

p <- 0.5
starting_value <- 6
mi <- 0
ma <- 10
trials <- 10
#track <- gamble(starting_value, mi, ma, p)
simlist <- replicate(trials, gamble(starting_value, mi, ma, p)) 

end_sims <- vector()
counter <- 1
for (i in 1:trials) {
  if (simlist[[i]][length(simlist[[i]])] == 0 | simlist[[i]][length(simlist[[i]])] == 10) {
    end_sims[counter] <- i
    counter <- counter + 1
  }
}

Answer 3

这是一种不同的方法。第一个想法是一次做很多试验。因此，我们有

gamble0 <-
    function(n_trials, k, n, p)
{
    ## create n_trials simulations
    stakes <- rep(k, n_trials)
    trials <- seq_len(n_trials)

    repeat {
        ## bet on all trials still in play, and update
        bet <- sample(c(1, -1), length(trials), TRUE, prob=c(1-p, p))
        stakes[trials] <- stakes[trials] + bet

        ## only continue to follow those trials that have not terminated
        trials <- trials[(stakes[trials] > 0L) & (stakes[trials] < n)]
        if (length(trials) == 0)
            break
    }
    stakes
}

结果是一个结果向量，计算速度很快，因为我们允许 R 进行 矢量化 计算（例如，调用一次 sample() 以生成 length(trials) 结果，而不是调用它 @ 987654324@次）。

> n <- 100000
> system.time(answer <- gamble0(n, 6, 10, .5))
   user  system elapsed 
  0.336   0.000   0.338 
> table(answer) / n
answer
      0      10 
0.39973 0.60027 

要在每次模拟中累积曲目，请使用list() 来跟踪仍在播放的每个曲目和试验。一旦我们记录了所有轨道的结果，通过创建轨道和试验的单个向量（通过unlist()）并使用split()重新拆分轨道，将迭代列表转换为轨道列表基于轨迹的矢量。

gamble2 <-
    function(n_trials, k, n, p)
{
    ## lists to hold tracks
    tracks <- trials <- list()
    ## initial conditions
    i <- 1L
    stakes <- rep(k, n_trials)
    trial <- seq_len(n_trials)
    repeat {
        ## store current tracks
        tracks[[i]] <- stakes
        trials[[i]] <- trial
        ## still more to do?
        idx <- (stakes > 0L) & (stakes < n)
        if (!any(idx))
            break
        ## update tracks that are still in play
        bet <- sample(c(1, -1), sum(idx), TRUE, c(1 - p, p))
        stakes <- tracks[[i]][idx] + bet
        trial <- trials[[i]][idx]
        ## increment step
        i <- i + 1L
    }
    ## reshape results from list-of-iterations to list-of-tracks
    tracks <- unlist(tracks, use.names = FALSE)
    trials <- unlist(trials, use.names = FALSE)
    tracks <- split(tracks, trials)
    ## report results
    list(iterations = i, tracks = tracks)
}

这是相对较快的，并且可以被操纵来调查属性，例如，

> n_trials <- 100000
> system.time(answer <- gamble2(n_trials, 6, 10, .5))
   user  system elapsed 
  2.172   0.000   2.172 
> tracks0 <- unlist(answer$tracks, use.names=FALSE)
> last <- cumsum(lengths(answer$tracks))
> table(tracks0[last]) / n_trials

      0      10 
0.39794 0.60206 
> hist(lengths(answer$tracks))

(gamble1()，自从被删除后，试图变得过于聪明，使用环境来存储迭代；R 在增长向量和列表方面变得更好，所以这种聪明是不必要的；这是也与@Gregor 的避免增长向量的建议相关——通过索引超过末尾 x[i] 或 x[[i]] 来增长向量现在在 R 中具有合理的性能。