Calcul Scientifique TP3 Approximation polynomiale corrig?

October 25, 2019

1 TP3 : Approximation polynomiale

1.1 Objectif
L’objectif de ce TP 3 est de :
- Observer les données à modéliser, - Choisir le modèle adéquat, ajustant ces données, - Estimer
les paramètres du modèle choisi au sens des moindres carrées.

1.2 Rappel du cours

1.2.1 Ajustement de données
On considère (n + 1) points, ( xi , yi )0≤i≤n , d’abscisses deux à deux distinctes.
On note par X = ( xi )0≤i≤n ∈ M(n+1),1 (R) et Y = (yi )0≤i≤n ∈ M(n+1),1 (R).
- La problématique consiste à juster un modèle défini par une fonction f (.; Λ) paramétrée par Λ =
(λ j )0≤ j≤ p ∈ M( p+1),1 (R) aux données ( xi , yi )0≤i≤n : approcher f ( xi ; Λ) par yi , ∀i ∈ {0, · · · , n}.
- La fonction f : Rn+1 → Rn+1 pourrait être linéaire ou non-linéaire. - Le choix du modèle
dépendra du phénomène étudié.
- Nous nous focalisons dans ce TP sur les fonctions polynomiales :
p
f (t; Λ) = ∑ λj tj, t ∈ R.
j =0

- Deux cas possibles qui se présentent : - p = n : il s’agit alors d’une interpolation polynomiale -
p < n : il s’agit d’une approximation polynomiale

1.2.2 Estimation du vecteur Λ au sens des moindres carrées

Le vecteur Λ∗ , optimal au sens des moindres carrées, pour lequel on a un meilleur ajustement
du modèle polynomiale $ f(t;Λ) = ∑ _{ j = 0} pλ_jtj, ; t ∈ R, $auxdonnes( xi , yi )0≤i≤n , est celui
minimisant la fonctionnelle
n
F (Λ; X ) = ∑ ( f ( x i ; Λ ) − y i )2 .
i =0

Celà revient à dire que :

Λ∗ = argmin F (Λ; X ).
Λ ∈ R p +1

1
 En effet, on note l’erreur résiduelle ei = f ( xi ; Λ) − yi , ∀i ∈ {0, · · · , n}, l’erreur commise par le
modèle, F (Λ; X ) sera ainsi exprimée comme suit :
n ( )2 n
∑ ∑ ei2 =∥ ε ∥2
2
F (Λ; X ) = f ( xi ; Λ ) − yi =
i =0 i =0

avec ε = (ei )0≤ j≤n , et ∥ . ∥2 désigne la norme euclidienne.

p
∑ λ j xi − yi , ∀i
j
En développeant l’expression de l’erreur ei = ∈ {0, · · · , n}, : un système de
j =0
n + 1 équations linéaires sera généré selon la forme matricielle suivante :
   p    
e0 1 x0 x02 · · · x0 λ0 y0
 ..   .. ..  ..   .. 
 .  . .    
    .   . 
 ei  = 1 x i x 2 · · · x p   λ2  −  y i 
   i i    
 ..   .. ..  ..   .. 
 .  . .   .  .
p
en 1 xn xn · · · xn
2 λp yn
| {z } | {z }| {z } | {z }
ε A Λ Y
2 2
$⇒ $F (Λ; X ) =∥ ε ∥2 =∥ AΛ − Y ∥2 .
La vecteur Λ∗ minimisant F (Λ; X ) satisfait la condition suivante :

∂F (Λ∗ ; X )
▽ F ( Λ ∗ ; X ) = 0 p +1 : = 0, 0 ≤ i ≤ p.
∂λi

où 0 p+1 désigne le vecteur nul de R p+1 , ▽ F (Λ∗ ; X ) correspond au gradient de F (.; X ) appliqué à
∂F (Λ∗ ; X )
Λ∗ et présente la dérivée partielle de F (., X ) par rapport à λi , appliquée à Λ∗ .
∂λi
En développant l’expression des dérivées partielles, il s’en suit :

▽ F (Λ∗ ; X ) = 0 p+1 ⇒ 2 t A( AΛ∗ − Y ) = 0 p+1

⇒ Λ∗ = (t AA)−1 t
AY
Cas de p=1
 n n n n

 ∑ ∑ yi − ∑ xi ∑ xi yi  xi2
 i =0 i =0 
 i =0 i =0 
 n ( n ) 
 ∗  
2
 ( n + 1) ∑ x i − ∑ x i
2

λ0  
 i = 0 i = 0 
Λ∗ =   =  
 
λ1∗  n n n 
 − ∑ x i ∑ y i + ( n + 1) ∑ x i y i 
 
 i =0 i =0 i =0 
 n ( n ) 
 2 
(n + 1) ∑ xi2 − ∑ xi
i =0 i =0

2
 2 Programmation
Ecrire une fonction, que nous appelons LSA pour ‘Least Squared Approximation’ qui, pour un
jeux de données (X,Y) et un ordre de modèle polynomiale p, retourne le vecteur Λ∗ de coefficients
estimés du modèle sous-jacent, au sens des moindres carrées.
[1]: import numpy as np
import [Link] as plt
[2]: def LSA(X,Y,p):
n=len(X)
A=[Link]((n,p+1))
for i in [Link](1,p+1):
A[:,i]=X[:,0]**i
Lambda=[Link]([Link]().dot(A)).dot([Link]().dot(Y))
return Lambda
Ecrire une fonction Poly(t,coeff) qui évalue un polynôme de coefficients coeffen t.
[3]: def Poly(t,coeff):
[Link](len(t),1)
n=len(coeff)
P=coeff[0]*[Link]((len(t),1))
for i in [Link](1,n):
P+=coeff[i]*t**i
return P

2.1 Exercice
1. Déterminer l’expression de la droite optimale, au sens des moindres carrées, qui ajuste les
points suivants :

i 0 1 2 3 4
xi 1 2 3 4 5
yi 0.9 1.5 3.5 4.2 4.9

2. Représenter ensuite graphiquement le résultat de l’ajustement.

[4]: X=[Link]([[1,2,3,4,5]]).transpose()
Y=[Link]([[0.9,1.5,3.5,4.2,4.9]]).transpose()
Lambda=LSA(X,Y,1)
[5]: t=[Link](-1,6,.01)
t=[Link](-1,1)
[Link](figsize=(20,10))
[Link](X,Y,'ro',t,Poly(t,Lambda),'b--',linewidth=3,markersize=12)
[Link]('t',fontsize=30)
[Link](fontsize=20)
[Link](fontsize=20)

3
 [Link](('Mesures','Approximation'),fontsize=30, loc = 0)
[Link](True)

3 Application 1
1. On considère les données générées suivantes :
n=100
[Link](1)
noise=[Link](0,0.2,100)
X1=[Link](-1,1,n)
Y1=X1+noise
Y2=X1**2+noise
Y3=X**3+noise

2. Tracer les graphes de Yi, i=0,1,2, en fonction de X1.

3. Observer les graphiques et justifier pour chaque cas le choix du modèle ajustant les données.

4. Estimer les paramètres des modèles proposés. Comparer les aux valeurs avec lesquelles les
nous avons généré les données, en terme d’erreur quadratique.

5. Changer le nombre de points n=10ˆ5. Observer les nouveaux résultats et conclure.

Question 1
[6]: n=100
[Link](1)
noise=[Link](0,0.2,n)

4
 noise=[Link](-1,1)
X1=[Link](-1,1,n)
X1=[Link](-1,1)
Y1=X1+noise
Y2=X1**2+noise
Y3=X1**3+noise
Question 2
[7]: # Graphique 1
[Link](figsize=(20,10))
[Link](X1,Y1,'ro',linewidth=3,markersize=12)
[Link]('X1',fontsize=30)
[Link]('Y1',fontsize=30)
[Link](fontsize=20)
[Link](fontsize=20)
[Link](True)

[8]: # Graphique 2
[Link](figsize=(20,10))
[Link](X1,Y2,'ro',linewidth=3,markersize=12)
[Link]('X1',fontsize=30)
[Link]('Y2',fontsize=30)
[Link](fontsize=20)
[Link](fontsize=20)
[Link](True)

5
[9]: #Graphique 3
[Link](figsize=(20,10))
[Link](X1,Y3,'ro',linewidth=3,markersize=12)
[Link](True)
[Link]('X1',fontsize=30)
[Link]('Y3',fontsize=30)
[Link](fontsize=20)
[Link](fontsize=20)
[9]: (array([-1.5, -1. , -0.5, 0. , 0.5, 1. , 1.5]),
<a list of 7 Text yticklabel objects>)

6
 Question 4
[10]: Lambda1=LSA(X1,Y1,1)
Lambda2=LSA(X1,Y2,2)
Lambda3=LSA(X1,Y3,3)

# Les coefficients des modèles polynomiaux avec lesquels les données ont été␣
,→générées

Coeff_Y1=[Link]([[0,1]]).transpose()
Coeff_Y2=[Link]([[0,0,1]]).transpose()
Coeff_Y3=[Link]([[0,0,0,1]]).transpose()

E1=[Link](Lambda1-Coeff_Y1,2)**2
E2=[Link](Lambda2-Coeff_Y2,2)**2
E3=[Link](Lambda3-Coeff_Y3,2)**2

print(E1,E2,E3)

0.0012413931877733055 0.0034087916857334587 0.017180574118298892

Question 5
[11]: n=10**5
[Link](1)
noise=[Link](0,0.2,n)
noise=[Link](-1,1)
X1=[Link](-1,1,n)
X1=[Link](-1,1)
Y1=X1+noise
Y2=X1**2+noise
Y3=X1**3+noise

Lambda1=LSA(X1,Y1,1)
Lambda2=LSA(X1,Y2,2)
Lambda3=LSA(X1,Y3,3)

# Les coefficients des modèles polynomiaux avec lesquels les données ont été␣
,→générées

Coeff_Y1=[Link]([[0,1]]).transpose()
Coeff_Y2=[Link]([[0,0,1]]).transpose()
Coeff_Y3=[Link]([[0,0,0,1]]).transpose()

E1=[Link](Lambda1-Coeff_Y1,2)**2
E2=[Link](Lambda2-Coeff_Y2,2)**2
E3=[Link](Lambda3-Coeff_Y3,2)**2

print(E1,E2,E3)

7
 1.9305760506954906e-06 2.551335611566981e-05 4.44488540921642e-05

4 Application 2

[12]: import pandas as pd # importation du module pandas

[13]: # lecture des données
boston = pd.read_csv('[Link]')

4.1 Notes
Description des données boston

:Number of Instances: 506

:Number of Attributes: 13 numeric/categorical predictive

:Median Value (attribute 14) is usually the target

:Attribute Information (in order):

- CRIM per capita crime rate by town
- ZN proportion of residential land zoned for lots over 25,000 [Link].
- INDUS proportion of non-retail business acres per town
- CHAS Charles River dummy variable (= 1 if tract bounds river; 0 otherwise)
- NOX nitric oxides concentration (parts per 10 million)
- RM average number of rooms per dwelling
- AGE proportion of owner-occupied units built prior to 1940
- DIS weighted distances to five Boston employment centres
- RAD index of accessibility to radial highways
- TAX full-value property-tax rate per $10,000
- PTRATIO pupil-teacher ratio by town
- B 1000(Bk - 0.63)^2 where Bk is the proportion of blacks by town
- LSTAT % lower status of the population
- MEDV Median value of owner-occupied homes in $1000's

Exploration des objets de type dataframe

1. Executer et commenter les instructions suivantes :
[14]: print(type(boston))
#boston
#boston[0:10]
#[Link]()
#[Link]()
#boston['CRIM']
#boston['CRIM'].values
#[Link]()
#Z=[Link]('CRIM',axis=1)

8
 <class '[Link]'>

2. Observer l’évolution des valeurs des maisons Y en fonction du nombre de chambres

Xqu’elles contiennent

[15]: X=boston['RM'].values
X=[Link](-1,1)
Y=boston['MEDV'].values
Y=[Link](-1,1)
[Link](figsize=(20,10))
[Link](X,Y,'ro',markersize=8)
[Link]('Value of House /1000 ($)',fontsize=30)
[Link]('Number of rooms',fontsize=30)
[Link]()

3. Approcher les données (X,Y) par un modèle polynomiale. Donner son ordre.

[16]: [Link](figsize=(20,10))
[Link](X,Y,'ro',X, Poly(X,LSA(X,Y,1)),'b--', linewidth=3,markersize=8)
[Link]('Value of House /1000 ($)',fontsize=30)
[Link]('Number of rooms',fontsize=30)
[Link](fontsize=20)
[Link](fontsize=20)
[16]: (array([-10., 0., 10., 20., 30., 40., 50., 60.]),
<a list of 8 Text yticklabel objects>)

9
 4. Procéder de la sorte pour un modèle polynomiale de dégré 2

[17]: [Link](figsize=(20,10))
[Link](X,Y,'ro',X, Poly(X,LSA(X,Y,2)),'b^', linewidth=3,markersize=8)
[Link]('Value of House /1000 ($)',fontsize=30)
[Link]('Number of rooms',fontsize=30)
[Link](fontsize=20)
[Link](fontsize=20)
[17]: (array([ 0., 10., 20., 30., 40., 50., 60., 70.]),
<a list of 8 Text yticklabel objects>)

10
 5. Comment peut on choisir entre les modèles polynomiaux?
6. Ecrire un programme qui trace l’erreur quadratique commise par l’ajustement en fonction
de p, le degré du prolynôme d’approximation

[18]: p=5
Err=[Link]((p,1))
for i in [Link](0,p):
Err[i]=[Link](Poly(X,LSA(X,Y,i))-Y,2)
[19]: [Link](figsize=(20,10))
[Link](Err,'ro', linewidth=3,markersize=8)
[Link]('Ordre du polynôme',fontsize=30)
[Link]('Erreur d\'ajustement',fontsize=30)
#[Link]('log')
[Link](fontsize=20)
[Link](fontsize=20)
[19]: (array([130., 140., 150., 160., 170., 180., 190., 200., 210., 220.]),
<a list of 10 Text yticklabel objects>)

11
 5 Application 3

[20]: # lecture des données

mpg = pd.read_csv('[Link]')

5.1 Notes
Description des données mpg
This dataset is a slightly modified version of the dataset provided in the StatLib library. In line
with the use by Ross Quinlan (1993) in predicting the attribute “mpg”, 8 of the original instances
were removed because they had unknown values for the “mpg” attribute. The original dataset is
available in the file “[Link]-original”.
“The data concerns city-cycle fuel consumption in miles per gallon, to be predicted in terms of
3 multivalued discrete and 5 continuous attributes.” (Quinlan, 1993)
Attribute Information:

1. mpg: continuous
2. cylinders: multi-valued discrete
3. displacement: continuous
4. horsepower: continuous
5. weight: continuous
6. acceleration: continuous
7. model year: multi-valued discrete
8. origin: multi-valued discrete
9. car name: string (unique for each instance)

Quinlan,R. (1993). Combining Instance-Based and Model-Based Learning. In Proceedings on

the Tenth International Conference of Machine Learning, 236-243, University of Massachusetts,
Amherst. Morgan Kaufmann

1. Explorer les données MPG

[21]: [Link]()
[Link](20)

<class '[Link]'>
RangeIndex: 406 entries, 0 to 405
Data columns (total 10 columns):
Unnamed: 0 406 non-null int64
mpg 398 non-null float64
cylinders 406 non-null float64
displacement 406 non-null float64
horsepower 400 non-null float64
weight 406 non-null float64
acceleration 406 non-null float64
model 406 non-null float64
origin 406 non-null float64

12
 car_name 406 non-null object
dtypes: float64(8), int64(1), object(1)
memory usage: 31.8+ KB

[21]: Unnamed: 0 mpg cylinders displacement horsepower weight \

0 0 18.0 8.0 307.0 130.0 3504.0
1 1 15.0 8.0 350.0 165.0 3693.0
2 2 18.0 8.0 318.0 150.0 3436.0
3 3 16.0 8.0 304.0 150.0 3433.0
4 4 17.0 8.0 302.0 140.0 3449.0
5 5 15.0 8.0 429.0 198.0 4341.0
6 6 14.0 8.0 454.0 220.0 4354.0
7 7 14.0 8.0 440.0 215.0 4312.0
8 8 14.0 8.0 455.0 225.0 4425.0
9 9 15.0 8.0 390.0 190.0 3850.0
10 10 NaN 4.0 133.0 115.0 3090.0
11 11 NaN 8.0 350.0 165.0 4142.0
12 12 NaN 8.0 351.0 153.0 4034.0
13 13 NaN 8.0 383.0 175.0 4166.0
14 14 NaN 8.0 360.0 175.0 3850.0
15 15 15.0 8.0 383.0 170.0 3563.0
16 16 14.0 8.0 340.0 160.0 3609.0
17 17 NaN 8.0 302.0 140.0 3353.0
18 18 15.0 8.0 400.0 150.0 3761.0
19 19 14.0 8.0 455.0 225.0 3086.0

acceleration model origin car_name

0 12.0 70.0 1.0 chevrolet chevelle malibu
1 11.5 70.0 1.0 buick skylark 320
2 11.0 70.0 1.0 plymouth satellite
3 12.0 70.0 1.0 amc rebel sst
4 10.5 70.0 1.0 ford torino
5 10.0 70.0 1.0 ford galaxie 500
6 9.0 70.0 1.0 chevrolet impala
7 8.5 70.0 1.0 plymouth fury iii
8 10.0 70.0 1.0 pontiac catalina
9 8.5 70.0 1.0 amc ambassador dpl
10 17.5 70.0 2.0 citroen ds-21 pallas
11 11.5 70.0 1.0 chevrolet chevelle concours (sw)
12 11.0 70.0 1.0 ford torino (sw)
13 10.5 70.0 1.0 plymouth satellite (sw)
14 11.0 70.0 1.0 amc rebel sst (sw)
15 10.0 70.0 1.0 dodge challenger se
16 8.0 70.0 1.0 plymouth 'cuda 340
17 8.0 70.0 1.0 ford mustang boss 302
18 9.5 70.0 1.0 chevrolet monte carlo
19 10.0 70.0 1.0 buick estate wagon (sw)

13
 2. Observer l’évolution de la consommation du carburant en fonction des poids des véhicules.

[22]: Weight=mpg['weight'].values
Weight=[Link](-1,1)
MPG=mpg['mpg'].values
MPG=[Link](-1,1)

[Link](figsize=(20,10))
[Link](Weight,MPG,'ro',linewidth=3,markersize=8)
[Link]('Weight',fontsize=30)
[Link]('MPG',fontsize=30)
[Link](fontsize=20)
[Link](fontsize=20)
[22]: (array([ 5., 10., 15., 20., 25., 30., 35., 40., 45., 50.]),
<a list of 10 Text yticklabel objects>)

3. Modéliser la relation MPG vs Weight par un modèle polynomiale. Argumenter votre choix.
4. Donner une estimation du poids d’un véhicule qui consomme 30 MPG. Argumenter votre
réponse.

[23]: # Nettoyage des données

L=[Link]([Link](MPG)==True)[0]
MPG=[Link](MPG,L)
MPG=[Link](-1,1)
Weight=[Link](Weight,L)
Weight=[Link](-1,1)

14
[24]: [Link](figsize=(20,10))
[Link](Weight,MPG,'ro',Weight, Poly(Weight,LSA(Weight,MPG,2)),'b*',␣
,→linewidth=3,markersize=8)

[Link]('Weight',fontsize=30)
[Link]('MPG',fontsize=30)
[Link](fontsize=20)
[Link](fontsize=20)
[24]: (array([ 5., 10., 15., 20., 25., 30., 35., 40., 45., 50.]),
<a list of 10 Text yticklabel objects>)

[25]: x=[Link](1500,6001)
[Link]
x=[Link](-1,1)
y=Poly(x,LSA(Weight,MPG,2))
y[[Link](y<=30)[0][0]]
x[[Link](y<=30)[0][0]] # Le poids du véhicule estimé
[25]: array([2184])

6 Références
[1] Kiusalaas, J. (2013). Numerical methods in engineering with Python 3. Cambridge university
press.
[2] Numpy Package
[3] Mathplotlib Package
[4] Jupyter markdowns
[5] Pandas Package

15