Saving the results of a simulation in a matrix

[swaheera]

• Imagine there is a coin such that if the current flip is heads, the next flip is heads with p=0.7 and tails with p=0.3 If the current flip is tails, the next flip is tails with p=0.7 and heads with p=0.3
• The coin always starts with heads
• We simulate 100 random flips of this coin. Then repeat this 1000 times to create the analysis file
• p1 = probability of flipping heads n steps after your kth heads
• p2 = probability of staring at heads and flipping heads on your nth flip
• From the analysis file, for different values of (n,k) we calculate the absolute value of p1-p2

Here is the code I wrote:

num_flips <- 1000
num_reps <- 10000

coin_flips <- matrix(nrow=num_reps, ncol=num_flips)

for (j in 1:num_reps) {
# Set first flip to be heads
coin_flips[j, 1] <- "H"


for (i in 2:num_flips) {
# If the last flip was heads
if (coin_flips[j, i - 1] == "H") {
# The next flip is heads with probability 0.7
coin_flips[j, i] <- ifelse(runif(1) < 0.7, "H", "T")
} else {
# If the last flip was tails
# The next flip is tails with probability 0.7
coin_flips[j, i] <- ifelse(runif(1) < 0.7, "T", "H")
}
}
}


results <- matrix(nrow=10, ncol=10)

# Loop over k and n
for (k in 1:10) {
for (n in 1:10) {
  
outcomes_P1 <- character(num_reps)
outcomes_P2 <- character(num_reps)
    
# Loop
for (j in 1:num_reps) {
# Find the kth head
indices_of_kth_heads <- which(coin_flips[j, ] == "H")[k]
        
       
outcomes_P1[j] <- coin_flips[j, indices_of_kth_heads + n]
        
        
outcomes_P2[j] <- coin_flips[j, n]
}
    

P1 <- sum(outcomes_P1 == "H") / length(outcomes_P1)
P2 <- sum(outcomes_P2 == "H") / length(outcomes_P2)
    
# Absolute difference between P1 and P2
    results[k, n] <- abs(P1 - P2)
}
}

The results look like this:

print(results)

[,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10]
[1,] 0.316 0.115 0.038 0.001 0.011 0.009 0.007 0.014 0.006 0.002
[2,] 0.306 0.091 0.010 0.009 0.006 0.006 0.027 0.004 0.004 0.008
[3,] 0.285 0.100 0.018 0.006 0.004 0.025 0.004 0.002 0.031 0.049
[4,] 0.296 0.086 0.022 0.006 0.007 0.005 0.005 0.010 0.054 0.063
[5,] 0.291 0.105 0.041 0.004 0.016 0.028 0.008 0.038 0.072 0.056
[6,] 0.297 0.085 0.010 0.026 0.017 0.012 0.030 0.023 0.050 0.044
[7,] 0.274 0.069 0.007 0.008 0.032 0.038 0.019 0.049 0.056 0.060
[8,] 0.282 0.066 0.021 0.030 0.050 0.019 0.029 0.031 0.061 0.043
[9,] 0.284 0.103 0.062 0.040 0.025 0.027 0.027 0.031 0.049 0.054
[10,] 0.309 0.126 0.055 0.007 0.036 0.008 0.027 0.024 0.050 0.050

I think the code is not working as intended. The 0.316 in the top left corner should be zero (assuming that is n=k=1).

How can I fix this code?