#!/usr/bin/env python
# coding: utf-8

# <a href="https://colab.research.google.com/github/JLaurentin/notebooks_BCPST2_EDEN/blob/master/td0BCPST2b_LesPrincipauxModules.ipynb" target="_parent"><img src="https://colab.research.google.com/assets/colab-badge.svg" alt="Open In Colab"/></a>

# # Python : Les principaux modules utilisés en BCPST
# 
# <img src="https://avatars3.githubusercontent.com/u/7388996?v=3&s=200">
# 
# Ce notebook est proposé dans le cadre d'une préparation au ** Cours de mathématiques & informatique ** en B.C.P.S.T.2 - Externat des enfants nantais proposé par [J. Laurentin](<j.laurentin@worldonline.fr>).

# La  bibliothèque standard de Python possède plus de deux cents bibliothèques au sein desquels vous trouverez des fonctions "prêtes à l'emploi", regroupées par thèmes qu'on appelle des "modules".  
# Les principaux modules dont nous aurons besoin sont les modules `math`, `numpy`, `scipy`, `matplotlib` et `random`. Il est possible d'accéder à la documentation correspondante en exécutant `help(nom_module)`.   
# Nous aurons l'occasion d'en recontrer bien d'autres au cours de l'année et ils seront introduits au fur et à mesure des besoins.  
# N'hésitez pas, à jeter un oeil à la description des [modules disponibles](http://docs.python.org/py3k/library/index.html) ainsi qu'à ceux développés par la communauté [Python](http://pypi.python.org/pypi)

# **Attention !** Si vous faites ce notebook par "petites touches" (ce que je vous recommande), vous pouvez bien sûr le sauvegarder pour le reprendre là ou vous l'avez laissé. Pour autant, les bibliothèques dont vous avez besoin doivent **à nouveau** être importées. Il faudra donc ré-exécuter toutes les lignes d'importation de bibliothèque à chaque nouvelle cession du notebook.

# **COMMENCER PAR ENREGISTRER CE NOTEBOOK SUR VOTRE DRIVE**

# ## 1. Le module math

# In[ ]:


from math import *


# On n'hésitera pas à consulter la documentation associée à ce module.  
# En premier lieu on se familiarisera avec celle qui énumère l'ensemble des [fonctions disponibles](https://docs.python.org/3/library/math.html) :  Elles sont réparties entre plusieures rubriques, *théorie des nombres*, *fonctions puissances et logarithmiques*, *fonctions trigonométriques*, *conversions angulaires* (`math.degrees()` et `math.radians()`) mais aussi quelques constantes usuelles :  
# - `pi`, `e`, `inf` et `nan`  
# 
# On accèdera à l'aide de chaque fonction en exécutant : `help(nom_fonction)`.

# In[ ]:


help(cos)


# On utilisera ce module, notamment pour le calcul des principales fonctions mathématiques au sein d'une liste.  
# *Remarque :* Si le module `numpy` est par ailleurs importé, il est inutile d'utiliser le module `math` car on retrouve dans numpy les mêmes fonctions, sans compter qu'il n'est pas utile de démultiplier les modules importés.
# 
# A titre d'exemple, on pourra utiliser le module `math` pour déterminer l'ensemble des termes d'une suite définie de façon explicite.  
# Ainsi, si $u_n=sin(\frac{1}{n})e^{n/2}$, $\forall n\in\mathbb{N}^*$, alors la liste `Lu` des $10$ premiers termes de la suite vaut :

# In[ ]:


Lu = [sin(1/n)*exp(n) for n in range(1,11)]
print(Lu)


# **Les modules associés complémentaires :** Le module math est particulièrement adapté à l'utilisation sous forme de simple *calculette* de Python. Dans ces conditions, il est bon, si nécessaire, de l'associer avec deux modules qui permettent de travailler avec les nombres complexes ou encore avec les fractions.

# ### a. le module cmath

# Ce module permet de travailler sur les nombres complexes.  
# Pour autant, il est possible de s'en passer pour les opérations d'arithmétique élémentaire.  
# Ainsi, si $z_1=2+3i$ et $z_2=1+i$  qui seront notés comme en physique $z_1=2+3j$ et $z_2=1+j$, il n'y a besoin d'aucun module complémentaire pour déterminer leur module, leur conjugué, leur partie réelle, imaginaire ou encore pour les additionner ou les multiplier.  
# On notera qu'il y a deux façons de définir les nombres complexes : $z = a+b.j$ ou z=complex(a,b).  
# **Exemple :**

# In[ ]:


z1 = 2+5.j ; z2 = complex(2,1)
print(z1,z2)
print('le module de z1 vaut : ',abs(z1))
print(abs(z1)==sqrt(13))
print('le conjugué de z1 vaut : ',z1.conjugate())
print('la partie réelle de z1 : ',z1.real)
print('la partie imaginaire de z1 :',z1.imag)
print('la somme z1+z2 vaut : ',z1+z2)
print('le produit z1.z2 vaut :',z1*z2)


# *Entraînez-vous en modifiant les valeurs de z1 et z2 !*

# Pour autant, il n'est pas possible de calculer l'argument du nombre complexe (appelé *phase* dans Python), de passer de la forme arithmétique à la forme en polaire, ou réciproquement.  
# Pour cette raison, il faut importer le module cmath.

# In[ ]:


from cmath import *
z2 = complex(1,1)
print(z2)
print("l'argument de z2 vaut : ",phase(z2))
print("le module et l'argument de z2 (forme polaire) sont :",polar(z2))
# Et si z = rho.exp(j theta) on obtiendra sa forme cartésienne en écrivant :
z3 = rect(2,pi/3)
print('la forme cartésienne de z3 est : ',z3)


# On se méfiera des problèmes d'approximation liés à la représentation binaire des nombres en machine.  
# ainsi, pour $z=e^{2i\pi/3}=-\cfrac{1}{2}+i\cfrac{\sqrt{3}}{2}$, on obtient :

# In[ ]:


t = complex(0,2*pi/3)
z = exp(t)
print(z)
print(z.real+1/2) # doit valoir 0
print(z.imag-sqrt(3)/2) # doit également valoir 0


# ### b. le module fractions

# Ce module permet de travailler avec les nombres rationnels, de les additionner, de les multiplier, de mettre en évidence des formes irréductibles.  
# Voici des exemples inspirés de la documentation Python pour ce module.  
# On notera la souplesse d'utilisation de la méthode Fraction qui permet de créer un nombre rationnel en prenant, en paramètre d'entrée, aussi bien des couples d'entiers (numérateur/dénominateur) que des flottants ou encore des chaînes de caractères...  

# In[ ]:


from fractions import *


# In[ ]:


print(Fraction(16, -10))
print(Fraction(123))
print(Fraction(123, 1))
print(Fraction(2.25))
print(Fraction('3/7'))
print(Fraction(' -3/7 '))
print(Fraction('1.414213'))
print(Fraction('-.125'))
print(Fraction('7e-6'))


# Il est aussi possible d'obtenir des approximations rationnelles de nombres irrationnels, à la précision souhaitée :

# In[ ]:


print(Fraction(pi).limit_denominator(1000))


# Mais il y a aussi quelques pièges... il faut en particulier pouvoir lever des problèmes d'approximations liées aux représentation sous forme binaire des nombres de type float...

# In[ ]:


from math import cos
print(Fraction('1.1'))
print(Fraction(1.1)) # 1.1 est de type float.
print(Fraction(1.1).limit_denominator())
print(Fraction(cos(pi/3)))
print(Fraction(cos(pi/3)).limit_denominator())


# ## 2. Le module numpy

# Ce module permet de travailler sur un environnement matriciel, les objets manipulés étant considérés au choix comme des tableaux (de type array) ou des matrices si on souhaite se restreindre à des calculs d'algèbre linéaire.  

# In[ ]:


import numpy as np


# Il est possible d'accéder à l'ensemble des fonctions proposées par Python au sein d'une bibiothèque.  
# A titre d'exemple, si ci-dessous, vous placez votre curseur juste après ``np.random.`` et que vous pressez la touche  **Tab** (tabulation, située en haut à gauche de votre clavier), vous allez voir la liste des façons de compléter votre code par une fonction du sous-module `np.random`.

# In[ ]:


np.random.


# Si vous tapez une parenthèse ouvrante suivie d'un touche **Tab** après une fonction quelconque du module, vous verrez apparaître l'aide sur cette fonction fournie par la documentation... 

# In[ ]:


np.random.rand(


# Pour ouvrir l'aide dans une fenêtre indépendante, au bas de votre écran, ajoutez un **?** juste après l'objet ou le nom de la méthode et exécutez la cellule à l'aide d'un **Shift+Entrée**:

# In[ ]:


get_ipython().run_line_magic('pinfo', 'np.random')


# Et maintenant un peu de pratique...

# In[ ]:


L = [[1,2],[3,4]]
T1 = np.array(L)
print(T1)


# In[ ]:


type(T1)


# In[ ]:


M1 = np.matrix(L)
print(M1)


# In[ ]:


type(M1)


# La manipulation des tableaux est proche de celle des listes et offre la même souplesse.  
# On perd évidemment la possibilité d'appliquer les méthodes associées aux listes (comme `append` ou `count`) mais, en contrepartie, dans le cas des tableaux ou des matrices, on gagne en facilité d'écriture pour accéder à certains éléments ou sous-matrices. Sans compter qu'il est possible de faire appel à des fonctions de la bibliothèque numpy pour multiplier, inverser des matrices ou encore résoudre des systèmes...

# ### 2.1. Création de matrices.

# On peut créer des **matrices lignes** à l'image des listes en utilisant la fonction `arange()` ou encore une matrice ligne de $N$ éléments, compris entre $a$ et $b$ grâce à la fonction `linspace(a,b,N)`.  
# *Exemples :*

# In[ ]:


print(np.arange(10))
print(np.arange(1,10))
print(np.arange(1,10,2))
print(np.linspace(1,10,10)) # Vérifier que a=1 et b=10 sont compris dans ce tableau.
print(np.linspace(1,10,20))


# Il est également possible de travailler en dimension deux.  
# Des matrices prédéfinies existent et peuvent être crées rapidement :

# In[ ]:


A1 = np.zeros((3,4))
A2 = np.ones((3,3))
A3 = np.diag(np.arange(3))
Id = np.diag(np.ones(3))
print(A1,'\n\n',A2,'\n\n',A3,'\n\n',Id) # le '\n' permet d'insérer un passage à la ligne 


# Ou elle peuvent être créées manuellement. On commencera pour ça par créer une liste qu'on transformera au choix en tableau (`np.array(L)`) ou matrice (`np.matrix(L)`) selon l'objectif poursuivi.  
# *Exemples :*

# In[ ]:


L = [[1,2,3],[4,5,6],[7,8,9]]
A = np.array(L)
print('A = ',A)


# In[ ]:


L=[[1,2],[3,4]]


# In[ ]:


M1 = np.matrix(L)
print('M1 = ',M1)


# ### 2.2. Manipulation des matrices.  
# Sauf mention contraire, les lignes de code qui suivent s'appliquent aussi bien à A qu'à M1.

# **Récupération d'une valeur :**

# In[ ]:


print(A)
A23 = A[1,2]
print('A23 vaut ',A23)


# Entraînez-vous avec d'autres valeurs...

# **Extraction d'une colonne :**

# In[ ]:


C = A[:,0] # première colonne
print(C)


# Ou d'une partie d'une colonne (ci-dessous 2è et 3è coefficients de la colonne 2)

# In[ ]:


A[1:3,1]


# **Extraction d'une ligne ou de plusieurs lignes :**

# In[ ]:


L2 = A[1,:]
print(L2)


# In[ ]:


L_12 = A[:2,:]
print(L_12)


# ### 2.3. Quelques opérations possibles sur les matrices.

# On trouve dans le module numpy l'ensemble des opérations algébriques sur les matrices :

# In[ ]:


print(A)


# In[ ]:


# Addition de matrices :
print(A+Id)


# Mais **attention** ! La matrice $A$ est de type `numpy.array`. Les opérations qui s'y réfèrent **ne sont pas** celles de l'algèbre linéaire mais sont définies "termes à termes".  
# Ainsi, $2*A$ multiplie chacun des termes de $A$ par $2$, $A+1$ ajoute $1$ à chacun des termes ou encore $A*Id$ ne calcule pas le produit matriciel mais retourne le produit termes à termes.  
# Entraîner vous !

# In[ ]:


print(A*Id)


# Dans ce cas comment faire le produit matriciel ??
# 1. Si on travaille avec des type `array` on exécutera `np.dot(A,B)` ou encore `A@B`.
# 2. Si on travaille avec des type `matrix` il suffit de faire `A*B`.

# In[ ]:


A = np.array([[1,2],[3,4]]) ; B=np.array([[2,0],[0,2]])
print('produit 1 :\n',np.dot(A,B))
print('produit 2 :\n',A@B)


# In[ ]:


M = np.matrix([[1,2],[3,4]]) ; N=np.matrix([[2,0],[0,2]])
print(M*N)


# Mais on peut faire bien d'autres choses que ça sur les matrices, notamment si on on travaille avec le module dédié consacré à l'algèbre linéaire, à savoir `linalg` (n'hésitez pas à utiliser `help(np.linalg)`)   
# Essayez dans la ligne de code ci-dessous, avec la matrice `M1` les fonctions suivantes :
# 1. `np.transpose(M)` ou `M.T`
# 2. `M1 = np.linalg.inv(M)` puis `M1*M` 
# 3. `r = np.linalg.matrix_rank(M)`. *Note :* $M$ étant inversible, son rang vaut $2$.
# 4. `X = np.linalg.solve(M,b)` avec `b=np.matrix([[1],[1]])` pour la résolution du système $MX=b$ : $\begin{cases} x+2y&=1\\3x+4y&=1\end{cases}$

# In[ ]:


print(M1)


# ### 2.4. Les applications possibles en statistiques :

# Une série statistique étant entrée sous forme de tableau, la bibliothèque numpy permet de calculer l'ensemble des paramètres de position et de dispersion mais également, dans le cas de séries doubles, les coefficients de corrélation ou encore les coefficients de la droite de régression.  
#   
# Prenons l'exemple suivant, tiré de *Mathématiques et statistiques pour les sciences de la nature*, G. Biau, J. Droniou, M. Herslich (*E.D.P. Sciences*, 2010), p. 250

# "Une équipe de chimistes a mis au point un alliage de fer et de carbone capable de résister à des conditions extrêmes et souhaite modéliser sa résistance en fonction de la teneur en carbone".  
# Dix teneurs en carbonnes ont été testées (exprimées en nombre d'unités pour 100 :  
# $x_i$ : 0.3,0.7,1,1.1,1.2,2.5,2.7,3,5,5.2  
# pour lesquelles la "charge de rupture" du nouvel alliage (en tonnes) a été déterminée :  
# $y_i$ : 0.8,1.77,2.13,2.42,2.87,5.34,5.5,6.41,10.31,10.47  
# 1. Déterminer la moyenne, la médiane, le premier et troisième quartile, la variance de chacune des séries.
# 2. Calculer le coefficient de corrélation ainsi que les coefficients de la droite de régression.

# In[ ]:


# Pensez à importer le module numpy si vous ne l'avez pas encore fait !
TX = np.array([0.3,0.7,1,1.1,1.2,2.5,2.7,3,5,5.2])
TY = np.array([0.8,1.77,2.13,2.42,2.87,5.34,5.5,6.41,10.31,10.47])
print('Pour X, les principaux paramètres de position et de dispersions sont :')
print('le minimum de X vaut : ',np.min(TX),' et son maximum vaut : ',np.max(TX))
print('la moyenne de X vaut : ',np.mean(TX))
print('la médiane de X vaut : ',np.median(TX))
print('le premier quartile vaut :',np.percentile(TX,25))
print('le troisième quartile vaut :',np.percentile(TX,75))
print('la variance de X vaut : ',np.var(TX))
print("l'écart-type de X vaut : ",np.std(TX))
print('Pour la série multivariée (X,Y) on obtient les paramètres suivants :')
print('La coefficient de corrélation du couple (X,Y) vaut : ',np.corrcoef(TX,TY)[0,1]) # matrice symétrique
print('Les coefficients a et b de la droite de régression y = ax+b sont ',np.polyfit(TX,TY,1))


# ### 2.5. Le calcul polynômial.

# Le module numpy permet de travailler simplement sur les polynômes. Toutes les opérations au programme sont possibles : détermination du degré, multiplication d'un polynôme par un scalaire, addition et somme de deux polynômes, recherche des racines, dérivation.  
# **Attention :** Un polynôme $P=\displaystyle\sum_{k=0}^na_kX^k=a_0+a_1X+\cdots+a_nX^n$ est caractérisé, sous Python, par la liste de ses coefficients **d'indices décroissants**.  
# A titre d'exemple, nous allons travailler ci-dessous avec $P=X^2-2X+1$, $Q=X^2+2X+3$ et $R=X^2+1$.  

# In[ ]:


# votre module numpy est importé ?
import numpy as np


# In[ ]:


P = np.poly1d([1,-2,1])
print(np.poly1d(P))


# In[ ]:


Q = np.poly1d([1,2,3],variable='Z')
print(np.poly1d(Q))


# On peut aussi **définir un polynôme à partir de ses racines** si on ajoute un second paramètre égale à True.  
#   
# *Par exemple :*

# In[ ]:


P2 = np.poly1d([1,-2,1],True)
print(np.poly1d(P2))


# Il s'agit bien du polynôme $P_2(X)=X^3-3X+2$ dont $-2$ est racine simple car $P_2(-2)=0$ et $1$ est racine double car $P_2(1)=P_2'(1)=0$, soit $P_2(X)=(X-1)^2(X-2)$.

# Pour le calcul de $P(a)$ où $a\in\mathbb{R}$, il suffit d'appeler $P(a)$ mais **surtout pas** $P[a]$ qui désigne, pour peu que $a$ soit un entier, la valeur du coefficient d'indice $a$ du polynôme $P$ interprété comme une liste.

# In[ ]:


print(Q[1]) # deuxième coefficient de Q
print(Q(1)) # valeur de Q(1)
print(np.polyval(Q,1)) # même chose que précédemment
print(np.polyval(Q,np.arange(11))) # Pour évaluer le polynôme Q en 0,1,2,...,10. Résultat affiché sous forme de liste.


# In[ ]:


print(P.order) #degré de P
print(P.c) # coefficients de P - sous forme de liste.
print(P.r) # racines de P


# Essayez avec les autres polynômes.  
# Les racines complexes sont bien évidement retournées ! Le vérifier avec $R(X)=X^2+1$

# Quant à la **dérivation**, il suffit de faire appel à la méthode 'deriv' qui dérive à l'ordre indiqué.

# In[ ]:


P = np.poly1d([2,1,-3,4])
Pp1=np.poly1d.deriv(P,m=1)
Pp2=np.poly1d.deriv(P,m=2)
Pp3=np.poly1d.deriv(P,m=3)
print(np.poly1d(P))
print(np.poly1d(Pp1))
print(np.poly1d(Pp2))
print(np.poly1d(Pp3))


# ## 3. Le module matplotlib.pyplot

# Faire des calculs, obtenir des listes ou des tableaux de résultats est une chose mais que ce soit pour émettre des hypothèses, appuyer vos commentaires ou discuter les valeurs obtenues, il est indispensable de pouvoir proposer des représentations graphiques.  
# On importera pour ça le module `matplotlib.pyplot` qui offre une large palette d'affichages possibles et se prête à de nombreux champs d'applications (analyse,statistiques,traitement d'images, etc.)
# 

# In[ ]:


# Ces lignes configurent matploblib pour montrer les figures "encapsulées" dans le notebook JUPYTER
# au lieu d'ouvrir une nouvelle fenêtre pour chaque figure.  
# CETTE CELLULE NE DOIT PAS ÊTRE EXECUTEE sous Pyzo ou sur colab.research.google.com où on se contentera, dans la cellule suivante, 
# d'importer matplotlib.pyplot
get_ipython().run_line_magic('matplotlib', 'inline')
#%matplotlib notebook


# In[ ]:


import matplotlib.pyplot as plt
# Assurez-vous d'avoir importé le module numpy as np...!


# ### 3.1. Les représentations graphique en analyse :

# #### a) Les suites numériques

# la suite numérique de réels $u_n = \cfrac{1}{3}$

# Donnons une représentation des dix premiers termes de la suite géométrique : $u_{n}=3\left(\cfrac{2}{3}\right)^n$, $\forall n\in\mathbb{N}$ :

# In[ ]:


U = [3*(1/3)**n for n in range(10)]
plt.plot(np.arange(10),U,'ro')
plt.title('Représentation des 10 premiers termes de $(u_n)$')
plt.xlabel('n')
plt.grid()
plt.ylabel('$u_n$')


# Pour comprendre la deuxième ligne de la cellule qui précède, on se reportera à la documentation de la fonction `plot()` en exécutant `help(plt.plot)` dans la cellule ci-dessous.  
# 
# **Exercice 1 :** Modifiez ensuite la précédente cellule pour tracer les termes de la suite avec des carrés verts à la place des ronds rouges.

# In[ ]:


help(plt.plot)


# **Exercice 2 :** Reprenez le script précédent et, en vous appuyant sur vos T.D. de BCPST1, tracer les premiers termes de suites arithmétiques, géométriques, arithmético-géométriques, récurrentes linéaires d'ordre 2 dont vous avez obtenus la forme explicite. Faire des hypothèses graphiques sur leur convergence ou leur divergence que vous validerez par un calcul de limite.

# #### b) Les représentations graphiques de fonctions d'une variable réelle.

# Rappelons en préalable qu'il est possible de définir de façon simple et rapide de telles fonctions grâce au mot clé `lambda`.  
# *Exemple avec* $f_1$ *définie sur* $\mathbb{R}$ par $f_1(x)=e^x-1$ :

# In[ ]:


f1 = lambda x:np.exp(x)-1


# In[ ]:


print(f1(1))


# In[ ]:


X = np.linspace(0,2,10)
f1(X)


# **Attention !** Si vous aviez défini `f1` par `f1 = lambda x:exp(x)-1` vous auriez eu dans le calcul de `f1(X)` un message d'erreur car `X` est de type `np.array` et la fonction exponentielle utilisée doit également appartenir au module numpy.

# Il est désormais possible de tracer l'allure de la courbe représentant `f1` sur l'intervalle $[a,b]$ de son choix :

# In[ ]:


a,b = -1,2
X = np.linspace(a,b,100)
plt.plot(X,f1(X),'b')
plt.grid()
plt.title('Graphe de $f_1$')
plt.xlabel('x')
plt.ylabel('y')


# Il est possible de tracer le graphe de deux fonctions dans une même fenêtre et leur affecter des caractéristiques de couleur, d'épaisseur de trait, de style, etc. différents.  
# Deux méthodes sont possibles que nous allons illustrer en cherchant à déterminer les graphiquement les valeurs de $x\in\mathbb{R}$ telles que $e^x-1=2x$ :

# In[ ]:


a,b = -1,2
X = np.linspace(a,b,100)
plt.plot(X,f1(X),'b',X,2*X,'r--')
plt.grid()
plt.title('Graphe de $f_1$')
plt.xlabel('x')
plt.ylabel('y')


# Le problème de la rédaction précédente est qu'elle ne permet pas de légender les courbes tracées.  
# On préfèrera donc :

# In[ ]:


a,b = -1,2
X = np.linspace(a,b,100)
plt.plot(X,f1(X),'b',label='$y=e^x-1$') # le label est une façon d'étiqueter la courbe de f1
plt.plot(X,2*X,'r--',lw=2,label='$y=2x$') # lw pour "linewidth"
plt.grid()
plt.title('Graphe de $f_1$')
plt.legend(loc='best') # d'autres choix de placements de la légendes peuvent être faits (cf help(plt.legend))
plt.xlabel('x')
plt.ylabel('y')


# Il est aussi possible de découper une fenêtre graphique en autant de sous-fenêtres qu'on le souhaite à l'aide de la fonction `subplot(nl,nc,i)` où `nl` et `nc` désignent respectivement le nombre de lignes et de colonnes de la fenêtre principale découpée comme un tableau et `i` indique le numéro de la sous-fenêtre dans laquelle on se place.  
# 
# **Exemple 1** : Deux sous-fenêtres, l'une au dessus de l'autre à l'aide de `subplot(2,1,i)`.

# In[ ]:


x = np.linspace(0,3)

y1 = np.cos(2*np.pi*x)*np.exp(-x)
y2 = np.cos(2*np.pi*x)

plt.subplot(211) # Forme abrégée de 'plt.subplot(2,1,1)'
plt.plot(x, y1, 'yo-')
plt.title("Exemple avec subplot 2 lignes, 1 colonne")
plt.ylabel('Oscillations amorties')

plt.subplot(212) # Forme abrégée de 'plt.subplot(2,1,2)'
plt.plot(x, y2, 'r.-')
plt.xlabel('time (s)')
plt.ylabel('Non amorties')


# **Exemple 2** On cherche, dans quatre fenêtres distinctes, à confronter le graphe de la fonction $x\longmapsto e^x$ avec la partie régulière de ses développements limités au voisinage de $0$ à l'ordre $1$, $2$, $3$ et $4$.

# In[ ]:


x = np.linspace(-2,2,50)
y = np.exp(x)
P1 = lambda t:1+t
P2 = lambda t:1+t+t**2/2
P3 = lambda t:1+t+t**2/2+t**3/6
P4 = lambda t:1+t+t**2/2+t**3/6+t**4/24
couleur = ['r','g','b','m']
plt.subplots(figsize=(12, 8)) # permet de définir la taille de la fenêtre principale
for i in range(1,5):
    plt.subplot(2,2,i)
    plt.plot(x,y,'k--',lw=3)
    plt.plot(x,eval('P'+str(i))(x),couleur[i-1])
    plt.title("D.L. au voisinage de 0 à l'ordre "+str(i))


# **Exercice :** Rappeler le développement limité à l'ordre $3$ au voisinage de $0$ de $f:x\longmapsto\ln(1+x)$ et compléter le code ci-dessous pour tracer dans trois fenêtres **consécutives** alignées le graphe de $f$ confronté à celui des parties régulières de ses développements limités à l'ordre $1$ (première fenêtre), à l'ordre $2$ (deuxième fenêtre) et à l'ordre $3$ (troisième fenêtre).

# In[ ]:


x = np.linspace(-0.8,2)
y = # à compléter
P1 = lambda t:t
P2 = lambda # à compléter
P3 = lambda # à compléter
plt.subplots(figsize=(12, 4)) # entraînez-vous en supprimant cette ligne ou en modifiant les valeurs.
plt.subplot(131)
plt.plot(x,y,'k',x,P1(x),'r')
plt.ylabel('y')
plt.subplot(# à compléter)
plt.plot(x,y,'k',x,P2(x),'g')
plt.ylabel('y')
plt.subplot(# à compléter)
plt.plot(x,y,'k',x,P3(x),'b')
plt.ylabel('y')
plt.tight_layout() # mettre cette ligne en commentaire et relancer pour vérifier que les textes se chevauchent.


# Puis on peut compléter ses graphes par du texte, à l'intérieur ou à l'extérieur de la figure, ou bien même colorer certaines partie d'un graphique comme le montre ci-dessous un exemple de représentation d'une aire sous la courbe associées à un calcul intégrale ()

# In[ ]:


from matplotlib.patches import Polygon
# d'après Pyx - Python Graphic Package 
# exemple repris dans http://jeffskinnerbox.me/notebooks/running-code-in-the-ipython-notebook.html

func = lambda x:(x-3)*(x-5)*(x-7)+85
 
a, b = 2, 9 # bornes du calcul intégral
x = np.arange(0, 10, 0.01) # tableau des abscisses de 0 à 10 par pas de h = 0.01
y = func(x)
plt.plot(x, y, linewidth=1,label='$(x-3)*(x-5)*(x-7)+85$')
plt.legend(loc='best')
# Traçons en ombré l'aire sous la courbe grâce à la construction d'un polygone :
ix = np.arange(a, b, 0.01) # le tableau des abscisses concernées
iy = func(ix)
verts = [(a,0)] + list(zip(ix,iy)) + [(b,0)] # création des couples reliés par le polygone.
poly = Polygon(verts, facecolor='0.8', edgecolor='k')
plt.subplot(111).add_patch(poly) # On colle dans la fenêtre le polygone dont la surface est grisée

plt.text(0.5 * (a + b), 30,
     r"$\int_a^b f(x)\mathrm{d}x$", horizontalalignment='center',
     fontsize=16) # pour mettre le texte de son choix dans le graphe

plt.axis([0,10, 0, 180])
plt.figtext(0.9, 0.05, 'x') # Pour placer le texte à l'extérieur de la figure
plt.figtext(0.05, 0.8, 'y')


# Le calcul, quant à lui, se fait aisément grâce au module `scipy.integrate` de la manière suivante :

# In[ ]:


from scipy import integrate
func = lambda x:(x-3)*(x-5)*(x-7)+85
a,b = 2,9
I,e = integrate.quad(func,a,b)
print("L'intégrale cherchée vaut ",I," et l'erreur commise vaut :",e)


# **Double échelle des ordonnées :** Il est possible de mettre en place un double système d'axes.  
# Imaginons qu'on dispose de précipitations et de températures qu'on souhaite représenter sur une même figure :

# In[ ]:


T = [12,15,9,8,10,11.5,13,15,11,9.3]
P = [3,5,0,0,4,11,9,1,0,5]
t = np.arange(1,11) 
plt.ylabel("Précipitations", color = "b")
plt.ylim(0,12.0)
plt.bar(t,P,color="b")
plt.twinx() # second système d'axes :
plt.ylabel("Températures", color="r")
plt.ylim(0,20)
plt.plot(t, T, "o-r")


# #### c) Les nuages de points (scatter).

# In[ ]:


x = np.random.randn(500) # 500 tirages selon la loi normale centrée réduite.
y = np.random.randn(500) # idem
plt.figure(figsize=(10,5))
plt.scatter(x, y, s=100, c="PowderBlue") # taille des disques constante égale à 100
plt.show()
sizes = np.abs((100*np.random.randn(500))) # tailles au hasard entre 0 et 100
colors = np.random.randn(500)
plt.figure(figsize=(10,5))
plt.scatter(x, y, s=sizes, c=colors)


# ### 3.2. Les représentations graphiques en statistiques :

# Reprendons l'exemple 2.4 pour lequel dix teneurs en carbonnes avaient été testées (exprimées en nombre d'unités pour 100 :  
# $x_i$ : 0.3,0.7,1,1.1,1.2,2.5,2.7,3,5,5.2  
# et pour lesquelles la "charge de rupture" du nouvel alliage (en tonnes) avait été déterminée :  
# $y_i$ : 0.8,1.77,2.13,2.42,2.87,5.34,5.5,6.41,10.31,10.47 

# In[ ]:


TX = np.array([0.3,0.7,1,1.1,1.2,2.5,2.7,3,5,5.2])
TY = np.array([0.8,1.77,2.13,2.42,2.87,5.34,5.5,6.41,10.31,10.47])


# ceci est un exemple

# Il est possible de tracer les histogrammes des fréquences, des fréquences relatives et des fréquences cumulées de chacun de ces tableaux mais aussi de représenter une "boite à moustaches" qui permet de visualiser notamment la médiane, le premier et le troisième quartile d'une série statistique.

# In[ ]:


plt.hist(TX,color="RoyalBlue") # par défaut les fréquences sont réparties en 10 classes


# In[ ]:


plt.hist(TX,4,normed=True,color="PowderBlue") # fréq. relatives, réparties en 4 classes (bornes sont dans le 2nd tableau donné en sortie)
# ATTENTION : les calculs sont faits de telle façon que la somme des aires soit égale à 1


# In[ ]:


plt.hist(TX,4,normed=True,cumulative=True) # histogramme des fréquences cumulées réparties en 4 classes.


# In[ ]:


plt.boxplot(TX)


# La boite à moustache représente par un trait rouge horizontale la **médiane** de la série statistique tandis que les bornes inférieures et supérieures de la boite représentent respectivement le **premier et le troisième quartile** de cette même série.  
# On appelle *moustaches* les traits en pointillés terminés par un court segment horizontale. Ils indiquent l'amplitude de la série en indiquant son minimum et son maximum sous réserve qu'ils soient contenus entre deux bornes associés par défaut à 1.5 fois l'intervalle inter-quartile. Par exemple, si $Q_3$ désigne le troisième quartile, alors la borne supérieure des moustaches représente $(Q_3+1.5(Q_3-Q_1)$.  
# Toute valeur au-delà de ces bornes sera représentée individuellement par des disques au delà des moustaches. Elles sont considérées comme des valeurs exceptionnelles ou atypiques et apparaissent de façon distincte pour aider la discussion de l'amplitude de la série.   
# On notera qu'il n'existe pas de telles valeurs sur la série des teneurs en carbone.  
# 
# *Note :*  Ce type de représentation est particulièrement utile pour comparer deux séries de données.  
# Imaginons qu'on cherche à comparer la charge de rupture en fonction de différents alliage avec la variation de charge de rupture pour un alliage dit *témoin*, valeurs enregistrées dans un tableau `TY_temoin`. On écrira :

# In[ ]:


TY_temoin = np.array([4.8,3.77,1.24,3.56,2.87,5.04,4.5,6.13,3.81,4.74])
plt.boxplot([TY,TY_temoin])


# **Que pouvons-nous conclure ?** On commence par noter que les médianes sont sensiblement égales, proche de $4$.  
# Assez naturellement, l'amplitude de la série de données associées à l'alliage témoin est bien moindre.     
# Il semble par ailleurs que $1.24$ puisse être considérée comme une valeur aberrante (erreur expérimentale, de lecture, etc.). Il s'agira d'en discuter en observant les raisons de ce résultat mais, dans la discution comparée de ces deux séries, on retiendra une amplitude du témoin comprise entre $2.87$ et $6.13$.

# **Etude de régression linéaire simple :** Nous considérons désormais que les différents alliages représentés par `TX` constitue la *variable explicative* qui pilote le phénomène *charge de rupture* représenté par `TY`, la *variable expliquée*.  
# Il est clair que `TX` et `TY` ne sont pas indépendantes et il est naturel de vouloir estimer une relation entre la teneur en carbone et la charge de rupture.  C'est même l'enjeu de notre protocole expérimental !  
# Pour ça, on commence par tracer un nuage des points $M_i=(x_i,y_i)$  :

# In[ ]:


plt.plot(TX,TY,'ro')


# La forme "allongée" du nuage de points suggère une relation linéaire entre les deux séries.  
# On détermine donc les coefficients de la droite de régression comme indiqué au paragraphe 3. avant de tracer la droite sur la représentation graphique précédente :

# In[ ]:


a,b = np.polyfit(TX,TY,1)
print(a,b)
plt.plot(TX,TY,'ro',TX,a*TX+b,'b--')


# La preuve graphique d'une relation linéaire de la forme $y_i=1.98 x_i+0.3$ est par ailleurs confirmée par la donnée du coefficient de corrélation :

# In[ ]:


print('le coefficient de corrélation vaut : ',np.corrcoef(TX,TY)[0,1])


# ** Les courbes avec barres d'erreurs : **

# In[ ]:


print(0.2*np.random.rand(10))


# In[ ]:


x = np.linspace(0, 3, 10) # abscisses : 10 pts dans [0, 3]
y = np.exp(-x) # ordonnées
e1 = 0.2*np.random.rand(len(x)) # erreurs au hasard dans ]0,0.2[
e2 = 0.2*np.random.rand(len(x)) # idem
# Taille de l'image à afficher :
fig = plt.figure(figsize=(10,8))
# barre d"erreur Y symétrique
plt.axis([0, 3, 0, 1.2])
plt.subplot(311)
plt.errorbar(x, y, yerr=e1, fmt="go-", ecolor="b")
# barre d"erreur X symétrique
plt.axis([0, 3, 0, 1.2])
plt.subplot(312)
plt.errorbar(x, y, xerr=e1, fmt="go-", ecolor="b")
# barre d"erreur Y asymétrique
plt.axis([0, 3, 0, 1.2])
plt.subplot(313)
plt.errorbar(x, y, yerr=[e1,e2], fmt="go-", ecolor="b")


# ** Tracé d'un camembert : **

# In[ ]:


# distribution des fréquences des voyelles dans un texte rédigé en français
# commprenant 764 voyelles.
x = [115, 319, 120, 96, 108, 6]
plt.figure(figsize=(5,5))
L = ["A", "E","I", "O", "U", "Y"]
plt.pie(x, labels=L)
plt.show()
# Labels dans les 'portions':
plt.figure(figsize=(5,5))
plt.pie(x, labels=L, autopct="%.2f %%")
plt.show()
# Explosion de certaines portions:
e = [0.1, 0, 0.2, 0, 0.1, 0]
plt.figure(figsize=(5,5))
plt.pie(x, labels=L, autopct="%.2f %%", explode=e)
plt.show()


# ## 4. Le module random

# Cette bibliothèque, associée à `scipy.stats` est particulièrement utile pour toutes les questions qui concernent les dénombrements et les probabilités (cf. paragraphe 5.).
# 
# Commençons importer la bibiothèque `random` :

# In[ ]:


import random as rdm


# La fonction la plus utile est, de loin, la fonction `rdm.random()` qui simule un nombre réel pris dans l'intervalle $[0,1[$ selon une loi uniforme.

# In[ ]:


print(rdm.random())
print([rdm.random() for k in range(10)]) # les appels successifs seront supposés indépendants.


# *Remarque :* On peut vérifie en faisant 10000 appels successifs à cette fonction que les fréquences relatives de chaque intervalle de longueur $0.1$ sont sensiblement identiques. 

# In[ ]:


plt.hist([rdm.random() for k in range(10000)],normed=True)


# Il est également utile, notamment en dénombrement, de savoir simuler le tirage au hasard dans une liste d'entiers compris entre $a$ et $b$ (par exemple un tirage d'une boule dans une urne composée de boules numérotées entre $a$ et $b$).  
# On fera alors :

# In[ ]:


a,b = 1,6
print(rdm.randint(a,b)) # peut être interprété comme la simulation d'un lancer de dé.
print([rdm.randint(a,b) for k in range(10)]) # simulation de 10 lancers indépendants d'un dé équilibré


# ** Attention ** Noter que, exceptionnellement dans le cas de la fonction `randint(a,b)` les valeurs de $a$ et de $b$ sont incluses, autrement dit, cette fonction retourne un entier dans l'ensemble $\{a,a+1,\cdots,b-1,b\}$.

# **Application :** Une urne est composée de boules numérotées entre $1$ et $10$. On effectue deux tirages avec remise dans cette urne. Estimer la probabilité que la somme des numéros obtenus soit égale à $5$ 

# In[ ]:


m = 1000 # on va recommencer m fois l'expérience
LS = [] # initialisation de la liste des sommes obtenues à chacune de ces expériences
for k in range(m):
    S = 0
    for tirage in range(2):
        S= S+rdm.randint(1,10)
    LS.append(S) # on ajoute à la liste LS toute nouvelle somme de deux numéros obtenus.
print('la proportion de somme égale à 5 vaut : ',LS.count(5)/m)
print("la probabilité d'obtenir 5 vaut quant à elle : ",str(1/25))


# Pourquoi une probabilité de $1/25$ ?  
# Il suffit d'écrire que :  
# $\mathbb{P}(X+Y=5)=\mathbb{P}((X=1)\cap(Y=4))+\mathbb{P}((X=2)\cap(Y=3))+\mathbb{P}((X=3)\cap(Y=2))+\mathbb{P}((X=4)\cap(Y=1))$  
# Or les tirages sont indépendants, donc $\mathbb{P}((X=i)\cap(Y=j))=\mathbb{P}(X=i)\cdot \mathbb{P}(Y=j)=\cfrac{1}{10}\cdot\cfrac{1}{10}=\cfrac{1}{100}$.  
# *Conclusion :* $\mathbb{P}(X+Y=5)=\cfrac{4}{100}=\cfrac{1}{25}$...

# **Important !** On peut aussi définir un ensemble (sous forme de liste) et demander de choisir un élément de cet ensemble.

# In[ ]:


couleurs = ['coeur','pique','trèfle','carreau']
print(rdm.choice(couleurs))
print([rdm.choice(couleurs) for k in range(10)])


# **Application : **  On dispose d'un jeu de $32$ cartes dans lequel on prélève 1 carte.  
# Estimer la probabilité d'obtenir un Roi de coeur, d'obtenir un coeur :

# In[ ]:


hauteurs = [1,7,8,9,10,'Valet','Dame','Roi']
couleurs = ['coeur','pique','trèfle','carreau']
univers=[str(v)+' '+c for v in hauteurs for c in couleurs]
m = 10000
nbeAsP = 0
nbCoeurs = 0
for n in range(m):
    if rdm.choice(univers) == '1 pique':
        nbAsP = nbAsP + 1
    if rdm.choice(couleurs) == 'coeur':
        nbCoeurs += 1

print("la probabilité d'obenir un As de pique est estimée à : ",nbAsP/m)
print('La probabilité exacte vaut : ',1/32)
print("la probabilité d'obenir un coeur est estimée à : ",nbCoeurs/m)
print('La probabilité exacte vaut : ',1/4)


# ## 5. Le module scipy.stats

# Python permet de travailler simplement sur l'ensemble des lois de probabilité au programme des deux années de BCPST.  
# On citera pour mémoire les lois discrètes qui sont vues en première année :
# * La loi uniforme
# * La loi de bernoulli
# * La loi binomiale
# * La loi hypergéométrique  
# 
# Si $X$ est une variable aléatoire qui suit l'une de ces lois, alors on doit savoir :  
# 1) Déterminer $\mathbb{P}(X=k)$, $\forall k\in X(\Omega)$.  
# 2) Expliciter et tracer la fonction de répartition définie par : $F_X(x)=\mathbb{P}(X\leq x)$, $\forall x\in\mathbb{R}$. 
# 3) Connaître l'espérance et la variance (sauf pour la loi hypergéométrique).  
# 
# Il suffit, pour chacun de ces trois points, d'importer la loi qui nous intéresse grâce à :  
# `from scipy.stats import nom_de_la_loi`  
# Alors :  
# 1) `nom_de_la_loi.pmf(k)` pour *Probability Mass Function* permet de calculer $\mathbb{P}(X=k)$.  
# 2) `nom_de_la_loi.cdf(x)` pour *Cumulative Density Function* permet de calculer $F_X(x)$.  
# 3) `nom_de_la_loi.stats(..,moments='mv')` retourne la `m`-oyenne et la `v`-ariance.  
# 
# On lira avec profit la documentation associée au module [`scipy.stats`](http://docs.scipy.org/doc/scipy/reference/stats.html) en insistant notamment sur la rubrique *Discrete distributions* de la table des matières. 

# **Exemple : ** La loi binomiale de paramètres $n=5$ et $p=0.4$.  
# Déterminons $\mathbb{P}(X=k)$ pour $k\in \{0,\cdots,5\}$ ainsi que sa loi de probabilité :

# In[ ]:


from scipy.stats import binom
n,p = 5,0.4
k = 3
print('P(X = ',k,') = ',binom.pmf(k,n,p)) # attention à l'ordre des arguments (cf. la documentation)
# n'hésitez pas à essayer avec d'autres valeurs de k puis d'autres valeurs de n et de p !


# Donner une représentation graphique de la loi de probabilité de $X\hookrightarrow\mathcal{B}(n,p)$ :

# In[ ]:


n,p = 5,0.4
Xom = range(n+1) # X(Omega)={0,...,n}
LoiX = [binom.pmf(k,n,p) for k in Xom]
plt.plot(Xom,LoiX, 'ro')
plt.vlines(Xom,0,LoiX, colors='k',lw=2) # faire help(plt.vlines) pour l'emploi de cette primitive.
plt.xlim(xmin=-0.5,xmax=n+0.5) # Pas indispensable... juste si on veut modifier le tracer des bornes des abscisses


# Donner une représentation graphique de la fonction de répartition de $X\hookrightarrow\mathcal{B}(n,p)$ :

# In[ ]:


n,p = 5,0.4
Xom = range(n+1) # X(Omega)={0,...,n}
largeur = 1
plt.bar(Xom,binom.cdf(Xom,n,p),width = largeur,color='r',alpha = 0.6)
plt.title('Fonction de répartition de la loi binomiale de paramètres '+str(n)+' et '+str(p))


# Enfin pour l'**espérance** et la **variance** de la loi binomiale, on écrira :

# In[ ]:


n,p = 5,0.4
binom.stats(n,p,moments='mv')


# On vérifie qu'on a bien $\mathbb{E}(X)=np=5\times 0.4=2$ et $\mathbb{V}(X)=npq=5\times 0.4\times 0.6=1.2$

# ** A vous de vous entraîner avec d'autres paramètres ou avec les autres lois au programme !**

# ## 6. Le module IPython.display

# Anecdotique mais peut être utile si on souhaite importer ou lire des pages HTML ou encore des vidéos.  
# Par exemple, pour introduction au module `Turtle` :

# In[ ]:


from IPython.display import YouTubeVideo
# Pour une petite introduction au module `Turtle`.
YouTubeVideo('fmyeDwQeyQY')


# In[ ]: