Evoluation d’une population (correction)

Links: notebook, html ., PDF, python, slides ., presentation ., GitHub

Evolution d’une population à partir des tables de mortalités et d’une situation initiale.

%pylab inline
import matplotlib.pyplot as plt
plt.style.use('ggplot')
Populating the interactive namespace from numpy and matplotlib
from jyquickhelper import add_notebook_menu
add_notebook_menu()
run previous cell, wait for 2 seconds

Exercice 1 : pyramides des âges

from actuariat_python.data import population_france_year
population = population_france_year()
df = population
df.head(n=3)
naissance age hommes femmes ensemble
0 2016 0 382585 364392 746977
1 2015 1 390810 373164 763974
2 2014 2 402728 386467 789195
hommes = df["hommes"]
femmes = df["femmes"]
somme = hommes - femmes

Je reprends ici le code exposé à Damien Vergnaud’s Homepage en l’adaptant un peu avec les fonctions de matplotlib via l’interface pyplot. Puis j’ajoute la différence par âge. On commence souvent par la gallerie pour voir si un graphe ou juste une partie est similaire à ce qu’on veut obtenir.

from matplotlib import pyplot as plt
from numpy import arange
plt.style.use('ggplot')
fig, ax = plt.subplots(figsize=(8,8))
ValH = ax.barh(arange(len(hommes)),hommes,1.0,label="Hommes",color='b',linewidth=0,align='center')
ValF = ax.barh(arange(len(femmes)),-femmes,1.0,label="Femmes",color='r',linewidth=0,align='center')
diff, = ax.plot(somme,arange(len(femmes)),'y',linewidth=2)
ax.set_title("Pyramide des âges")
ax.set_ylabel("Ages")
ax.set_xlabel("Habitants")
ax.set_ylim([0,110])
ax.legend((ValH[0],ValF[0],diff),('Hommes','Femmes','différence'))
<matplotlib.legend.Legend at 0x1e44c4a0c88>
../_images/seance4_projection_population_correction_7_1.png

Le même en utilisant la fonction insérée dans le module actuariat_python.

from actuariat_python.plots import plot_population_pyramid
plot_population_pyramid(df["hommes"], df["femmes"], figsize=(8,4))
<matplotlib.axes._subplots.AxesSubplot at 0x1e44bc09fd0>
../_images/seance4_projection_population_correction_9_1.png

Exercice 2 : calcul de l’espérance de vie

Le premier objectif est de calculer l”espérance de vie à l’âge t à partir de la table de mortalité. On récupère cette table.

from actuariat_python.data import table_mortalite_france_00_02
df=table_mortalite_france_00_02()
import pandas
pandas.concat([ df.head(n=3), df.tail(n=3) ])
Age Homme Femme
0 0 100000.0 100000.0
1 1 99511.0 99616.0
2 2 99473.0 99583.0
115 117 0.0 0.0
116 118 0.0 0.0
117 120 0.0 0.0

On note P_t la population l’âge t. La probabilité de mourir à la date t+d lorsqu’on a l’âge t correspond à la probabilité de rester en vie à jusqu’à l’âge t+d puis de mourir dans l’année qui suit :

m_{t+d} = \frac{P_{t+d}}{P_t}\frac{P_{t+d} - P_{t+d+1}}{P_{t+d}}

.

L’espérance de vie s’exprime :

\mathbb{E}(t) = \sum_{d=1}^\infty d m_{t+d} = \sum_{d=1}^\infty d \frac{P_{t+d}}{P_t}\frac{P_{t+d} - P_{t+d+1}}{P_{t+d}} = \sum_{d=1}^\infty d \frac{P_{t+d} - P_{t+d+1}}{P_{t}}

On crée une matrice allant de 0 à 120 ans et on pose \mathbb{E}(120)=0. On utilise le module numpy.

import numpy
hf = df[["Homme", "Femme"]].as_matrix()
hf = numpy.vstack([hf, numpy.zeros((8,2))])
hf.shape
(126, 2)
nb = hf.shape[0]
esp = numpy.zeros ((nb,2))
for t in range(0,nb):
    for i in (0,1):
        if hf[t,i] == 0:
            esp[t,i] = 0
        else:
            somme  = 0.0
            for d in range(1,nb-t):
                if hf[t+d,i] > 0:
                    somme += d * (hf[t+d,i] - hf[t+d+1,i]) / hf[t,i]
            esp[t,i] = somme
esp[:1]
array([[ 75.00752,  82.48832]])

Enfin, on dessine le résultat avec matplotlib :

from matplotlib import pyplot as plt
plt.style.use('ggplot')
h = plt.plot(esp)
plt.legend(h, ["Homme", "Femme"])
plt.title("Espérance de vie")
<matplotlib.text.Text at 0x1e44ca6b2b0>
../_images/seance4_projection_population_correction_17_1.png

Le calcul implémenté ci-dessus n’est pas le plus efficace. On fait deux boucles imbriquées dont le coût global est en O(n^2) mais surtout on effectue les mêmes calculs plusieurs fois. Pour le réduire à un coût linéaire O(n), il faut s’intéresser à la quantité :

P_{t+1} \mathbb{E}(t+1) - P_t \mathbb{E}(t) = \sum_{d=1}^\infty d (P_{t+d+1} - P_{t+d+2}) - \sum_{d=1}^\infty d (P_{t+d} - P_{t+d+1})

L’implémentation devra utiliser la fonction numpy.cumsum et cette astuce Pandas Dataframe cumsum by row in reverse column order?.

# à suivre

numpy et pandas ont plusieurs fonction en commun dès qu’il s’agit de parcourir les données. Il existe aussi la fonction DataFrame.cumsum.

Exercice 3 : simulation de la pyramide l’année suivante

L’objectif est d’estimer la population française à l’année T. Si P(a,T) désigne le nombre de personnes d’âge a en T, on peut estimer P(a,T+1) en utilisant la probabilité de mourir m(a) :

P(a+1, T+1) = P(a,T) * (1 - m(a))

On commence par calculer les coefficients m(a) avec la table hf obtenue lors de l’exercice précédent tout en gardant la même dimension (on aura besoin de la fonction nan_to_num :

mortalite = (hf[:-1] - hf[1:]) / hf[:-1]
mortalite = numpy.nan_to_num(mortalite)  # les divisions nulles deviennent nan, on les remplace par 0
mortalite = numpy.vstack([mortalite, numpy.zeros((1,2))])
m = mortalite
c:Python36_x64libsite-packagesipykernel_launcher.py:1: RuntimeWarning: invalid value encountered in true_divide
  """Entry point for launching an IPython kernel.

La population a été obtenue lors de l’exercice 1, on la convertit en un objet numpy :

pop = population[["hommes","femmes"]].as_matrix()
pop = numpy.vstack( [pop, numpy.zeros((m.shape[0] - pop.shape[0],2))])
pop.shape
(126, 2)

Ensuite on calcule la population en 2016 :

pop_next = pop * (1-m)
pop_next = numpy.vstack([numpy.zeros((1,2)), pop_next[:-1]])
pop_next[:5]
array([[      0.        ,       0.        ],
       [ 380714.15935   ,  362992.73472   ],
       [ 390660.76242827,  373040.38118375],
       [ 402618.6873624 ,  386385.50208369],
       [ 405546.26293667,  387030.90400956]])
pop[:5]
array([[ 382585.,  364392.],
       [ 390810.,  373164.],
       [ 402728.,  386467.],
       [ 405636.,  387097.],
       [ 416074.,  396365.]])
from actuariat_python.plots import plot_population_pyramid
plot_population_pyramid(pop_next[:,0], pop_next[:,1])
<matplotlib.axes._subplots.AxesSubplot at 0x1e44c903400>
../_images/seance4_projection_population_correction_29_1.png

Exercice 4 : simulation jusqu’en 2100

Il s’agit de répéter l’itération effectuée lors de l’exercice précédent. Le plus est de recopier le code dans une fonction et de l’appeler un grand nombre de fois.

def iteration(pop, mortalite):
    pop_next = pop * (1-mortalite)
    pop_next = numpy.vstack([numpy.zeros((1,2)), pop_next[:-1]]) # aucune naissance
    return pop_next

popt = pop
for year in range(2016, 2051):
    popt = iteration(popt, mortalite)

plot_population_pyramid(popt[:,0], popt[:,1], title="Pyramide des âges en 2050")
<matplotlib.axes._subplots.AxesSubplot at 0x1e44cbbb9b0>
../_images/seance4_projection_population_correction_31_1.png

Exercice 5 : simulation avec les naissances

Dans l’exercice précédent, la seconde ligne de la fonction iteration correspond à cas où il n’y a pas de naissance. On veut remplacer cette ligne par quelque chose proche de la réalité :

  • les naissances sont calculées à partir de la population féminines et de la table de fécondité
  • on garde la même proportion homme/femme que celle actuellement observée
ratio = pop[0,0] / (pop[0,1] + pop[0,0])
ratio
0.51217775112218988

Il y a un peu plus de garçons qui naissent chaque année.

from actuariat_python.data import fecondite_france
df=fecondite_france()
df.head()
age 2005 2015
3 15.0 0.1 0.0
4 16.0 0.2 0.2
5 17.0 0.5 0.5
6 18.0 1.0 1.1
7 19.0 2.0 2.1
from matplotlib import pyplot as plt
df.plot(x="age", y=["2005","2015"])
<matplotlib.axes._subplots.AxesSubplot at 0x1e44d2e3ef0>
../_images/seance4_projection_population_correction_36_1.png

On convertit ces données en une matrice numpy sur 120 lignes comme les précédentes. On se sert des méthodes fillna et merge.

ages = pandas.DataFrame(dict(age=range(0,120)))
merge = ages.merge(df, left_on="age", right_on="age", how="outer")
fecondite = merge.fillna(0.0)
fecondite[13:17]
age 2005 2015
13 13 0.0 0.0
14 14 0.0 0.0
15 15 0.1 0.0
16 16 0.2 0.2
mat_fec = fecondite[["2015"]].as_matrix() / 100 # les chiffres sont pour 100 femmes
mat_fec.shape
(120, 1)
mat_fec.sum()
1.899

Si la matrice pop a plus de ligne que la matrice mat_fec, on doit compléter la seconde avec autant de lignes nulles que la précédente.

if mat_fec.shape[0] < pop.shape[0]:
    zeros = numpy.zeros((pop.shape[0] - mat_fec.shape[0], mat_fec.shape[1]))
    mat_fec = numpy.vstack([mat_fec, zeros])
mat_fec.sum()
1.899

Il faut maintenant coder une fonction qui calcule le naissances pour l’année suivantes.

def naissances(pop, fec):
    # on suppose que pop est une matrice avec deux colonnes homme, femme
    # et que fec est une matrice avec une colonne fécondité
    n = pop[:,1] * fec[:,0]
    return n.sum()

nais = naissances(pop, mat_fec)
nais
770829.48399999994

Et on reprend la fonction iteration et le code de l’exercice précédent :

def iteration(pop, mortalite, fec, ratio):
    pop_next = pop * (1-mortalite)
    nais = naissances(pop, fec)
    row = numpy.array([[nais*ratio, nais*(1-ratio)]])
    pop_next = numpy.vstack([row, pop_next[:-1]]) # aucune naissance
    return pop_next

popt = pop
for year in range(2016, 2051):
    popt = iteration(popt, m, mat_fec, ratio)

plot_population_pyramid(popt[:,0], popt[:,1], title="Pyramide des âges en 2050")
<matplotlib.axes._subplots.AxesSubplot at 0x1e44d5a5940>
../_images/seance4_projection_population_correction_46_1.png

On va plus loin et on stocke la population dans un vecteur :

total = [[2015, pop[:,0].sum(),pop[:,1].sum()]]
popt = pop
for year in range(2016, 2101):
    popt = iteration(popt, m, mat_fec, ratio)
    total.append([year, popt[:,0].sum(),popt[:,1].sum()])

plot_population_pyramid(popt[:,0], popt[:,1], title="Pyramide des âges en 2101")
<matplotlib.axes._subplots.AxesSubplot at 0x1e44d6137f0>
../_images/seance4_projection_population_correction_48_1.png
df = pandas.DataFrame(data=total, columns=["année","hommes","femmes"])
df.plot(x="année", y=["hommes", "femmes"], title="projection population française")
<matplotlib.axes._subplots.AxesSubplot at 0x1e44d51b9b0>
../_images/seance4_projection_population_correction_49_1.png

Le code suivant permet de combiner les deux graphes sur la même ligne avec la fonction subplots :

from matplotlib import pyplot as plt
fig, ax = plt.subplots(1, 2, figsize=(14,6))
plot_population_pyramid(popt[:,0], popt[:,1], title="Pyramide des âges en 2050", ax=ax[0])
df.plot(x="année", y=["hommes", "femmes"], title="projection population française", ax=ax[1])
<matplotlib.axes._subplots.AxesSubplot at 0x1e44e202630>
../_images/seance4_projection_population_correction_51_1.png

Autre représentation d’une pyramide avec le module pygal

Le code suivant utilise le module pygal et le notebook Interactive plots in IPython notebook using pygal ou celui-ci pygal.ipynb.

from actuariat_python.data import population_france_year
population = population_france_year()
population.head(n=3)
naissance age hommes femmes ensemble
0 2016 0 382585 364392 746977
1 2015 1 390810 373164 763974
2 2014 2 402728 386467 789195
from IPython.display import SVG
import pygal
pyramid_chart = pygal.Pyramid(human_readable=True, legend_at_bottom=True)
pyramid_chart.title = 'Population française en T+1'
pyramid_chart.x_labels = map(lambda x: str(x) if not x % 5 else '', population["age"])
pyramid_chart.add("hommes", population["hommes"])
pyramid_chart.add("femmes", population["femmes"])
SVG(pyramid_chart.render())
../_images/seance4_projection_population_correction_54_0.svg