1A.e - Enoncé 12 décembre 2017 (1)#

Correction du premier énoncé de l’examen du 12 décembre 2017. Celui-ci mène à l’implémentation d’un algorithme qui permet de retrouver une fonction $f$ en escalier à partir d’un ensemble de points $(X_i, f(X_i))$ .

Q1 - échantillon aléatoire#

Générer un ensemble aléatoire de 1000 nombres $(X_i,Y_i)$ qui vérifie :

$X_i$ suit une loi uniforme sur $[0,16]$
$Y_i = \mathbb{1}_{[\sqrt{X_i}] \mod 2}$ où $[A]$ est la partie entière de $A$ .

On pourra se servir de la fonction random du module random.

[2]:

import random

X = [random.random() * 16 for i in range(0, 1000)]
Y = [int(x**0.5) % 2 for x in X]

Q1 - dessiner le nuage de points - donnée#

Le code suivant vous est donné afin de vérifier vos réponses.

[3]:

%matplotlib inline

[4]:

import matplotlib.pyplot as plt

plt.plot(X, Y, ".")

[4]:

[<matplotlib.lines.Line2D at 0x1f893f9e518>]

../../_images/practice_exams_td_note_2018_1_5_1.png

Q2 - tri#

Trier les points selon les $X$ .

[5]:

nuage = [(x, y) for x, y in zip(X, Y)]
nuage.sort()
nuage[:5]

[5]:

[(0.014962888038782651, 0),
 (0.020462778257442693, 0),
 (0.022310859639508962, 0),
 (0.03078728731371605, 0),
 (0.03153252863972433, 0)]

Q3 - moyenne#

On suppose que les $Y$ sont triés selon les $X$ croissants. Calculer la moyenne des différences entre $Y$ et la moyenne $m$ des $Y$ (en valeur absolue) sur un intervalle $[i,j]$ , $j$ exclu. Ecrire une fonction def somme_diff(nuage, i, j) qui exécute ce calcul qui correspond à $\sum_{k=i}^{j-1} |Y_k - m|$ avec $m = (\sum_{k=i}^{j-1} Y_k) / (j-i)$ .

[6]:

def somme_diff(xy, i, j):
    m = sum(e[1] for e in xy[i:j]) / (j - i)
    return sum(abs(e[1] - m) for e in xy[i:j])


somme_diff(nuage, 0, 5), somme_diff(nuage, 0, len(nuage))

[6]:

(0.0, 476.2380000000092)

Q4 - distance#

Soit $i,j$ deux entiers, on coupe l’intervalle en deux : $i,k$ et $k,j$ . On calcule la somme_diff sur ces deux intervalles, on fait la somme des différences (en valeurs absolues) de ces moyennes par rapport à la valeur sur le plus grand intervalle. On écrit la fonction def difference(nuage, i, j, k):.

[7]:

def difference(nuage, i, j, k):
    m1 = somme_diff(nuage, i, k)
    m2 = somme_diff(nuage, k, j)
    m = somme_diff(nuage, i, j)
    return abs(m1 + m2 - m)


difference(nuage, 0, len(nuage), 100)

[7]:

18.56022222223197

Q5 - fonction comme paramètre#

Le langage Python permet de passer une fonction à une autre fonction en tant qu’argument. Un exemple :

[8]:

def fct(x, y):
    return abs(x - y)

def distance_list(list_x, list_y, f):
    return sum(f(x, y) for x, y in zip(list_x, list_y))

distance_list([0, 1], [0, 2], fct)

[8]:

Ecrire la fonction précédente en utilisant la fonction fct.

[9]:

def somme_diff(xy, i, j, f):
    m = sum(e[1] for e in xy[i:j]) / (j - i)
    # On a modifié les fonctions précédentes pour calculer
    # une fonction d'erreur "custom" ou définie par l'utilisateur.
    return sum(f(e[1], m) for e in xy[i:j])


def difference(nuage, i, j, k, f):
    m1 = somme_diff(nuage, i, k, f)
    m2 = somme_diff(nuage, k, j, f)
    m = somme_diff(nuage, i, j, f)
    return abs(m - m1) + abs(m - m2)


difference(nuage, 0, len(nuage), 100, fct)

[9]:

494.7982222222412

Q6 - optimiser#

On veut déterminer le $i$ optimal, celui qui maximise la différence dans l’intervalle $[i,j]$ . On souhaite garder la fonction fct comme argument. Pour cela, implémenter la fonction def optimise(nuage, i, j, f):.

[10]:

def optimise(nuage, i, j, f):
    mx = -1
    ib = None
    for k in range(i + 1, j - 1):
        d = difference(nuage, i, j, k, f)
        if ib is None or d > mx:
            mx = d
            ib = k
    if ib is None:
        # Au cas où l'intervalle est vide, on retourne une coupure
        # égale à i.
        ib = i
        mx = 0
    return ib, mx


optimise(nuage, 0, len(nuage), fct)

[10]:

(565, 711.6476814159435)

[11]:

import matplotlib.pyplot as plt

x = nuage[552][0]
plt.plot(X, Y, ".")
plt.plot([x, x], [0, 1])

[11]:

[<matplotlib.lines.Line2D at 0x1f893ea57b8>]

../../_images/practice_exams_td_note_2018_1_18_1.png

Le premier point de coupure trouvé (le trait orange) correspond à un des bords d’un des escaliers.

Q7 - optimisation encore#

Recommencer sur les deux intervalles trouvés. La question était juste histoire que le résultat à la question précédente est reproductible sur d’autres intervalles.

[12]:

optimise(nuage, 0, 68, fct), optimise(nuage, 68, len(nuage), fct)

[12]:

((1, 0.0), (565, 618.0710615624871))

[13]:

import matplotlib.pyplot as plt

x = nuage[58][0]
x2 = nuage[552][0]
plt.plot(X, Y, ".")
plt.plot([x, x], [0, 1])
plt.plot([x2, x2], [0, 1])

[13]:

[<matplotlib.lines.Line2D at 0x1f8943ece80>]

../../_images/practice_exams_td_note_2018_1_22_1.png

Q8 - fonction récursive#

Pouvez-vous imaginer une fonction récursive qui produit toutes les séparations. Ecrire la fonction def recursive(nuage, i, j, f, th=0.1):.

[14]:

def recursive(nuage, i, j, f, th=0.1):
    k, mx = optimise(nuage, i, j, f)
    if mx <= th:
        return None
    r1 = recursive(nuage, i, k, f, th=th)
    r2 = recursive(nuage, k, j, f, th=th)
    if r1 is None and r2 is None:
        return [k]
    elif r1 is None:
        return [k] + r2
    elif r2 is None:
        return r1 + [k]
    else:
        return r1 + [k] + r2


r = recursive(nuage, 0, len(nuage), fct)
r

[14]:

[68, 242, 565]

[15]:

import matplotlib.pyplot as plt

plt.plot(X, Y, ".")
for i in r:
    x = nuage[i][0]
    plt.plot([x, x], [0, 1])

../../_images/practice_exams_td_note_2018_1_25_0.png

Q9 - coût#

Quel est le coût de la fonction optimize en fonction de la taille de l’intervalle ? Peut-on mieux faire (ce qu’on n’implémentera pas).

Tel qu’il est implémenté, le coût est en $O(n^2)$ , le coût peut être linéaire en triant les éléments dans l’ordre croissant, ce qui a été fait, ou $n\ln n$ si on inclut le coût du tri bien qu’on ne le fasse qu’une fois. Voyons plus en détail comment se débarrasser du coût en $O(n^2)$ . Tout d’abord la version actuelle.

[16]:

def somme_diff_abs(xy, i, j):
    m = sum(e[1] for e in xy[i:j]) / (j - i)
    return sum(abs(e[1] - m) for e in xy[i:j])


def difference_abs(nuage, i, j, k):
    m1 = somme_diff_abs(nuage, i, k)
    m2 = somme_diff_abs(nuage, k, j)
    m = somme_diff_abs(nuage, i, j)
    return abs(m1 + m2 - m)


def optimise_abs(nuage, i, j):
    mx = -1
    ib = None
    for k in range(i + 1, j - 1):
        d = difference_abs(nuage, i, j, k)
        if ib is None or d > mx:
            mx = d
            ib = k
    if ib is None:
        ib = i
        mx = 0
    return ib, mx


%timeit optimise_abs(nuage, 0, len(nuage))

503 ms ± 21.1 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

L’instruction suivante permet de voir où le programme passe la majeure partie de son temps.

[17]:

# %prun optimise_abs(nuage, 0, len(nuage))

La fonction sum cache une boucle, avec la boucle for dans la fonction optimise, cela explique le coût en $O(n^2)$ . Le fait qu’à chaque itération, on passe une observation d’un côté à l’autre de la coupure puis on recalcule les moyennes… Il y a deux façons d’optimiser ce calcul selon qu’on tient compte du fait que les valeurs de $Y$ sont binaires ou non.

Dans le premier cas, il suffit de compter les valeurs 0 ou 1 de part et d’autres de la coupure (histogrammes) pour calculer la moyenne. Lorsque $k$ varie, il suffit de mettre à jour les histogrammes en déduisant et en ajoutant le $Y_k$ aux bons endroits.

[18]:

def histogramme_y(xy, i, j):
    d = [0, 0]
    for x, y in xy[i:j]:
        d[y] += 1
    return d


def somme_diff_histogramme(d):
    m = d[1] * 1.0 / (d[0] + d[1])
    return (1 - m) * d[1] + m * d[0]


def optimise_rapide(nuage, i, j):
    # On calcule les histogrammes.
    d1 = histogramme_y(nuage, i, i + 1)
    d2 = histogramme_y(nuage, i + 1, j)
    d = d1.copy()
    d[0] += d2[0]
    d[1] += d2[1]

    m = somme_diff_histogramme(d)
    m1 = somme_diff_histogramme(d1)
    m2 = somme_diff_histogramme(d2)
    mx = -1
    ib = None
    for k in range(i + 1, j - 1):
        d = abs(m1 + m2 - m)
        if ib is None or d > mx:
            mx = d
            ib = k
        # On met à jour les histogrammes. On ajoute d'un côté, on retranche de l'autre.
        y = nuage[k][1]
        d1[y] += 1
        d2[y] -= 1
        m1 = somme_diff_histogramme(d1)
        m2 = somme_diff_histogramme(d2)
    if ib is None:
        ib = i
        mx = 0
    return ib, mx


# On vérifie qu'on obtient les mêmes résultats.
optimise_rapide(nuage, 0, len(nuage)), optimise_abs(nuage, 0, len(nuage))

[18]:

((565, 235.4096814159292), (565, 235.40968141593424))

C’est carrément plus rapide et cela marche pour toute fonction fct.

[19]:

%timeit optimise_rapide(nuage, 0, len(nuage))

1.63 ms ± 150 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

Si on ne suppose pas que les $Y_i$ sont binaires et qu’ils sont quelconques, les histogrammes contiendront plus de deux éléments. Dans ce cas, il faut conserver deux tableaux triés des $Y_i$ , de part et d’autres de la coupure. Lorsqu’on bouge la coupure $k$ , cela revient à déplacer $Y_k$ d’un tableau à l’autre ce qui se fera par recherche dichotomique donc en $O(\ln n)$ . La mise à jour de la moyenne des valeurs absolues est immédiate si la fonction fct=abs(x-y) et pas forcément immédiate dans le cas général. Lorsque c’est une valeur absolue, il faut utiliser quelques résultats sur la régression quantile.

Q10 - autre nuage de points#

Comment l’algorithme se comporte-t-il lorsque tous les points sont distincts ?

[20]:

import random

X2 = list(range(10))
Y2 = X2

[21]:

import matplotlib.pyplot as plt

plt.plot(X2, Y2, ".")

[21]:

[<matplotlib.lines.Line2D at 0x1f8944d3ba8>]

../../_images/practice_exams_td_note_2018_1_38_1.png

[22]:

nuage2 = [(x, y) for x, y in zip(X2, Y2)]
nuage2.sort()

[23]:

r = recursive(nuage2, 0, len(nuage2), fct)
len(r), r

[23]:

(5, [2, 3, 5, 7, 8])

[24]:

import matplotlib.pyplot as plt

plt.plot(X2, Y2, ".")
for i in r:
    x = nuage2[i][0]
    plt.plot([x, x], [0, 10])

../../_images/practice_exams_td_note_2018_1_41_0.png

La fonction va placer un point dans chaque intervalle, il y aura quasiment autant de points de coupures que de points. Presque autant car l’implémentation a quelques effets de bords comme la boucle for k in range(i+1,j-1): qui ne considère pas les extrémités d’un intervalle comme de potentiels points de coupures.

[25]:

Notebook on github