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_2015
population = population_france_2015()
df = population
df.head(n=3)
naissance age hommes femmes ensemble
8 2015 0 391371 374179 765550
9 2014 1 403204 385442 788646
10 2013 2 405502 386831 792333
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 0xafba320>
../_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 0xac51e48>
../_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 100000
1 1 99511 99616
2 2 99473 99583
109 109 2 20
110 110 1 9
111 111 0 4

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
(120, 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.48836]])

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 0xb7f36a0>
../_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 en 2016

L’objectif est d’estimer la population française en 2016. Si P(a,2015) désigne le nombre de personnes d’âge a en 2015, on peut estimer P(a,2016) en utilisant la probabilité de mourir m(a) :

P(a+1, 2016) = P(a,2015) * (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

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
(120, 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.        ],
       [ 389457.19581   ,  372742.15264   ],
       [ 403050.02956457,  385314.31382509],
       [ 405391.93441436,  386749.4253236 ],
       [ 412291.77032761,  391786.09193267]])
pop[:5]
array([[ 391371.,  374179.],
       [ 403204.,  385442.],
       [ 405502.,  386831.],
       [ 412383.,  391853.],
       [ 416626.,  399632.]])
from actuariat_python.plots import plot_population_pyramid
plot_population_pyramid(pop_next[:,0], pop_next[:,1])
<matplotlib.axes._subplots.AxesSubplot at 0xb034c88>
../_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 0xb8f7a90>
../_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.5112285285089152

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.8 0.8
4 16 2.3 2.3
5 17 5.8 5.0
6 18 11.9 9.5
7 19 22.3 17.2
from matplotlib import pyplot as plt
df.plot(x="age", y=["2005","2015"])
<matplotlib.axes._subplots.AxesSubplot at 0xcfab908>
../_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.8 0.8
16 16 2.3 2.3
mat_fec = fecondite[["2015"]].as_matrix() / 1000 # les chiffres sont pour 1000 femmes
mat_fec.shape
(120, 1)

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
801294.71230000001

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 0xbc4a6a0>
../_images/seance4_projection_population_correction_43_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 0xd057828>
../_images/seance4_projection_population_correction_45_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 0xd3d8080>
../_images/seance4_projection_population_correction_46_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 0xd7ed0f0>
../_images/seance4_projection_population_correction_48_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_2015
population = population_france_2015()
population.head(n=3)
naissance age hommes femmes ensemble
8 2015 0 391371 374179 765550
9 2014 1 403204 385442 788646
10 2013 2 405502 386831 792333
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 2015'
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_51_0.svg