The simulation problem

La tortue est posée au point \((0,0)\).

A l’étape 1, la tortue se déplace de \(+u\) avec \(u\) tiré au hasard uniformément parmi \((0,1),(1,0),(0,-1),(-1,0)\)
A l’étape elle se redéplace par le même procédé aléatoire.
On répète ce déplacement \(n\) fois.
On note \(N_{n}\) le nombre de fois où la tortue revient à un point déjà visité dans le passé.
Pour \(n=100\), puis 1000 puis 10000 , utilisez la méthode de monté Carlo pour avoir une idée graphique de la loi de \(N_{n}\).

\(\textbf{Position and movement}\)

\(N_n\) is the number of times the turtle revisits a point it has already visited in the first \(n\) steps of its random walk. The turtle starts at a point \((x_0, y_0)\) on a two-dimensional grid and it moves to a neighbourhood point chosen uniformly at random from the set \(\{(x_0+1, y0), (x_0-1, y0), (x_0, y_0+1), (x_0, y_0-1)\}\) at each step. Let \((x_i, y_i)\) be the position of the turtle after i steps. Then, Nn is defined as:

\[N_n = \#\{i, j | 0 \leq i < j \leq {n-1} \& (x_i, y_i) = (x_j, y_j)\}\]

Therefore, \(N_n\) is a random variable that depends on the random walk of the turtle and its distribution gives us information about the probability of the turtle revisiting points in the first \(n\) steps of its walk.

We will present the random walk algorithm in 1 dimension. Here the turtle is at position \(0\). It can move to \({-1,1}\) after 1 step and move to either \(-2\) or \(0\) if the first step was \(-1\).

1-D random walk algorithm

no_steps <- 1000
number_Walks <- 300

# We create a matrix to store the position of each step
positions <- matrix(0, ncol = no_steps + 1, nrow = number_Walks)

# We create a loop to calculate the position of each step

for (r in 1:number_Walks)
{
u <- 0 # The initial position at u= 0

for (i in 1:no_steps) 
{
step <- runif(1, -1, 1) # Generate a uniform random number between -1 and 1
u <- u + step
positions[r, i + 1] <- u # The new position of 1 random walk
}
}

ran_walk_1D2 <-function(number_Walks,no_steps)
{
# We create a matrix of uniform generated variables U([-1,1])
mat_grid <- matrix(runif(number_Walks * no_steps, -1, 1), ncol = no_steps)

# We calculate the position of each steps here
positions2 <- apply(cbind(rep(0, number_Walks), mat_grid), MARGIN = 2, FUN = cumsum)

return(positions2)
}

simul1d_1 <- ran_walk_1D2(number_Walks = 1, no_steps = 100)
simul1d_100 <- ran_walk_1D2(number_Walks = 100, no_steps = 100)
simul1d_1000 <- ran_walk_1D2(number_Walks = 100, no_steps = 1000)
simul1d_10000 <- ran_walk_1D2(number_Walks = 100, no_steps = 10000)

Plot of the 1-D random walk

print("A random walk with 100 simulations")

## [1] "A random walk with 100 simulations"

plot(simul1d_1, type = "l", xlab = "Steps", ylab = "Position", main = "Random Walk Path", col = "dark green", fg = "dark green")

Multiple 1-D random walk

matplot(simul1d_100, type = "l", col = 1:number_Walks, xlab = "Number of steps", ylab = "Position")

matplot(simul1d_1000, type = "l", col = 1:number_Walks, xlab = "Number of steps", ylab = "Position")

matplot(simul1d_10000, type = "l", col = 1:number_Walks, xlab = "Number of steps", ylab = "Position")

Probability distribution generated from multiple random walks

hist_100 <- hist(simul1d_100, breaks = 80, freq = FALSE,col = "dark green", fg = "dark green",
                  xlab = "Number of steps", ylab = "Length of random walk",
                   main = "Distribution of the random variables generated")

hist_1000 <- hist(simul1d_1000, breaks = 80, freq = FALSE,col = "dark green", fg = "dark green",
                  xlab = "Number of steps", ylab = "Length of random walk",
                   main = "Distribution of the random variables generated")

hist_10000 <- hist(simul1d_10000, breaks = 80, freq = FALSE,col = "dark green", fg = "dark green",
                   xlab = "Number of steps", ylab = "Length of random walk",
                   main = "Distribution of the random variables generated")

cdf_100 <- plot(ecdf(ran_walk_1D2(100,100)), main = "Empirical cumulative distribution with 100 steps ")

cdf_1000 <- plot(ecdf(ran_walk_1D2(100,1000)), main = "Empirical cumulative distribution with 1000 steps")

cdf_10000 <- plot(ecdf(ran_walk_1D2(100,10000)), main = "Empirical cumulative distribution 10000 steps")

From the graphical results show the histogram and the cumulative density function, we can see that the random walk converges closer to the normal distribution with a higher number of simulations \(M\). Now, we will use monte carlo to have an idea of the distribution of the positions revisited from the turtle. We have considered the chi-squared distribution and poisson distribution for comparison of the distribution that will be generated.

2-D random walk algorithm

Random_Walk_2D <- function(n_sim) 
  
{
pos_desp <- c(-1, 0, 1) # possible displacement values

n_dup <-numeric(n_sim) # initialize a vector to store the number of duplicated positions in each simulation

for (i in 1:n_sim) 
  {
 Alea <- sample(pos_desp, size = 2*n_sim, replace = TRUE) # generate a random sample of displacement values with replacement
 desp <- matrix(Alea, ncol = 2)                           # create a 2D matrix for the displacements with two columns, one for each                                                             dimension
 pos <- apply(X = desp, MARGIN = 2, FUN = cumsum)         # apply the cumsum function to the columns to get step N+1
  
return(pos)

}
}

set.seed(1234)
disp_2D_mc <- function(n_sim) 
  
{
  
pos_desp <- c(-1, 0, 1) # possible displacement values

n_dup <-numeric(n_sim) # initialize a vector to store the number of duplicated positions in each simulation

for (i in 1:n_sim) 
  {
 Alea <- sample(pos_desp, size = 2*n_sim, replace = TRUE) # generate a random sample of displacement values with replacement
 desp <- matrix(Alea, ncol = 2)                           # creates a 2D matrix for the displacements with two columns, one for each                                                             dimension
 pos <- apply(X = desp, MARGIN = 2, FUN = cumsum)         # apply the cumsum function to the columns to get step N+1
  
 n_dup[i] <- sum(duplicated(pos))                         # Counts the sum of duplicated positions
}
prob_dup <- sum(n_dup>0)/n_sim                            # calculate the probability of duplicated positions (N_n/n)

Num_dup <-  max(n_dup)
hist_dup <- hist(n_dup, main = paste("Histogram of duplicated points (Simulations =", n_sim, ")"), 
                                  xlab = "Number of duplicated points", ylab = "Frequency", breaks =n_sim/10,xlim = c(0,Num_dup),col = "dark green", fg = "dark green", )

E_X <- mean(n_dup)
x <- seq(0, max(n_dup))
Poiss <- dpois(x, E_X) #Density of poisson distribution
Chi2 <- dchisq(x, E_X) #Density of chi-square distribution

lines(x,Poiss*sum(n_dup), col = "red") 
lines(x,Chi2*sum(n_dup), col = "blue")
legend("bottomleft", legend = c("Simulated distribution", "Poisson distribution", "Chi Square distribution"), 
       lty = c(1, 1), col = c("black", "red", "blue"), cex = 0.7)

return(n_dup)
print("Histogram of the Simulated duplicated positions\n", hist_dup,
      "\n Number of duplicated positions of the turtle\n", Num_dup,
      ) 
      
}

Probability distribution generated from 100 Simulations having duplicated positions

#simul2d_100 <-disp_2D_mc(100)
simul2d_100.cdf <-plot(ecdf(disp_2D_mc(100)), main = "Empirical cumulative distribution 100 steps",col = "dark green", fg = "dark green")

Probability distribution generated from 1000 Simulations having duplicated positions

#simul2d_1000 <-disp_2D_mc(1000)
simul2d_1000.cdf <-plot(ecdf(disp_2D_mc(1000)), main = "Empirical cumulative distribution 1000 steps",col = "dark green", fg = "dark green")

Probability distribution generated from 10000 Simulations having duplicated positions

#simul2d_10000 <-disp_2D_mc(10000)
simul2d_10000.cdf <-plot(ecdf(disp_2D_mc(10000)), main = "Empirical cumulative distribution 10000 steps" )

We have successfully simulated the random walk in 2 dimension. From the same analysis in the 1 dimension results, we can see that the distribution of the positions duplicated from the Monte-Carlo simulation converges this time to a poisson distribution with the parameter \(\lambda\) for \(P_{o}(\lambda)\) as \(n -> +\inf\) with the mean of the simulation. The result makes sense with the poisson distribution property as it represents occurrences for non regular events.

Project: Simulation and R

A.V Hurtado Vanessa and Narain Ritish

28/03/2023