Analyse de la sensibilité de la sortie d'un automate cellulaire

Analyse de la sensibilité d'un automate cellulaire

# Fonction qui calcule le nombre de cases occupées (=1) en fonction d'une variable X
# qui est un data.frame qui contient les probabilités de mourir à chaque pas de temps (p) et 
# le nombre de cases (/9) pour obtenir un décès d'étouffement.
# La variable t donne les temps pour lesquels on retourne la taille de la population
# Si t est une valeur, on retourne juste la valeur pour ce temps, si t est un vecteur, 
# on retourne toutes les valeurs pour les temps t

ac <- function(X, t=t) {

  Nfinal <- NULL

  for (grandtotal in 1:nrow(X)) {

    penc <- X[grandtotal, 1]
    v <- X[grandtotal, 2]*9
 
  m <- matrix(data = 0, ncol=12, nrow=12)
  for (i in 1:10) {
    col <- sample(1:10, 1)+1
    row <- sample(1:10, 1)+1
    m[row, col] <- 1
  }

  N <- NULL

  for (tps in 1:max(t)) {
    m1 <- m
 
    for (col in 2:11) {
      for (row in 2:11) {
        voisins <- sum(m[(row-1):(row+1), (col-1):(col+1)])-m[row, col]
        if (m[row, col] == 1) {
          if (runif(1)>1-penc) m1[row, col] <- 0
          if ((voisins == 0) | (voisins > v)) m1[row, col] <- 0
        }
        if (m[row, col] == 0) {
          if (voisins >= 2) m1[row, col] <- 1
        }
      }
    }
    m <- m1
    N <- c(N, sum(m))
  }

  Nfinal <- c(Nfinal, N[t])

  }

  return(Nfinal)
}

# Je retourne les valeurs de la taille de la population pour t=1 à 10 avec une probalité de mourir de 0.3
ac(X=data.frame(p=0.3, v=6/9), t=1:10)
# Je retourne les valeurs de la taille de la population pour t=1000 avec une probabilité de 
# mourir de 0.3 et un nombre de case occupées pour mourir d'étouffement de 6 ou plus
ac(X=data.frame(p=0.3, v=6/9), t=1000)
# Je retourne les valeurs de la taille de la population pour t=1000 avec 
# une probabilité de mourir de 0.4 et de 0.3 et un nombre de case occupées 
# pour mourir d'étouffement de 6 ou 7 (ou plus)
ac(X=data.frame(p=c(0.4, 0.3), v=c(6, 7)/9), t=1000)

# Si j'ai une probabilité de mourir de 0, le système est périodique
# mais est différent selon l'état initial
plot(x = 1:100, ac(X=data.frame(p=0, v=6/9), t=1:100), type="l", bty="n", las=1)
# J'introduis de la stochasticité; c'est plus amusant !
plot(x = 1:100, ac(X=data.frame(p=0.3, v=6/9), t=1:100), type="l", bty="n", las=1)


library("sensitivity")
# Je génère deux tableaux identiques avec des nombres aléatoires
# p est la probabilité de mourir
# v est le nombre de voisins /9 (pour rester dans l'intervalle 0-1) maximum sinon on meurt
X1 <- data.frame(p=runif(100, 0.1, 0.6), v=runif(100, 0.5, 0.9))
X2 <- data.frame(p=runif(100, 0.1, 0.6), v=runif(100, 0.5, 0.9))


xjansen100 <- soboljansen(model = ac, X1 = X1, X2 = X2,
                                           nboot = 1000, t=100)
# La conclusion est que les variables agissent surtout en synergie l'une de l'autre
# ce qui est normal car les deux jouent sur la mortalité.
# v a un impact plus important que p qui ne joue presque pas seul mais beaucoup avec v.
plot(xjansen100)