I am trying to simulate a walk around a square such that the probability of walking to a vertex adjacent in the square is p/2 going left and p/2 going right, then $1-p$ going diagonally. I've written some code to simulate this and made a function to calculate the number of occurrences of a subset of vertices. When using the starting vertex as the subset, it's getting close to the value of 0.25 that theory tells me I should expect.
move_func <- function(x, n){
pl <- c(x)
for(i in 1:n){
k <- cbind(replicate(p*1000,c(1,0)),replicate(p*1000, c(0,1)),replicate((1-
p)*2000, c(1,1))) #produces matrix of vectors in proportion to probabilities
x <- (x + k[,sample(ncol(k), 1)])%%2 # samples from matrix to get 2 x 1
#vector representing a movement around the square
pl <- cbind(pl, x) #adds new movement to matrix of verticies visited
}
return(pl) # returns final matrix will all verticies visited on the walk
}
test_in_H <- function(x, H){
sum(apply(H,2,identical, x = x)) # tests if vector x is in subset H (a
# matrix)
}
pl <- move_func(x, n)
print(mean(replicate(100, sum(apply(pl, 2, test_in_H, H = H ))/n))) #repeats
#the simulation 100 times and finds the average
}
However, it varies a lot around this even when a lot of steps are taken in the walk and it's repeated many times.
Does anyone know if I'm making a mistake in my code, or if it's simply going to vary quite a bit due to it being a simulation.