Devoir 1 : Dynamique de populations

Question 0 : Intuition

Par intuition, il nous semble que pour les 10 trajectoires l’évolution de P(t) en partant de l’état initial (n1(1), n2(1)) = (1, 1), dans les presque demi de trajectoires, P(t) vont tendre vers 1, et dans les autres demi de trajectoires, P(t) vont tendre vers 0. La raison de notre intuition est que tout d’abord, le premier immigrant choisit la ville 1 avec la probabilité ½, si il a choisi la ville 1, P(1) = ⅔, donc probablement le deuxième immigrant va choisir la ville 1, si c’est le cas, P(2) ¾. Donc avec l’augmentation de la population de la ville 1, la possibilité de recevoir le nouvel immigrant est de plus en plus grande. A l’inverse, si il a choisit la ville 2, P(1) = ⅓, donc probablement le deuxième immigrant va choisir la ville 2, si c’est le cas, P(2) = ¼. Donc avec l’augmentation de la population de la ville 2, la possibilité de recevoir le nouvel immigrant pour la ville 1 est de plus en plus petite.

Par intuition, il nous semble que dans l’histogramme de la valeur limite de P(t) pour 1000 trajectoires, la fréquence de la valuer 0 sera presque 500 et la fréquence de la valeur 1 sera preseque 500.

Par intuition, il nous semble que pour les états initiaux (n1(1), n2(1)) 𝞊{(1, 1), (10, 10), (100, 100),} la vitesse d’évolution de P(t) sera de plus en plus lente, mais les histogrammes des limites de trajectoires seront les mêmes. Et pour les états initiaux (n1(1), n2(1)) 𝞊{(1, 1), (1, 10), (1, 100),}, dans les histogrammes des limites de trajectoires, la fréquence de la valuer 0 sera presque 1000. En conclusion, par intuition, on pense que pour la ville dont la population est plus grande la possibilité P(t) va être de plus en plus grande, et pour la ville dont la population est moins grande la possibilité P(t) va être de plus en plus petite.

Question 1 : Premières observations de la dynamique

D'abord, on fixe la valeur maximale de t à 100,

rm(list = ls())

set.seed(2)

trajectoire <- 10

n1_ini <- 1
n2_ini <- 1

t <- 100

p <- rep(0, t + 1)
p[1] <- n1_ini/(n1_ini + n2_ini)

i <- 0
while (i < trajectoire) {
    j <- 0

    n1 <- n1_ini
    n2 <- n2_ini

    while (j < t) {

        if (runif(1) < n1/(n1 + n2)) {
            n1 <- n1 + 1
        } else {
            n2 <- n2 + 1
        }

        j <- j + 1
        p[j + 1] <- n1/(n1 + n2)

    }
    if (i != 0) {
        par(new = TRUE)
    }
    plot(p, type = "l", ylim = c(0, 1), col = i, main = "Courbes de P(t)", col.main = "blue", 
        xlab = "t", ylab = "P(t)")
    i <- i + 1
}

plot of chunk unnamed-chunk-1

On a trouvé que p(t) change beaucoup entre 0 et 100, et il a l'air que p(100) est possible d'égaler à tous les valeurs entre 0 et 1.

Pour comprendre mieux, on a choisi une valeur assez grand 50000 comme la valeur maximale de t, et on a fixé le nombre de trajectoires à 20.

rm(list = ls())

set.seed(3)

trajectoire <- 20

n1_ini <- 1
n2_ini <- 1

t <- 50000

p <- rep(0, t + 1)
p[1] <- n1_ini/(n1_ini + n2_ini)

i <- 0
while (i < trajectoire) {
    j <- 0

    n1 <- n1_ini
    n2 <- n2_ini

    while (j < t) {

        if (runif(1) < n1/(n1 + n2)) {
            n1 <- n1 + 1
        } else {
            n2 <- n2 + 1
        }

        j <- j + 1
        p[j + 1] <- n1/(n1 + n2)

    }
    if (i != 0) {
        par(new = TRUE)
    }
    plot(p, type = "l", ylim = c(0, 1), col = i, main = "Courbes de P(t)", col.main = "blue", 
        xlab = "t", ylab = "P(t)")
    i <- i + 1
}

plot of chunk unnamed-chunk-2

C'est évident que p(t) n'est pas bien stable entre 0 et 5000, quand t>20000 p(t) reste stable, et ça nous démontre encore une fois, c'est possible que p(t) sera enfin égal à tous les valeurs entre 0 et 1.

Question 2 : Loi statistique de la limite

Selon la dernière question, on a trouvé que p(t) sera très stable quand t > 20000, donc dans les programme suivants, on a choisi la valeur de t comme 20000 pour caleculer la limite de p(t).

rm(list = ls())

set.seed(2)

trajectoire <- 1000

n1_ini <- 1
n2_ini <- 1

t <- 20000

p_limit <- rep(0, trajectoire)

i <- 0
while (i < trajectoire) {
    j <- 0

    n1 <- n1_ini
    n2 <- n2_ini

    while (j < t) {

        if (runif(1) < n1/(n1 + n2)) {
            n1 <- n1 + 1
        } else {
            n2 <- n2 + 1
        }

        j <- j + 1
    }
    p_limit[i + 1] <- n1/(n1 + n2)
    i <- i + 1
}

hist(p_limit, , main = "Histogram of p_limit, n1=1, n2=1")

plot of chunk unnamed-chunk-3

Selon la histogramme, on a trouvé que c'est possible que p_limit tombe dans tous les intervalles ,et les probabilités que p(t) tombe dans chaque intervalle sont prèsque égales. On pense que ça suit la loi uniforme.

Question 3 : Influence de la répartition initiale

rm(list = ls())

set.seed(2)

trajectoire <- 1000

n1_ini <- 10
n2_ini <- 10

t <- 20000

p_limit <- rep(0, trajectoire)

i <- 0
while (i < trajectoire) {
    j <- 0

    n1 <- n1_ini
    n2 <- n2_ini

    while (j < t) {

        if (runif(1) < n1/(n1 + n2)) {
            n1 <- n1 + 1
        } else {
            n2 <- n2 + 1
        }

        j <- j + 1
    }
    p_limit[i + 1] <- n1/(n1 + n2)
    i <- i + 1
}

hist(p_limit, main = "Histogram of p_limit, n1=10, n2=10")

plot of chunk unnamed-chunk-4

rm(list = ls())

set.seed(2)

trajectoire <- 1000

n1_ini <- 100
n2_ini <- 100

t <- 20000

p_limit <- rep(0, trajectoire)

i <- 0
while (i < trajectoire) {
    j <- 0

    n1 <- n1_ini
    n2 <- n2_ini

    while (j < t) {

        if (runif(1) < n1/(n1 + n2)) {
            n1 <- n1 + 1
        } else {
            n2 <- n2 + 1
        }

        j <- j + 1
    }
    p_limit[i + 1] <- n1/(n1 + n2)
    i <- i + 1
}

hist(p_limit, main = "Histogram of p_limit, n1=100, n2=100")

plot of chunk unnamed-chunk-5

Comparé avec la histogramme dans la question 2, p_limit a une tendance à tendre vers 0.5. quand la valeur initiale de n1 et n2 est plus grand, il est plus évident.

rm(list = ls())

set.seed(2)

trajectoire <- 1000

n1_ini <- 1
n2_ini <- 10

t <- 20000

p_limit <- rep(0, trajectoire)

i <- 0
while (i < trajectoire) {
    j <- 0

    n1 <- n1_ini
    n2 <- n2_ini

    while (j < t) {

        if (runif(1) < n1/(n1 + n2)) {
            n1 <- n1 + 1
        } else {
            n2 <- n2 + 1
        }

        j <- j + 1
    }
    p_limit[i + 1] <- n1/(n1 + n2)
    i <- i + 1
}

hist(p_limit, main = "Histogram of p_limit, n1=1, n2=10")

plot of chunk unnamed-chunk-6

rm(list = ls())

set.seed(2)

trajectoire <- 1000

n1_ini <- 1
n2_ini <- 100

t <- 20000

p_limit <- rep(0, trajectoire)

i <- 0
while (i < trajectoire) {
    j <- 0

    n1 <- n1_ini
    n2 <- n2_ini

    while (j < t) {

        if (runif(1) < n1/(n1 + n2)) {
            n1 <- n1 + 1
        } else {
            n2 <- n2 + 1
        }

        j <- j + 1
    }
    p_limit[i + 1] <- n1/(n1 + n2)
    i <- i + 1
}

hist(p_limit, main = "Histogram of p_limit, n1=1, n2=100")

plot of chunk unnamed-chunk-7

Dans ces deus histogrammes, p_limit a une tendance à tendre vers 0. Quand n2 est plus grand, il est plus évident.

Question 4 : Avez-vous changé d'avis?

Après d'avoir vu le résultat, on a trouvé que c'est différent que notre intuition.

Des autres Hypothèse: