贝叶斯统计：模拟在R中，上述实验10,000次，每次注意最长运行的长度

Question

我试图找到10000次模拟中掷硬币30次的最长运行的平均值。我需要在R中模拟上述的实验10,000次，并且每次记录最长运行的长度。

这里是我到目前为止的代码：

coin <- sample(c("H", "T"), 10000, replace = TRUE) 
table(coin) 
head(coin, n = 30) 
rle(c("H", "T", "T", "H", "H", "H", "H", "H", "T", "H")) 
coin.rle <- rle(coin) 
str(coin.rle) 

如何找到最长的平均10,000模拟？

Answer 1

所有的硬币翻转的彼此独立（即一个翻转的结果不会影响另一翻转）的。因此，我们可以一次翻转所有模拟的所有硬币，然后以这种方式进行格式化，这样可以更简单地总结每个30次翻转试验。这是我将如何处理这个问题。

# do all of the flips at once, this is okay because each flip 
# is independent 
coin_flips <- sample(c("heads", "tails"), 30 * 10000, replace = TRUE) 

# put them into a 10000 by 30 matrix, each row 
# indicates one 'simulation' 
coin_matrix <- matrix(coin_flips, ncol = 30, nrow = 10000) 

# we now want to iterate through each row using apply, 
# to do so we need to make a function to apply to each 
# row. This gets us the longest run over a single 
# simulation 
get_long_run <- function(x) { 
    max(rle(x)$length) 
} 

# apply this function to each row 
longest_runs <- apply(coin_matrix, 1, get_long_run) 

# get the number of simulations that had a max run >= 7. Divide this 
# by the number of simulations to get the probability of this occuring. 
sum(longest_runs >= 7)/nrow(coin_matrix) 

您应该得到18-19％之间的东西，但每次尝试这种模拟时都会有所不同。

Answer 2

我认为以下是你所追求的。

n_runs <- 10000 
max_runs <- numeric(n_runs) 
for (j in 1:n_runs) { 
coin <- sample(c("H", "T"), 30, replace = TRUE) 
max_runs[j] <- max(rle(coin)$length) 
} 
mean(max_runs) 

对于代码最好的说明检查的coin一小部分（如coin[20]）和它的rle（rle(coin[20])）。计算每段运行的长度，所以max(rle(coin)$length)给出最大运行。

编辑：下面可能会更快

len <- 30 
times <- 10000 

flips <- sample(c("H", "T"), len * times, replace = TRUE) 
runs <- sapply(split(flips, ceiling(seq_along(flips)/len)), 
        function(x) max(rle(x)$length)) 
mean(runs) # average of max runs 
sum(runs >= 7)/ times # number of runs >= 7