# Simulation of a reinforced random walk

Hello!
I tried to stimulate an edge reinforced random walk for ants :

``````K = 20 #number of trajectory
N = 500 # time of one simulation

reinforcement = function(alpha,c,Zn){
return((1+(Zn/c))^alpha)
}

alpha = 10
c = 0.1
for (i in 1:K){
chose_right= c()
nbr_right = 0
for (j in 1:N){
proba = reinforcement(alpha,c,nbr_right)/(reinforcement(alpha,c,nbr_right)+reinforcement(alpha,c,j-nbr_right))
wetake_right = runif(1)<proba
chose_right[j]=wetake_right
nbr_right=nbr_right+wetake_right
}
plot(cumsum(nbr_right))
}

``````

The idea is that, in a bifurcation, the ants have to choose between the left and right direction.

The probability to choose the right direction is given by : Xi = 1 for the right direction
Xi = 0 for the left

Zn is number of time the right path have been chosen after the passage of n ants;

c represent the quantity of pheromones the ants depose when they select a path;

alpha is a positive parameter ;

I want to draw a graph which shows that the ants select(or not) a direction over time and this selection is rapid or slow depending on the values for the parameters alpha and c .

How to see that?

Does someone have an idea?

Thank you for you help!

Hello @Kevser_Cifci and welcome to RStudio Community        I spent a bit of time looking at your question and I must say it is fairly interesting. Before I share my code; however, I must say that your code to compute reinforcement is quite different from the formula you provided. Maybe I am missing something? In my code, I try to replicate the formula instead, but you can always change that part if you want.

I rewrote the entire code so I could understand the process a bit more. Also, even after thinking a lot about it, I still do not understand the role of `K`. My code below does not incorporate it, but I think it still captures the essence of your question. I may be wrong, please let me know what you think when you see and peruse the code.

Regarding the visualizations, I made 4 of them. First, I choose a fixed value of `c`, then I run the process for increasing values of `alpha`. I do it 500 times then I plot the `Zn` lines by color-coding them. Yellow means a low `alpha` and red means a high `alpha`. I was expecting yellow lines to be at the bottom (resp. at the top) and red lines to be at the top (resp. at the bottom), but there is no pattern at all. This means that increasing the positive parameter `alpha` does not necessarily cause ants to choose right more (or less) often. Hence, the process is indeed a random walk.

As a way to show this better, I also take "a slice" of the visualization at `n = 500` (represented by `i` in the plot above) and show it in a scatterplot. Again, there is no clear pattern.

Finally, I made two similar plots by choosing a fixed value of `alpha` as well as increasing values of `c`. The results are similar in that no pattern is apparent.

I would like to end by reiterating that the formula I used to calculate the probability is different from what you use in your code. My formula is an implementation of the formula you provided in your question.

Here is the full code:

``````# Load packages ----

library(dplyr)
library(ggplot2)
library(purrr)
library(stringr)
library(tibble)
library(tidyr)

# Define simulation function ----

run_simulation <- function(N, alpha, c){
Zn <- 0
Zi <- c()

for(n in seq_len(N)){
p <- (c + Zn)^alpha / ((c + Zn)^alpha + (c + n - Zn)^alpha)
take_right <- rbinom(n = 1, size = 1, prob = p)
Zi[n] <- take_right
Zn <- Zn + take_right
}

cumsum(Zi)
}

# SIMULATE OVER ALPHA ------------------------------------

## Simulate over the alpha parameter for a fixed value of c (c = 100)

alpha_values <- seq(from = 0.001, to = 0.5, length.out = 250)

alpha_sims <- imap_dfc(
alpha_values,
function(alpha, i) run_simulation(N = 500, alpha = alpha, c = 100) %>% enframe(name = NULL)
) %>%
mutate(i = seq_len(nrow(.))) %>%
rename_with(.fn = ~ c(paste0("alpha_", alpha_values), "i")) %>%
pivot_longer(cols = -i, names_to = "alpha", values_to = "Z") %>%
mutate(alpha = str_remove(string = alpha, pattern = "alpha_") %>% as.numeric()) %>%
arrange(alpha, i) %>%
select(alpha, i, Z)

## Visualize simulations

### Line charts of the evolution of Z (Z1 -> Z500) for all values of alpha

ggplot(data = alpha_sims) +
geom_line(mapping = aes(x = i, y = Z, group = alpha, color = alpha)) +
scale_color_gradient(low = "yellow", high = "red") +
theme_bw()

ggsave(filename = here::here("visualizations/kevser_cifci_1.png"))

### Comparison of Z500 for all values of alpha

alpha_sims %>%
filter(i == 500) %>%
ggplot(mapping = aes(x = alpha, y = Z, color = alpha)) +
geom_point() +
scale_color_gradient(low = "yellow", high = "red") +
theme_bw()

ggsave(filename = here::here("visualizations/kevser_cifci_2.png"))

# SIMULATE OVER C ------------------------------------

## Simulate over the C parameter for a fixed value of alpha (alpha = 0.005)

c_values <- seq(from = 1, to = 500, length.out = 250)

c_sims <- imap_dfc(
c_values,
function(c, i) run_simulation(N = 500, alpha = 0.005, c = c) %>% enframe(name = NULL)
) %>%
mutate(i = seq_len(nrow(.))) %>%
rename_with(.fn = ~ c(paste0("c_", c_values), "i")) %>%
pivot_longer(cols = -i, names_to = "c", values_to = "Z") %>%
mutate(c = str_remove(string = c, pattern = "c_") %>% as.numeric()) %>%
arrange(c, i) %>%
select(c, i, Z)

## Visualize simulations

### Line charts of the evolution of Z (Z1 -> Z500) for all values of c

ggplot(data = c_sims) +
geom_line(mapping = aes(x = i, y = Z, group = c, color = c)) +
scale_color_gradient(low = "yellow", high = "red") +
theme_bw()

ggsave(filename = here::here("visualizations/kevser_cifci_3.png"))

### Comparison of Z500 for all values of c

c_sims %>%
filter(i == 500) %>%
ggplot(mapping = aes(x = c, y = Z, color = c)) +
geom_point() +
scale_color_gradient(low = "yellow", high = "red") +
theme_bw()

ggsave(filename = here::here("visualizations/kevser_cifci_4.png"))

``````
1 Like

Hello!

Thank you for your reply! It's very interesting and I think your simulation function is better than what I've tried to code.

In theory, there are 4 cases :

• when "c < 1" and "alpha>1 " : the ants select fastly a direction.
• when "c>1" and "alpha>1" : the ants select slowly a direction.
• when "c<1" and "alpha<1" : after a big number of passages, both directions are chosen equally
• when "c>1" and "alpha <1" : after a small number of passages, both directions are chosen equally.

I don't know if this theory is visible when we plot Xi or Zi according to the number of passages?

I think you're idea of plotting the process according to the values of c and alpha is interesting even if it doesn't show the results I want.

I'm very sorry for my late answer.

Kevser

This topic was automatically closed 21 days after the last reply. New replies are no longer allowed.

If you have a query related to it or one of the replies, start a new topic and refer back with a link.