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 premier ordre pour la programmation numérique
Largement utilisée dans le monde universitaire, la finance et l’industrie.
Mature, rapide, stable et en développement continu.
Nous avons déjà vu du code impliquant NumPy dans les cours précédents.
Dans ce cours, nous allons entamer une discussion plus systématique de
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 importations 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#
Le problème essentiel que NumPy résout est le traitement rapide de tableaux.
La structure la plus importante que NumPy définit est un type de données tableau, formellement appelé numpy.ndarray.
Les tableaux NumPy alimentent une très grande proportion de l’écosystème scientifique de 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 quelque peu aux listes natives de Python, sauf que
Les données doivent être homogènes (tous les éléments du même type).
Ces types doivent être l’un des types de données (
dtypes) fournis par NumPy.
Les plus importants de ces dtypes sont :
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 plat — ni vecteur ligne ni vecteur colonne.
z.shape
(10,)
Ici, le tuple de forme n’a qu’un seul élément, qui est la longueur du tableau (les tuples à un seul élément se terminent par une virgule).
Pour lui donner une dimension supplémentaire, nous pouvons modifier l’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 cas, pour créer le tableau 2x2, nous pourrions aussi passer 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.
Une fonction apparentée est np.empty, qui crée des tableaux en mémoire qui peuvent ensuite être remplis de données
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 établir une grille de nombres régulièrement espacés, 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]])
Voir aussi np.asarray, qui remplit une fonction 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 plat, l’indexation est la même que pour les 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.
Les colonnes et les lignes peuvent être extraites comme suit
z[0, :]
array([1, 2])
z[:, 1]
array([2, 4])
Les tableaux NumPy d’entiers peuvent aussi être utilisés pour extraire des é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 de dtype bool peut être utilisé pour extraire 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.
Une remarque en passant : tous les éléments d’un tableau peuvent être fixés à un même nombre en utilisant la notation de 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 disposent de méthodes utiles, toutes soigneusement 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 qui vaut la peine d’être connue 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 scalaire est similaire
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. En même temps, le broadcasting avancé est relativement complexe et certains des détails ci-dessous peuvent être survolés en première lecture.
Dans les opérations élément par élément, les tableaux peuvent ne pas avoir 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.
Voici une représentation visuelle pour montrer pourquoi ce broadcasting ne peut pas être exécuté :
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é (muté) 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 conduire à des erreurs chez les débutants en 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é conduit au comportement suivant (qui peut être choquant pour les programmeurs MATLAB…)
rng = np.random.default_rng()
a = rng.standard_normal(3)
a
array([ 0.82706403, 1.36260144, -1.67770637])
b = a
b[0] = 0.0
a
array([ 0. , 1.36260144, -1.67770637])
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).
Il a donc des droits égaux pour apporter des modifications à ce tableau.
C’est en fait le comportement par défaut le plus sensé !
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.99303188, 0.91306177, -0.58199588])
b = np.copy(a)
b
array([-1.99303188, 0.91306177, -0.58199588])
Maintenant b est une copie indépendante (appelée copie profonde)
b[:] = 1
b
array([1., 1., 1.])
a
array([-1.99303188, 0.91306177, -0.58199588])
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 élimine le besoin de boucles explicites élément par élément telles que
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, et combiner celles-ci avec les ufuncs donne un très large ensemble de fonctions rapides opérant élément par élément.
z
array([1, 2, 3])
(1 / np.sqrt(2 * np.pi)) * np.exp(- 0.5 * z**2)
array([0.24197072, 0.05399097, 0.00443185])
Toutes les fonctions définies par l’utilisateur n’agissent pas élément par élément.
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([-1.04361928, -1.46393278, -0.16153239, -0.16967398])
np.where(x > 0, 1, 0) # Insère 1 si x > 0 est vrai, sinon 0
array([0, 0, 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([0, 0, 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#
En règle générale, les comparaisons sur les tableaux sont effectuées é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. ])
Bien sûr, nous pouvons — et faisons souvent — effectuer ceci 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.982)
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 qui sont construits par-dessus 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é.
Examinons un exemple pour voir cela en action.
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 tenir compte de l’efficacité.
É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, plutôt qu’une quelconque forme de boucle Python.
(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
Ce code fait le travail
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\).
Tirez une variable aléatoire uniforme \(U\) sur \([0, 1]\) et renvoyez le \(i\) tel que \(U \in I_i\).
La probabilité de tirer \(i\) est la longueur de \(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 le flux pour un exemple simple, tel que q = [0.25, 0.75]
Il est utile de dessiner les intervalles sur papier.
Votre exercice consiste à l’accélérer en utilisant 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 n’est pas évidente, mais si vous prenez votre temps et la lisez lentement, vous comprendrez.
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.
Pour gérer cela, une option consiste à calculer Q chaque fois que la méthode draw
est appelée.
Mais c’est inefficace par rapport à un calcul unique de Q.
Une meilleure option consiste à utiliser des descripteurs.
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.
Votre tâche consiste à
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.
En substance, nous avons simplement repris ce code de QuantEcon et ajouté 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.0109 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.7256 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