Statistiques bayésiennes: Simuler en R, l'expérience décrite ci-dessus 10 000 fois et noter à chaque fois la longueur de la plus longue course

Question

J'essaie de trouver la moyenne de la plus longue série de 10 000 simulations en retournant 30 fois une pièce. J'ai besoin de simuler en R, l'expérience décrite plus haut 10 000 fois et noter à chaque fois la longueur du plus long parcours.

Voici mon code à ce jour:

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) 

Comment trouver la moyenne la plus longue piste de 10.000 simulations?

Answer 1

Toutes les flips de pièces de monnaie sont indépendants les uns des autres (à savoir, le résultat d'une bascule n'influence pas une autre FLIP). Pour cette raison, nous pouvons retourner toutes les pièces pour toutes les simulations à la fois et ensuite formater de telle sorte qu'il sera plus simple de résumer chaque 30 essai de retournement. Voici comment j'aborderais cela.

# 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) 

Vous devriez obtenir quelque chose entre 18-19%, mais cela va varier un peu à chaque fois que vous essayez cette simulation.

Answer 2

Je pense que ce qui suit fait ce que vous cherchez.

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) 

Pour une explication du code mieux pour inspecter un petit segment de coin (comme coin[20]) et le rle de celui-ci (rle(coin[20])). Les longueurs sont calculées pour chaque segment d'exécution, de sorte que le max(rle(coin)$length) donne l'exécution maximale.

EDIT: Après pourrait être plus rapide

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