10. NumPy#
« Soyons clairs : le travail de la science n’a rien à voir avec le consensus. Le consensus est l’affaire de la politique. La science, au contraire, ne requiert qu’un seul chercheur qui se trouve avoir raison, ce qui signifie qu’il ou elle dispose de résultats vérifiables par référence au monde réel. En science, le consensus n’a aucune importance. Ce qui importe, ce sont des résultats reproductibles. » – Michael Crichton
En plus de ce qui est inclus dans Anaconda, ce cours nécessitera les bibliothèques suivantes :
!pip install quantecon
10.1. Vue d’ensemble#
NumPy est une bibliothèque de référence pour le calcul numérique en Python
Largement utilisée dans le monde universitaire, la finance et l’industrie.
Mature, rapide, stable et continuellement améliorée.
Nous avons déjà rencontré du code utilisant NumPy dans les chapitres précédents.
Dans ce chapitre, nous allons étudier de manière plus systématique
les tableaux NumPy et
les opérations fondamentales de traitement de tableaux fournies par NumPy.
(Pour une référence alternative, consultez la documentation officielle de NumPy.)
Nous utiliserons les instructions d’importation suivantes.
import numpy as np
import random
import quantecon as qe
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.axes3d import Axes3D
from matplotlib import cm
10.2. Tableaux NumPy#
L’un des principaux objectifs de NumPy est de permettre un traitement rapide des tableaux.
La structure la plus importante que NumPy définit est un type de données tableau, formellement appelé numpy.ndarray.
Les tableaux NumPy constituent la base d’une grande partie de l’écosystème scientifique Python.
10.2.1. Notions de base#
Pour créer un tableau NumPy ne contenant que des zéros, nous utilisons np.zeros
a = np.zeros(3)
a
array([0., 0., 0.])
type(a)
numpy.ndarray
Les tableaux NumPy ressemblent à certains égards aux listes Python, à deux différences importantes près :
Les données doivent être homogènes : tous les éléments doivent avoir le même type.
Ces types doivent être l’un des types de données (
dtypes) fournis par NumPy.
Parmi les types de données les plus couramment utilisés, on trouve :
float64 : nombre à virgule flottante sur 64 bits
int64 : entier sur 64 bits
bool : True ou False sur 8 bits
Il existe également des dtypes pour représenter les nombres complexes, les entiers non signés, etc.
Sur les machines modernes, le dtype par défaut pour les tableaux est float64
a = np.zeros(3)
type(a[0])
numpy.float64
Si nous voulons utiliser des entiers, nous pouvons le spécifier comme suit :
a = np.zeros(3, dtype=int)
type(a[0])
numpy.int64
10.2.2. Forme et dimension#
Considérons l’affectation suivante
z = np.zeros(10)
Ici, z est un tableau unidimensionnel : ce n’est ni un vecteur ligne ni un vecteur colonne.
z.shape
(10,)
Ici, le tuple représentant la forme du tableau ne contient qu’un seul élément : sa longueur (les tuples à un seul élément se terminent par une virgule).
Pour transformer ce tableau en tableau bidimensionnel, nous pouvons modifier son attribut shape
z.shape = (10, 1) # Convertit le tableau plat en vecteur colonne (bidimensionnel)
z
array([[0.],
[0.],
[0.],
[0.],
[0.],
[0.],
[0.],
[0.],
[0.],
[0.]])
z = np.zeros(4) # Tableau plat
z.shape = (2, 2) # Tableau bidimensionnel
z
array([[0., 0.],
[0., 0.]])
Dans le dernier exemple, nous aurions également pu créer directement le tableau \(2 \times 2\) en passant un tuple à la fonction zeros(), comme
dans z = np.zeros((2, 2)).
10.2.3. Création de tableaux#
Comme nous l’avons vu, la fonction np.zeros crée un tableau de zéros.
Vous pouvez probablement deviner ce que crée np.ones.
La fonction np.empty crée un tableau sans initialiser ses valeurs, afin qu’il puisse être rempli ultérieurement
z = np.empty(3)
z
array([0., 0., 0.])
Les nombres que vous voyez ici sont des valeurs parasites.
(Python alloue 3 morceaux contigus de mémoire de 64 bits, et le contenu existant de ces emplacements mémoire est interprété comme des valeurs float64)
Pour générer une suite de nombres régulièrement espacés sur un intervalle, utilisez np.linspace
z = np.linspace(2, 4, 5) # De 2 à 4, avec 5 éléments
Pour créer une matrice identité, utilisez soit np.identity, soit np.eye
z = np.identity(2)
z
array([[1., 0.],
[0., 1.]])
De plus, les tableaux NumPy peuvent être créés à partir de listes Python, de tuples, etc. en utilisant np.array
z = np.array([10, 20]) # ndarray à partir d'une liste Python
z
array([10, 20])
type(z)
numpy.ndarray
z = np.array((10, 20), dtype=float) # Ici 'float' est équivalent à 'np.float64'
z
array([10., 20.])
z = np.array([[1, 2], [3, 4]]) # Tableau 2D à partir d'une liste de listes
z
array([[1, 2],
[3, 4]])
La fonction np.asarray joue un rôle similaire, mais ne fait pas
de copie distincte des données déjà présentes dans un tableau NumPy.
Pour lire les données d’un tableau à partir d’un fichier texte contenant des données numériques, utilisez np.loadtxt — voir la documentation pour plus de détails.
10.2.4. Indexation de tableaux#
Pour un tableau unidimensionnel, l’indexation fonctionne comme pour les autres séquences Python :
z = np.linspace(1, 2, 5)
z
array([1. , 1.25, 1.5 , 1.75, 2. ])
z[0]
np.float64(1.0)
z[0:2] # Deux éléments, en commençant par l'élément 0
array([1. , 1.25])
z[-1]
np.float64(2.0)
Pour les tableaux 2D, la syntaxe d’indexation est la suivante :
z = np.array([[1, 2], [3, 4]])
z
array([[1, 2],
[3, 4]])
z[0, 0]
np.int64(1)
z[0, 1]
np.int64(2)
Et ainsi de suite.
On peut extraire les lignes et les colonnes de la manière suivante
z[0, :]
array([1, 2])
z[:, 1]
array([2, 4])
Un tableau d’indices entiers peut également servir à sélectionner plusieurs éléments
z = np.linspace(2, 4, 5)
z
array([2. , 2.5, 3. , 3.5, 4. ])
indices = np.array((0, 2, 3))
z[indices]
array([2. , 3. , 3.5])
Enfin, un tableau booléen peut être utilisé comme masque pour sélectionner des éléments
z
array([2. , 2.5, 3. , 3.5, 4. ])
d = np.array([0, 1, 1, 0, 0], dtype=bool)
d
array([False, True, True, False, False])
z[d]
array([2.5, 3. ])
Nous verrons ci-dessous pourquoi cela est utile.
Remarquons également que l’on peut attribuer une même valeur à tous les éléments d’un tableau à l’aide d’une tranche
z = np.empty(3)
z
array([2. , 3. , 3.5])
z[:] = 42
z
array([42., 42., 42.])
10.2.5. Méthodes des tableaux#
Les tableaux NumPy disposent de nombreuses méthodes utiles et optimisées
a = np.array((4, 3, 2, 1))
a
array([4, 3, 2, 1])
a.sort() # Trie a sur place
a
array([1, 2, 3, 4])
a.sum() # Somme
np.int64(10)
a.mean() # Moyenne
np.float64(2.5)
a.max() # Max
np.int64(4)
a.argmax() # Renvoie l'indice de l'élément maximal
np.int64(3)
a.cumsum() # Somme cumulée des éléments de a
array([ 1, 3, 6, 10])
a.cumprod() # Produit cumulé des éléments de a
array([ 1, 2, 6, 24])
a.var() # Variance
np.float64(1.25)
a.std() # Écart-type
np.float64(1.118033988749895)
a.shape = (2, 2)
a.T # Équivalent à a.transpose()
array([[1, 3],
[2, 4]])
Une autre méthode utile à connaître est searchsorted().
Si z est un tableau non décroissant, alors z.searchsorted(a) renvoie l’indice du
premier élément de z qui est >= a
z = np.linspace(2, 4, 5)
z
array([2. , 2.5, 3. , 3.5, 4. ])
z.searchsorted(2.2)
np.int64(1)
10.3. Opérations arithmétiques#
Les opérateurs +, -, *, / et ** agissent tous élément par élément sur les tableaux
a = np.array([1, 2, 3, 4])
b = np.array([5, 6, 7, 8])
a + b
array([ 6, 8, 10, 12])
a * b
array([ 5, 12, 21, 32])
Nous pouvons ajouter un scalaire à chaque élément comme suit
a + 10
array([11, 12, 13, 14])
La multiplication par un scalaire fonctionne de la même manière
a * 10
array([10, 20, 30, 40])
Les tableaux bidimensionnels suivent les mêmes règles générales
A = np.ones((2, 2))
B = np.ones((2, 2))
A + B
array([[2., 2.],
[2., 2.]])
A + 10
array([[11., 11.],
[11., 11.]])
A * B
array([[1., 1.],
[1., 1.]])
En particulier, A * B n’est pas le produit matriciel, c’est un produit élément par élément.
10.4. Multiplication matricielle#
Nous utilisons le symbole @ pour la multiplication matricielle, comme suit :
A = np.ones((2, 2))
B = np.ones((2, 2))
A @ B
array([[2., 2.],
[2., 2.]])
La syntaxe fonctionne avec des tableaux plats — NumPy devine intelligemment ce que vous voulez :
A @ (0, 1)
array([1., 1.])
Comme nous effectuons une post-multiplication, le tuple est traité comme un vecteur colonne.
10.5. Broadcasting#
(Cette section prolonge une excellente discussion du broadcasting fournie par Jake VanderPlas.)
Note
Le broadcasting est un aspect très important de NumPy. Cependant, les formes avancées de broadcasting sont relativement complexes et certains des détails ci-dessous peuvent être survolés en première lecture.
Lors d’une opération élément par élément, les tableaux n’ont pas nécessairement la même forme.
Lorsque cela se produit, NumPy étendra automatiquement les tableaux à la même forme chaque fois que possible.
Cette fonctionnalité utile (mais parfois déroutante) de NumPy s’appelle broadcasting.
L’intérêt du broadcasting est que
les boucles
forpeuvent être évitées, ce qui aide le code numérique à s’exécuter rapidement etle broadcasting peut nous permettre de mettre en œuvre des opérations sur des tableaux sans réellement créer certaines dimensions de ces tableaux en mémoire, ce qui peut être important lorsque les tableaux sont grands.
Par exemple, supposons que a soit un tableau \(3 \times 3\) (a -> (3, 3)), tandis que b est un tableau plat de trois éléments (b -> (3,)).
Lors de leur addition, NumPy étendra automatiquement b -> (3,) en b -> (3, 3).
L’addition élément par élément produira un tableau \(3 \times 3\)
a = np.array(
[[1, 2, 3],
[4, 5, 6],
[7, 8, 9]])
b = np.array([3, 6, 9])
a + b
array([[ 4, 8, 12],
[ 7, 11, 15],
[10, 14, 18]])
Voici une représentation visuelle de cette opération de broadcasting :
Qu’en est-il de b -> (3, 1) ?
Dans ce cas, NumPy étendra automatiquement b -> (3, 1) en b -> (3, 3).
L’addition élément par élément produira alors une matrice \(3 \times 3\)
b.shape = (3, 1)
a + b
array([[ 4, 5, 6],
[10, 11, 12],
[16, 17, 18]])
Voici une représentation visuelle de cette opération de broadcasting :
Dans certains cas, les deux opérandes seront étendus.
Lorsque nous avons a -> (3,) et b -> (3, 1), a sera étendu en a -> (3, 3), et b sera étendu en b -> (3, 3).
Dans ce cas, l’addition élément par élément produira une matrice \(3 \times 3\)
a = np.array([3, 6, 9])
b = np.array([2, 3, 4])
b.shape = (3, 1)
a + b
array([[ 5, 8, 11],
[ 6, 9, 12],
[ 7, 10, 13]])
Voici une représentation visuelle de cette opération de broadcasting :
Bien que le broadcasting soit très utile, il peut parfois sembler déroutant.
Par exemple, essayons d’additionner a -> (3, 2) et b -> (3,).
a = np.array(
[[1, 2],
[4, 5],
[7, 8]])
b = np.array([3, 6, 9])
a + b
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[62], line 7
3 [4, 5],
4 [7, 8]])
5 b = np.array([3, 6, 9])
6
----> 7 a + b
ValueError: operands could not be broadcast together with shapes (3,2) (3,)
La ValueError nous indique que les opérandes n’ont pas pu être broadcastées ensemble.
La représentation suivante montre pourquoi cette opération de broadcasting est impossible :
Nous pouvons voir que NumPy ne peut pas étendre les tableaux à la même taille.
C’est parce que, lorsque b est étendu de b -> (3,) à b -> (3, 3), NumPy ne peut pas faire correspondre b avec a -> (3, 2).
Les choses deviennent encore plus délicates lorsque nous passons à des dimensions supérieures.
Pour nous aider, nous pouvons utiliser la liste de règles suivante :
Étape 1 : Lorsque les dimensions de deux tableaux ne correspondent pas, NumPy étend celui qui a le moins de dimensions en ajoutant une ou plusieurs dimensions à gauche des dimensions existantes.
Par exemple, si
a -> (3, 3)etb -> (3,), alors le broadcasting ajoutera une dimension à gauche de sorte queb -> (1, 3);Si
a -> (2, 2, 2)etb -> (2, 2), alors le broadcasting ajoutera une dimension à gauche de sorte queb -> (1, 2, 2);Si
a -> (3, 2, 2)etb -> (2,), alors le broadcasting ajoutera deux dimensions à gauche de sorte queb -> (1, 1, 2)(on peut aussi voir ce processus comme un passage par l”Étape 1 deux fois).
Étape 2 : Lorsque les deux tableaux ont la même dimension mais des formes différentes, NumPy tentera d’étendre les dimensions où l’indice de forme est 1.
Par exemple, si
a -> (1, 3)etb -> (3, 1), alors le broadcasting étendra les dimensions de forme 1 dansaetbde sorte quea -> (3, 3)etb -> (3, 3);Si
a -> (2, 2, 2)etb -> (1, 2, 2), alors le broadcasting étendra la première dimension debde sorte queb -> (2, 2, 2);Si
a -> (3, 2, 2)etb -> (1, 1, 2), alors le broadcasting étendrabsur toutes les dimensions de forme 1 de sorte queb -> (3, 2, 2).
Étape 3 : Après les étapes 1 et 2, si les deux tableaux ne correspondent toujours pas, une
ValueErrorsera levée. Par exemple, supposonsa -> (2, 2, 3)etb -> (2, 2)Par l”Étape 1,
bsera étendu enb -> (1, 2, 2);Par l”Étape 2,
bsera étendu enb -> (2, 2, 2);Nous pouvons voir qu’ils ne correspondent pas l’un à l’autre après les deux premières étapes. Ainsi, une
ValueErrorsera levée
10.6. Mutabilité et copie de tableaux#
Les tableaux NumPy sont des types de données mutables, comme les listes Python.
En d’autres termes, leur contenu peut être modifié en mémoire après l’initialisation.
C’est pratique mais, combiné au modèle de nommage et de référence de Python, cela peut être source d’erreurs pour les personnes qui découvrent NumPy.
Dans cette section, nous passons en revue quelques points clés.
10.6.1. Mutabilité#
Nous avons déjà vu des exemples de mutabilité ci-dessus.
Voici un autre exemple de mutation d’un tableau NumPy
a = np.array([42, 44])
a
array([42, 44])
a[-1] = 0 # Change le dernier élément en 0
a
array([42, 0])
La mutabilité entraîne le comportement suivant, qui peut surprendre les utilisateurs de MATLAB.
rng = np.random.default_rng()
a = rng.standard_normal(3)
a
array([ 0.58571122, -0.41238965, 0.13644795])
b = a
b[0] = 0.0
a
array([ 0. , -0.41238965, 0.13644795])
Ce qui s’est passé, c’est que nous avons modifié a en modifiant b.
Le nom b est lié à a et devient simplement une autre référence au
tableau (le modèle d’affectation de Python est décrit plus en détail plus loin dans le cours).
Toute modification effectuée par l’intermédiaire de b affecte donc également le tableau accessible par a.
Ce comportement est généralement préférable, car il évite les copies inutiles.
Cela signifie que nous ne faisons circuler que des pointeurs vers les données, plutôt que de faire des copies.
Faire des copies est coûteux à la fois en termes de vitesse et de mémoire.
10.6.2. Faire des copies#
Il est bien sûr possible de faire de b une copie indépendante de a lorsque cela est nécessaire.
Cela peut être fait en utilisant np.copy
a = rng.standard_normal(3)
a
array([-1.35105448, 0.48191099, 0.13628037])
b = np.copy(a)
b
array([-1.35105448, 0.48191099, 0.13628037])
Maintenant b est une copie indépendante (appelée copie profonde)
b[:] = 1
b
array([1., 1., 1.])
a
array([-1.35105448, 0.48191099, 0.13628037])
Notez que la modification de b n’a pas affecté a.
10.7. Fonctionnalités supplémentaires#
Examinons quelques autres fonctionnalités utiles de NumPy.
10.7.1. Fonctions universelles#
NumPy fournit des versions des fonctions standard log, exp, sin, etc. qui agissent élément par élément sur les tableaux
z = np.array([1, 2, 3])
np.sin(z)
array([0.84147098, 0.90929743, 0.14112001])
Cela évite d’écrire explicitement une boucle parcourant les éléments un à un, comme dans l’exemple suivant
n = len(z)
y = np.empty(n)
for i in range(n):
y[i] = np.sin(z[i])
Parce qu’elles agissent élément par élément sur les tableaux, ces fonctions sont parfois appelées fonctions vectorisées.
Dans le jargon de NumPy, elles sont aussi appelées ufuncs, ou fonctions universelles.
Comme nous l’avons vu ci-dessus, les opérations arithmétiques usuelles (+, *, etc.) fonctionnent
également élément par élément. En combinant ces opérations aux ufuncs, on dispose d’un vaste ensemble d’opérations vectorisées et efficaces.
z
array([1, 2, 3])
(1 / np.sqrt(2 * np.pi)) * np.exp(- 0.5 * z**2)
array([0.24197072, 0.05399097, 0.00443185])
Les fonctions définies par l’utilisateur ne sont pas nécessairement vectorisées.
Par exemple, passer un tableau NumPy à la fonction f définie ci-dessous provoque une ValueError
def f(x):
return 1 if x > 0 else 0
La fonction NumPy np.where fournit une alternative vectorisée :
x = rng.standard_normal(4)
x
array([ 0.67312016, 1.38797249, -1.08844898, -0.68916089])
np.where(x > 0, 1, 0) # Insère 1 si x > 0 est vrai, sinon 0
array([1, 1, 0, 0])
Vous pouvez aussi utiliser np.vectorize pour vectoriser une fonction donnée
f = np.vectorize(f)
f(x) # Passe le même vecteur x que dans l'exemple précédent
array([1, 1, 0, 0])
Cependant, cette approche n’atteint pas toujours la même vitesse qu’une fonction vectorisée plus soigneusement conçue.
(Plus tard, nous verrons que JAX dispose d’une version puissante de np.vectorize qui peut générer, et génère généralement, du code hautement efficace.)
10.7.2. Comparaisons#
Les opérateurs de comparaison appliqués à des tableaux NumPy agissent généralement élément par élément
z = np.array([2, 3])
y = np.array([2, 3])
z == y
array([ True, True])
y[0] = 5
z == y
array([False, True])
z != y
array([ True, False])
La situation est similaire pour >, <, >= et <=.
Nous pouvons aussi effectuer des comparaisons avec des scalaires
z = np.linspace(0, 10, 5)
z
array([ 0. , 2.5, 5. , 7.5, 10. ])
z > 3
array([False, False, True, True, True])
Ceci est particulièrement utile pour l”extraction conditionnelle
b = z > 3
b
array([False, False, True, True, True])
z[b]
array([ 5. , 7.5, 10. ])
En pratique, cette sélection est souvent réalisée en une seule étape
z[z > 3]
array([ 5. , 7.5, 10. ])
10.7.3. Sous-packages#
NumPy fournit des fonctionnalités supplémentaires liées à la programmation scientifique via ses sous-packages.
Nous avons déjà vu comment nous pouvons générer des variables aléatoires en utilisant le
Generator aléatoire de NumPy.
z = rng.standard_normal(10000) # Génère des variables normales centrées réduites
y = rng.binomial(10, 0.5, size=1000) # 1 000 tirages de Bin(10, 0.5)
y.mean()
np.float64(4.929)
Un autre sous-package couramment utilisé est np.linalg
A = np.array([[1, 2], [3, 4]])
np.linalg.det(A) # Calcule le déterminant
np.float64(-2.0000000000000004)
np.linalg.inv(A) # Calcule l'inverse
array([[-2. , 1. ],
[ 1.5, -0.5]])
Une grande partie de ces fonctionnalités est aussi disponible dans SciPy, une collection de modules reposant sur NumPy.
Nous couvrirons les versions SciPy plus en détail bientôt.
Pour une liste complète de ce qui est disponible dans NumPy, consultez cette documentation.
10.7.4. Multithreading implicite#
Précédemment, nous avons abordé le concept de parallélisation via le multithreading.
NumPy tente d’implémenter le multithreading dans une grande partie de son code compilé.
Observons ce comportement à l’aide d’un exemple.
Le morceau de code suivant calcule les valeurs propres d’un grand nombre de matrices générées aléatoirement.
Son exécution prend quelques secondes.
n = 20
m = 1000
for i in range(n):
X = rng.standard_normal((m, m))
λ = np.linalg.eigvals(X)
Maintenant, examinons la sortie du moniteur système htop sur notre machine pendant que ce code s’exécute :
Nous pouvons voir que 4 des 8 processeurs fonctionnent à pleine vitesse.
Cela est dû au fait que la routine eigvals de NumPy divise proprement les tâches et
les distribue à différents threads.
10.8. Exercices#
Exercice 10.1
Considérons l’expression polynomiale
Plus tôt, vous avez écrit une fonction simple p(x, coeff) pour évaluer (10.1) sans chercher à optimiser les performances.
Écrivez maintenant une nouvelle fonction qui effectue le même travail, mais qui utilise des tableaux NumPy et des opérations sur tableaux pour ses calculs, sans utiliser de boucle Python explicite.
(Une telle fonctionnalité est déjà implémentée sous la forme de np.poly1d, mais pour les besoins de l’exercice, n’utilisez pas cette classe)
Indication
Utilisez np.cumprod()
Solution
La fonction suivante répond à la question
def p(x, coef):
X = np.ones_like(coef)
X[1:] = x
y = np.cumprod(X) # y = [1, x, x**2,...]
return coef @ y
Testons-le
x = 2
coef = np.linspace(2, 4, 3)
print(coef)
print(p(x, coef))
# À des fins de comparaison
q = np.poly1d(np.flip(coef))
print(q(x))
[2. 3. 4.]
24.0
24.0
Exercice 10.2
Soit q un tableau NumPy de longueur n avec q.sum() == 1.
Supposons que q représente une fonction de masse de probabilité.
Nous souhaitons générer une variable aléatoire discrète \(x\) telle que \(\mathbb P\{x = i\} = q_i\).
En d’autres termes, x prend des valeurs dans range(len(q)) et x = i avec probabilité q[i].
L’algorithme standard (de transformation inverse) est le suivant :
Divisez l’intervalle unité \([0, 1]\) en \(n\) sous-intervalles \(I_0, I_1, \ldots, I_{n-1}\) tels que la longueur de \(I_i\) soit \(q_i\).
Effectuez un tirage \(U\) selon la loi uniforme sur \([0, 1]\) et renvoyez le \(i\) tel que \(U \in I_i\).
La probabilité d’obtenir la valeur \(i\) est égale à la longueur de l’intervalle \(I_i\), qui est égale à \(q_i\).
Nous pouvons implémenter l’algorithme comme suit
from random import uniform
def sample(q):
a = 0.0
U = uniform(0, 1)
for i in range(len(q)):
if a < U <= a + q[i]:
return i
a = a + q[i]
Si vous ne voyez pas comment cela fonctionne, essayez de suivre pas à pas l’exécution de l’algorithme sur un exemple simple, tel que q = [0.25, 0.75]
Il est utile de dessiner les intervalles sur papier.
L’objectif de l’exercice est d’en proposer une version plus rapide à l’aide de NumPy, en évitant les boucles explicites
Indication
Utilisez np.searchsorted et np.cumsum
Si vous le pouvez, implémentez la fonctionnalité sous la forme d’une classe appelée DiscreteRV, où
les données pour une instance de la classe sont le vecteur de probabilités
qla classe dispose d’une méthode
draw(), qui renvoie un tirage selon l’algorithme décrit ci-dessus
Si vous le pouvez, écrivez la méthode de sorte que draw(k) renvoie k tirages de q.
Solution
Voici notre première tentative de solution :
from numpy import cumsum
class DiscreteRV:
"""
Generates an array of draws from a discrete random variable with vector of
probabilities given by q.
"""
def __init__(self, q, seed=None):
"""
The argument q is a NumPy array, or array like, nonnegative and sums
to 1
"""
self.q = q
self.Q = cumsum(q)
self.rng = np.random.default_rng(seed)
def draw(self, k=1):
"""
Returns k draws from q. For each such draw, the value i is returned
with probability q[i].
"""
return self.Q.searchsorted(self.rng.uniform(0, 1, size=k))
La logique peut sembler difficile au premier abord. Prenez le temps d’examiner chaque étape du code.
Il y a cependant un problème ici.
Supposons que q soit modifié après la création d’une instance de discreteRV, par exemple par
q = (0.1, 0.9)
d = DiscreteRV(q)
d.q = (0.5, 0.5)
Le problème est que Q ne change pas en conséquence, et Q correspond aux
données utilisées dans la méthode draw.
Une première solution consisterait à recalculer Q à chaque appel de la méthode draw.
Mais c’est inefficace par rapport à un calcul unique de Q.
Une solution plus efficace consiste à utiliser un descripteur Python.
Une solution de la bibliothèque quantecon utilisant des descripteurs et se comportant comme nous le souhaitons peut être trouvée ici.
Exercice 10.3
Rappelez-vous notre discussion précédente sur la fonction de répartition empirique.
L’objectif de cet exercice est de
Rendre la méthode
__call__plus efficace en utilisant NumPy.Ajouter une méthode qui trace la FDE sur \([a, b]\), où \(a\) et \(b\) sont des paramètres de la méthode.
Solution
Un exemple de solution est donné ci-dessous.
La solution ci-dessous adapte ce code de QuantEcon en lui ajoutant une méthode de tracé
"""
Modifies ecdf.py from QuantEcon to add in a plot method
"""
class ECDF:
"""
One-dimensional empirical distribution function given a vector of
observations.
Parameters
----------
observations : array_like
An array of observations
Attributes
----------
observations : array_like
An array of observations
"""
def __init__(self, observations):
self.observations = np.asarray(observations)
def __call__(self, x):
"""
Evaluates the ecdf at x
Parameters
----------
x : scalar(float)
The x at which the ecdf is evaluated
Returns
-------
scalar(float)
Fraction of the sample less than x
"""
return np.mean(self.observations <= x)
def plot(self, ax, a=None, b=None):
"""
Plot the ecdf on the interval [a, b].
Parameters
----------
a : scalar(float), optional(default=None)
Lower endpoint of the plot interval
b : scalar(float), optional(default=None)
Upper endpoint of the plot interval
"""
# === choisit un intervalle raisonnable si [a, b] n'est pas spécifié === #
if a is None:
a = self.observations.min() - self.observations.std()
if b is None:
b = self.observations.max() + self.observations.std()
# === génère le tracé === #
x_vals = np.linspace(a, b, num=100)
f = np.vectorize(self.__call__)
ax.plot(x_vals, f(x_vals))
plt.show()
Voici un exemple d’utilisation
fig, ax = plt.subplots()
rng = np.random.default_rng()
X = rng.standard_normal(1000)
F = ECDF(X)
F.plot(ax)
Exercice 10.4
Rappelez-vous que le broadcasting dans NumPy peut nous aider à effectuer des opérations élément par élément sur des tableaux ayant un nombre de dimensions différent sans utiliser de boucles for.
Dans cet exercice, essayez d’utiliser des boucles for pour reproduire le résultat des opérations de broadcasting suivantes.
Partie 1 : Essayez de reproduire cet exemple simple en utilisant des boucles for et comparez vos résultats avec l’opération de broadcasting ci-dessous.
rng = np.random.default_rng(123)
x = rng.standard_normal((4, 4))
y = rng.standard_normal(4)
A = x / y
Voici la sortie
print(A)
Partie 2 : Passez à la reproduction du résultat de l’opération de broadcasting suivante. Parallèlement, comparez les vitesses du broadcasting et de la boucle for que vous implémentez.
Pour cette partie de l’exercice, vous pouvez utiliser les fonctions tic/toc de la bibliothèque quantecon pour chronométrer l’exécution.
Assurons-nous que cette bibliothèque est installée.
!pip install quantecon
Nous pouvons maintenant importer le package quantecon.
rng = np.random.default_rng(123)
x = rng.standard_normal((1000, 100, 100))
y = rng.standard_normal(100)
with qe.Timer("Broadcasting operation"):
B = x / y
Broadcasting operation: 0.0113 seconds elapsed
Voici la sortie
print(B)
Solution
Solution de la Partie 1
rng = np.random.default_rng(123)
x = rng.standard_normal((4, 4))
y = rng.standard_normal(4)
C = np.empty_like(x)
n = len(x)
for i in range(n):
for j in range(n):
C[i, j] = x[i, j] / y[j]
Comparez les résultats pour vérifier votre réponse
print(C)
Vous pouvez aussi utiliser array_equal() pour vérifier votre réponse
print(np.array_equal(A, C))
True
Solution de la Partie 2
rng = np.random.default_rng(123)
x = rng.standard_normal((1000, 100, 100))
y = rng.standard_normal(100)
with qe.Timer("For loop operation"):
D = np.empty_like(x)
d1, d2, d3 = x.shape
for i in range(d1):
for j in range(d2):
for k in range(d3):
D[i, j, k] = x[i, j, k] / y[k]
For loop operation: 3.8515 seconds elapsed
Notez que la boucle for prend beaucoup plus de temps que l’opération de broadcasting.
Comparez les résultats pour vérifier votre réponse
print(D)
print(np.array_equal(B, D))
True