Introduction à Python

(Master 2 Mathématiques, modélisation et apprentissage)

  • Première version date de 1991 (Guido van Rossum, Pays-Bas).
  • On utilisera Python 3 dans ce cours. Python 2 ne sera plus maintenu à partir de janvier 2020.
  • Langage interprété, avec des usages très variés (calcul scientifique, web, interface graphique,...)
  • **Open Source**, en très forte croissance depuis quelques années, et langage le plus utilisé par les développeurs aujourd'hui.
  • Communauté d’utilisateurs très active (StackOverFlow.com)
  • Quelques softwares écrits en Python : BitTorrent, Dropbox...

Langage interprété vs compilé

  • Exemples de langages interprétés : Python, Matlab, Scilab, Octave, R
  • Exemples de langages compilés : C, C++, Java
  • Vitesse d’exécution : interprété < compilé

    Pourquoi, dans ce cas, considérer Python pour du calcul scientifique ?

  • Temps d’implémentation vs temps d’exécution : langage lisible et épuré ⇒ développement et maintenance rapides
  • Exécution en Python : rapide si les passages critiques sont exécutés avec un langage compilé : de nombreuses fonctions sont compilées et le code est interfaçable avec C/C++/FORTRAN

Fonctionnement

  1. Ouverture de l’environnement de développement (ex : Spyder, Jupyter Notebook). Non indispensable mais conseillé.
  2. Au choix :
    • Commande en ligne
    • Ecriture d’un script → exécution du script
    • Ecriture d’une fonction → chargement de la fonction → appel à la fonction

Tutoriel/aide dans la console : **nom_fonction?**

Les types et opérations sous Python

Voyons maintenant les principaux types et opérations sous Python.

Booléens

Deux valeurs possibles : False, True

In [ ]:
type(False)
In [ ]:
x = (1<2)
x=3
type(x)

Opérateurs de comparaison : ==, !=, >, >=, <, <=

In [ ]:
2 <=8 <15

Opérateurs logiques : not, or, and

In [ ]:
(3 == 3) or (9 > 24) 
In [ ]:
(9 > 24) and (3 == 3)

int, float, complex

In [ ]:
type(2**100)  
type(3.6)
type(3+2j)
In [ ]:
2*3      # produit
#2**3     # puissance
#20/3    # division flottante
#20//3   # division entière
#20%3    # modulo
#(9+5j).real
#(9+5j).imag
#abs(3+4j) 

les chaîne de caractères

In [ ]:
type('abc') 
c1 = 'L’eau vive'
c2 = ' est "froide" !'
c1+c2 # concatenation
#c1*2 #repetition 
#c1[2] 
#c1[-2]
#c1[2:5]
#len("abc") # longueur 
#"abcde".split("c") # scinde 
#"a−ha".replace('a','o') # 'o-ho' 
#'-'.join(['ci', 'joint']) # 'ci-joint'
#'abracadabra'.count('bra') # 2
#'PETIT'.lower() # 'petit'
#'grand'.upper() # 'GRAND' 

les listes

list : collection hétérogène, ordonnée et modifiable d’éléments séparés par des virgules, entourée de crochets.

In [ ]:
my_list=[4,7,3.7,'E',5,7] 
#type(my_list) # list
#len(my_list) # longueur
#my_list[0] # premier terme
#my_list[1:3] # [7, 3.7]
#[0,1] + [2,4] # [0,1,2,4] (concatenation)
#l = [2,4,5,9,1,6,4]
#k = [x for x in l if x<6] #(extract under condition)
#range(6)          # attention différent en Python 2 (liste) et 3 (boucle)
#range(2, 9, 2) 
#[x for x in range(3,6)]

opérations sur les listes

In [ ]:
nombres = [17, 38, 10, 25, 72]
nombres.sort()
#nombres.append(12)
#nombres.reverse()
#nombres.remove(38)
#print(nombres.index(17))
#nombres[0] = 11
#nombres[1:3] = [14, 17, 2]
#nombres.count(17) # 2 
nombres

**ATTENTION**, par défaut, en Python, les listes ne sont pas copiées.

In [ ]:
x = [4, 2, 10, 9, "toto"]
print(x)
y = x              # y: seulement un alias de x, pas une copie
y[2] = "tintin"    # change x(2) en "tintin"
print(x) 
#x = [4, 2, 10, 9, "toto"]
#y = x[:]   # On demande une copie
#y[2] = "tintin" # modifie y mais pas x
#print(x) #  [4, 2, 10, 9, "toto"]
#print(y) # [4, 2, "tintin", 9, "toto"]

les tuples

tuple : collection hétérogène, ordonnée, immuable.

In [ ]:
t=(5,7)
type(t) # tuple
t[0] # 5
#t[0]=2 # error: item assignment for tuple

Les boucles et opérateurs

if – [elif] – [else]

In [ ]:
x = 11
if x < 0:
    print("x est negatif")
elif x % 2:
    print("x est positif et impair")
else:
    print("x n'est pas negatif et est pair")
    print("Eh oui !")

while

In [ ]:
N = 0
x = 687687.567476
while (x > 0):
    x//=2
    N+=1
print("Approx. de log_2(x) : " + str(N-1))

for

In [ ]:
for lettre in "ciao":
    print(lettre)
    
for x in ["\n",2,'a', 3.14,"\n"]: 
    print(x) 

Les fonctions

Exercice

Ecrire une fonction f qui prend $n$ en entrée et calcule $$\sum_{k=1}^{n-1} k^2$$

In [ ]:
def f(n):
    print(sum([k**2 for k in range(n)])) 
In [ ]:
f(3)

Exercice

Ecrire une fonction f qui calcule le nième terme de la suite de Fibonacci initialisée par a et b, dont la relation de récurrence est $$ x_{n+2} = x_{n+1} + x_n.$$

In [ ]:
def fib(n,a=0,b=1): #0,1: val. par defaut     
     '''n-th term of Fibonacci sequence starting at a,b.'''#tutoriel de fib
     for i in range(n):
            z = a+b
            a=b
            b=z
     return a 
In [ ]:
fib(7,2,2)

Les Bibliothèques d’usage courant

  • NumPy : manipulation de tableaux numériques, fonctions mathématiques de base, simulation de variables aléatoires...
  • SciPy : fonctions mathématiques plus avancées (résolution d’équations, d’équations différentielles, calcul d’intégrales...)
  • Matplotlib : visualisation de données sous forme de graphiques scikit-learn : machine learning
  • SymPy : calcul symbolique

Calcul numérique = manipulations de nombres décimaux $\neq$ Calcul symbolique = manipulation d’expressions symboliques

Exemple : racines de $x^2 − x − 1 = 0$

  • calcul symbolique : $\frac{1+ \sqrt{5}}{2}$ , $\frac{1- \sqrt{5}}{2}$
  • calcul numérique : 1.618034, - 0.6180340

Import de bibliothèques ou de fonctions

  • import ma_bibliotheque
  • ma_bibliotheque.la_fonction(...)

  • import ma_bibliotheque as bibli # raccourci

  • bibli.la_fonction(...)

Moins précis (car la bibliothèque d’origine des fonctions n’est pas précisée à leur appel) : from ma_bibliotheque import la_fonction

  • from ma_bibliotheque import la_fonction
  • la_fonction (...)

  • from ma_bibliotheque import ∗

  • la_fonction (...)
  • Attention si on importe plusieurs bibliothèques ayant les mêmes fonctions...

La librairie Numpy

In [ ]:
import numpy as np  # toujours commencer par cette commande quand on veut utiliser la bibliothèque numpy
In [ ]:
b = np.array([[8,3,2,4],[5,1,6,0],[9,7,4,1]])
print(b)
#type(b) # numpy.ndarray       
#b.dtype # datatype: int
b.shape # (3,4)
#c = np.array([[8,2],[5,6],[9,7]], dtype=complex)
#c.dtype # datatype: complex
#c[0,0] # 8+0j

#More than 2 dimensions:
#d = np.array([[[8,3],[1,2]],[[5,1],[4,5]],[[9,7],[4,5]]])
#d.shape # (3, 2, 2)

 ##Reshaping:
#x = d.reshape(4,3) # tableau de taille (4,3)
#d.reshape(12,1) # tableau de taille (12,1)
#d.reshape(12,) # tableau unidimensionel de taille 12
#np.insert(np.arange(4,9),3,17) # 4,5,6,17,7,8

Attention, il existe aussi une classe numpy.matrix mais il est recommandé d'utiliser numpy.array.

Opérations sur les tableaux de nombres numpy

In [ ]:
X = np.arange(start=5,step=3,stop=16) # 5, 8,11,14
#A = np.ones((2,3)) # matrix filled with ones
#B = X.reshape(2,2)
#C = np.zeros((3,2)) # matrix filled with zeros
#D = np.eye(2) # identity matrix
#np.diag([1,2]) # diagonal matrix
#E = C+np.ones(C.shape) # addition: same as C+1
#F = B*D # entry-wise multiplication
#J=np.dot(B,D) # linear algebra product
#G = F.T # transpose matrix
#H = np.exp(G) #as most functions, exp is entry-wise(else use np.vectorize(my_function))
#x = np.array([4, 2, 1, 5, 1, 10])
#y=np.logical_and(x>=3, x<= 9, x!=1) # [T,F,F,T,F,F]
#x[y] 
#print(np.mean(np.random.randn(1000)>1.96)) 

algèbre linéaire avec numpy

In [ ]:
A = np.array([[2, 1, 1], [4, 3, 0]])
B = np.array([[1, 2], [12, 0]])
C = np.array([[1, 2], [12, 0], [-1, 2]])
D = np.array([[1, 2, -4], [2, 0, 0], [1, 2, 3]])
E = np.concatenate((A,B), axis=1)
F = np.concatenate((C,D), axis=1)
G = np.concatenate((E,F),axis=0)
H = np.random.randn(5,5)
I = H*G                       # produit terme à terme 
B5 = B**5                     # puissance terme à terme
B5 =np.linalg.matrix_power(B, 5)
Bm1 = np.linalg.inv(B)
dB = np.linalg.det(B)
x = np.linalg.solve(B,[3,12]) #résout B*x=[[3],[12]]
print(E)

analyse spectrale avec numpy

In [ ]:
A = np.array([[1, 2], [12, 3]])
x = np.linalg.eigvals(A) # eigenvalues
# eigenvalues and eigenvectors:
valp, vectp = np.linalg.eig(A) 

#Hermitian matrices methods:
S = np.array([[1, 2], [2, 3]])
y = np.linalg.eigvalsh(S) # eigenvalues
# eigenvalues and eigenvectors:
valp, vectp = np.linalg.eigh(S)
  
#Singular Value Decomposition
U,s,V=np.linalg.svd(A)           
Ap = np.matrix(U)*np.diag(s)*V
print(A-Ap)

copie de tableaux numpy

**Attention encore une fois !!!**

In [ ]:
x = np.array([[8,3,2],[5,1,0],[9,7,1]])
y = x
x[0,0]+=1
x[0,0]-y[0,0] # 0
z=x.copy()
x[0,0]+=1
x[0,0]-z[0,0] # 1

Génération de variables aléatoires discrètes avec numpy

  • import numpy.random as npr
  • my_sample = npr.ma_loi(paramètres, taille_du_tableau)
  • npr.randint(low=a,high=b,size=n) : v.a. unif. sur [ a, b[
  • npr.choice([a1,...,an],p=[p1,...,pn],size=n) : tirages indép. dans [a1,...,an] de loi [p1,...,pn]
  • npr.permutation(mon_urne) : permutation de mon_urne
  • npr.binomial(N,p,size=n)
  • npr.geometric(p,size=n)
  • npr.multinomial(n,tableau_des_probas,size=n)
  • npr.poisson(alpha,size=n)

Beaucoup d’autres exemples sur http://docs.scipy.org/doc/numpy/reference/routines.random.html

Génération de variables aléatoires continues avec numpy

  • import numpy.random as npr
  • my_sample = npr.ma_loi(paramètres, taille_du_tableau)
  • npr.rand(d1,d2,...) : tableau d1 x d2 x... de v.a.i. unif. sur [0, 1]
  • npr.uniform(low=a,high=b,size=n) : v.a.i. unif. sur [a, b[ (size=n peut être remplacé par size=(d1,d2,...), comme partout dans ce qui suit)
  • npr.randn(d1,d2,...) : tableau d1 x d2 x... de v.a.i. $\mathcal{N}(0,1)$
  • npr.multivariate_normal(mean=V,cov=C,size=n) : vecteurs aléatoires indépendants de loi $\mathcal{N} (V, C)$ rangés dans un tableau de taille $n\times N$,où $N$ est la taille de $V$ et $N\times N$ celle de $C$
  • npr.exponential(scale=s,size=n) : v.a.i. exponentielles de moyenne s

Beaucoup d’autres exemples sur http://docs.scipy.org/doc/numpy/reference/routines.random.html

Fonctions de numpy utiles en proba/stat

  • np.mean(x), np.std(x), np.percentile(x) : moyenne, écart-type et percentile d’un vecteur x (échantillon)
  • np.sum(x) somme des valeurs de x
  • np.cumsum(x) vecteur [x1,x1 +x2,...,x1 +···+xn] des sommes cumulées des coordonnées x1, . . . , xn de x
  • np.cov(x) matrice n × n de covariance des lignes du tableau x de taille n × p
  • scipy.stats : bibliothèque proposant densités, fonctions de répartition, quantiles, etc... de lois classiques. Cf http://docs.scipy.org/doc/scipy/reference/stats.html
  • matplotlib.pyplot bibliothèque d’affichage graphique

Boucles vs programmation matricielle

-> Eviter si possible les boucles en Python !!!

In [ ]:
from time import time
In [ ]:
n = int(1e7)
# Methode 1. Boucle for
t1 = time()
gamma1=sum([1./i for i in range(1,n+1)]) - np.log(n)
t2 = time()
temps1 = t2 - t1
# Methode 2. Numpy
t1 = time()
gamma2=np.sum(1. / np.arange(1,n+1)) - np.log(n)
t2 = time()
temps2 = t2 - t1
print("Facteur de gain: ", temps1/temps2)
In [ ]:
from timeit import timeit
N = 100
setup = """""
import numpy as np
n = int(1e5)
"""""
code_boucle = """
np.sum([1. / i for i in range(1, n)]) - np.log(n)
"""
time_boucle=timeit(code_boucle,setup=setup,number=N)

code_numpy = """
np.sum(1. / np.arange(1, n)) - np.log(n)
"""
time_numpy=timeit(code_numpy,setup=setup,number=N) 

print("Facteur : {}".format(time_boucle/time_numpy))

Affichage graphique avec matplotlib

  • import matplotlib.pyplot as plt
  • plt.plot(x,y) affiche la courbe affine par morceaux reliant les points d’abscisses x et d’ordonnées y (nombreuses options) pour x, y vecteurs de même dimension,
  • plt.hist trace un histogramme. Deux options pour les colonnes : bins= nombre de colonnes ou bins= abscisses des séparations des colonnes
  • plt.bar trace un diagramme en bâtons
  • plt.scatter(x,y) affiche le nuage de points d’abs. x et d’ord. y
  • plt.stem(x,y) affiche des barres verticales d’abs. x et hauteur y
  • plt.axis([xmin,xmax,ymin,ymax]) définit les intervales couverts par la figure
  • plt.axis(’scaled’) impose que les échelles en x et en y soient les mêmes
  • plt.show() affiche les fenêtres créées dans le script
  • plt.figure() crée une nouvelle fenêtre graphique
  • plt.title("mon titre") donne un titre à une figure
  • plt.legend(loc=’best’) affiche la légende d’un graphique (en position optimale)
  • plt.subplot subdivise la fenêtre graphique de façon à y afficher plusieurs graphiques

Exemple : représentation d’un échantillon de loi discrète

In [ ]:
import scipy.stats as sps
import matplotlib.pyplot as plt
%matplotlib inline
n, p, N = 20, 0.3, int(1e4)
B = np.random.binomial(n, p, N)
f = sps.binom.pmf(np.arange(n+1), n, p)
plt.hist(B,bins=n+1,density = 1,range=(-.5,n+.5),color="white",label="loi empirique")
plt.stem(np.arange(n+1),f,"r",label="loi theorique")
plt.legend()

Histogramme d’un échantillon de loi continue

Exercice

  1. Créer un vecteur E contenant 10000 réalisations indépendantes d'une loi $\mathcal{N}(0,1)$.
  2. Aficher sur le même graphique ce vecteur et la loi gaussienne théorique, obtenue grâce à la fonction norm.pdf de la bibliothèque scipy.
In [ ]:
E = np.random.randn(int(1e5))#echantillon
x = np.linspace(-4,4,1000)
f_x = sps.norm.pdf(x) #Densite gaussienne
plt.plot(x,f_x,"r",label="Theory")
#Affichage histo:
plt.hist(E,bins=50,density=1,label="Data")
plt.legend(loc='best')

La loi des grands nombres

Theorème, (Loi des Grands Nombres), Kolmogorov, 1929

Soient $(X_i)_{i\ge 1}$ copies indépendantes d'une même variable aléatoire $X$:\begin{eqnarray*} \mathbb{E}[|X|] <\infty&\implies& S_n:=\frac{X_1+\dots+X_n}{n}%\underset{n\to\infty}{\stackrel{p.s.}{\longrightarrow}} \underset{n\to\infty}{\longrightarrow} \mathbb{E}[X],\\ \mathbb{E}[|X|] =\infty&\implies& S_n:=\frac{X_1+\dots+X_n}{n}\;\textrm{ diverge}.\end{eqnarray*}

Exercice

  1. Tirer un vecteur $X$ dont les coordonnées sont $n$ réalisations indépendantes de la loi uniforme sur $[0,1]$. Afficher sur le même graphique la courbe $S_n$ en fonction de $n$ et une droite qui vaut $\mathbb{E}[X]$ pour tout $n$. Commentez le résultat.
In [ ]:
wn=10000
S=np.cumsum(np.random.rand(n))/np.arange(1,n+1)
plt.plot(range(1,n+1),S,'r',label="S_n")
plt.plot((1,n),(.5,.5),"b--",label="Esperance")
plt.ylabel('S_n')
plt.xlabel("n")
plt.legend(loc='best')
plt.title("LGN")
  1. On va maintenant se placer dans le cas d'une distribution dont l'espérance n'est pas forcément finie. Soient $s,U$ v.a. indépendantes, $U$ uniforme sur $[0,1]$ et $s=\pm 1$ avec probas 0.5,0.5. Alors la v.a. $X:=s U^{-1/\alpha}$, appelée ici $\alpha$-variable aléatoire, a pour densité $ (\alpha/2)\mathbf{1}_{|x|\ge 1} |x|^{-\alpha-1}$ et est d'espérance finie ssi $\alpha>1$. Afficher à nouveau $S_n$ en fonction de $n$.
In [ ]:
alpha, n = 0.9, 1000
U = np.random.rand(n)
X = (2*np.random.randint(0,2,n)-1)*U**(-1/alpha)
S=np.cumsum(X)/np.arange(1,n+1)
plt.title("LGN: convergence ou pas selon alpha")
plt.plot((0,n),(0,0),"b--")
plt.plot(S,"r",label="Moyenne empirique")
plt.xlabel("n")
plt.ylabel("S_n")
plt.legend(loc='best')
  1. On se replace dans le cas où les $X_i$ suivent une loi uniforme sur $[0,1]$. On peut se demander maintenant à quelle vitesse $S_n$ converge vers $\mu = \mathbb{E}[X]$ ? On cherche donc $\beta$ tel que ${n^\beta(S_n-\mu)\longrightarrow \ell\ne 0}$, i.e. $${S_n\approx \mu +\frac{\ell}{n^\beta}}.$$ Pour $\beta$ fixé, afficher sur un même graphique plusieurs courbes $(n, n^\beta(S_n-\mu))$ à l'aide d'une boucle for. Essayez par exemple $\beta =0.2$, $\beta = 0.5$, $\beta = 0.7$.
In [ ]:
beta = 0.8
n=1000
for k in range(10):
    S=np.cumsum(np.random.rand(n))/np.arange(1,n+1)
    T = (S - 0.5)*np.arange(1,n+1)**beta
    plt.plot(range(1,n+1),T)
    plt.plot((1,n),(0,0),"b--")

plt.ylabel('n^beta(S_n-mu)')
plt.xlabel("n")
plt.legend(loc='best')
plt.title("n^beta (erreur dans la LGN)")

Théorème de la limite centrale, Laplace, 1812, Lindeberg, 1920

$(X_i)_{i\ge 1}$ v.a. indépendantes de même loi d'espérance $\mu$ et d'écart-type $\sigma$: $$T_n = \frac{n^{1/2}}{\sigma}\left(\frac{X_1+\dots+X_n}{n}-\mu\right)\underset{n\to\infty}{\stackrel{loi}{\longrightarrow}} \mathcal{N}(0,1). $$

En d'autres termes, la moyenne empirique des $X_i$ $\approx$ moyenne théorique, avec une erreur aléatoire gaussienne d'ordre ${1/\sqrt{n}}$.

  1. Tirer aléatoirement $m$ vecteurs $(X_1,\dots, X_n)$ avec tous les $X_i$ i.i.d. de loi uniforme sur $[-1,1]$ et afficher l'histogramme des $m$ valeurs de $T_n$ ainsi obtenues. Afficher sur le même graphique la loi gaussienne centrée d'écart-type $\sigma$.
In [ ]:
n,m,sigma=1000,10000,3**(-1/2)
X=2*np.random.rand(m,n)-1
S=np.sum(X,axis=1)/(np.sqrt(n)*sigma)
M=max(np.abs(S))
x=np.linspace(-M,M,1000)
y=sps.norm.pdf(x)
plt.plot(x,y,'r',label="densite")
plt.hist(S,bins=int(round(m**(1./3)*M*.5)),density=1,histtype='step',label="Histogramme")
plt.legend(loc='best')
plt.title("TCL")

Interprétation

Généralement, une grande somme de petits aléas peu corrélés fluctue autour de sa moyenne selon une distribution gaussienne.

Exemples : anatomie (ex : taille des individus de sexe donné), QI, nombreuses mesures physiques, données économiques, etc...

Lecture et écriture dans un fichier externe

  • x=range(5)
  • mon_flux=open("my_data.txt","w") #w=write
  • mon_flux.write(str(x))
  • mon_flux.close()

Ecriture a la fin d'un fichier existant:

  • mon_flux=open("my_data.txt","a") #a=append
  • mon_flux.write("\n"+str(x+1))
  • mon_flux.close()

Lecture:

  • mon_flux=open("my_data.txt","r") #r=read
  • y=mon_flux.read()
  • print(y)

Voir aussi np.save, np.load, np.savetxt, csv.reader...

  • import numpy as np
  • x=np.random.rand(5,3)
  • np.save("my_npy_file.npy",x) #creation et ecriture
  • x2=np.load("my_npy_file.npy") #lecture

Lecture dans un fichier externe : exemple

In [ ]:
f=open("PopLynxRegionCanada_1821_1934.dat","r")
ytxt=f.readlines()                       # list de str
y=[int(row) for row in ytxt]             # convertit str en int
plt.plot(range(1821,1935),y,"r")
plt.title("Population de lynx")
plt.tight_layout()                       # pratique pour l’export
plt.show()

Python pour l’analyse : exemples

In [ ]:
import numpy as np
import scipy
import matplotlib.pyplot as plt
def f(x):
    return np.exp(x)+x
#zeros of f, computed starting at -0.2:
a=scipy.optimize.fsolve(f,-.2)
print(a,f(a))
#integral of f from 0 to 1:
b=scipy.integrate.quad(f,0,1)
print(b)
def g(y,t):
    return y
T=np.arange(start=0,stop=1,step=.001)
##solution, at T, of y’=g(y,t), y(T[0])=1:
y=scipy.integrate.odeint(g,1,T)
plt.plot(T,np.log(y),"r")
plt.show()

Bibliographie : les tutoriels/sites officiels