TP5-05DEC - Regression non linéaire

Encadré par Solym MANOU-ABI

Exercice 1 – Simulation et estimation d’une fonction de régression

Dans ce premier exercice, on souhaite étudier une relation de la forme \(Y_i = f(x_i) + \varepsilon_i\) en simulant des observations pour lesquelles la variance du bruit dépend de la valeur de \(x\). On prendra \(f(x) = x\) et un bruit de variance croissante via un facteur \(w_i = 1 + x_i/2\). Les instructions sont données dans le texte de l’énoncé. On suit la démarche pas à pas et on commente les résultats obtenus.

Génération des données

On commence par créer un jeu de données contenant les variables x, w et Y. La variable x est une suite allant de 1 à 30 par pas de 0,5. On fixe un graine (set.seed) pour rendre les résultats reproductibles. On définit ensuite w = 1 + x/2, on génère un bruit Z gaussien standard et on calcule Y = x + w * Z. Enfin, on rassemble tout dans un data.frame tab.

set.seed(42)                      # graine pour la reproductibilité
x <- seq(1, 30, by = 0.5)         # points de 1 à 30 avec un pas de 0,5
w <- 1 + x/2                      # facteur de variance croissante
Z <- rnorm(length(x))            # bruit gaussien de moyenne 0 et variance 1
Y <- x + w * Z                    # modèle Y = f(x) + w * Z avec f(x)=x
tab <- data.frame(x = x, w = w, Y = Y)
head(tab)

Observation : on observe que la variable w augmente avec x, ce qui signifie que l’ampleur du bruit est faible pour les petites valeurs de x et plus forte pour les grandes valeurs. Autrement dit, on crée intentionnellement une hétéroscédasticité, c’est‑à‑dire une variance non constante de Y en fonction de x.

Nuage de points (x, Y)

Nous représentons ensuite le nuage de points \((x_i, Y_i)\) pour visualiser la relation entre la variable explicative x et la variable réponse Y. On s’attend à voir un alignement approximatif le long de la droite de pente 1, mais avec une dispersion qui augmente lorsque x croît.

plot(tab$x, tab$Y,
     pch = 16, cex = 0.7, col = "steelblue",
     xlab = "x", ylab = "Y",
     main = "Exercice 1 – Nuage de points (x, Y)")
grid()

Observation : le nuage de points s’élargit clairement lorsque x augmente : les points situés à droite sont beaucoup plus dispersés que ceux situés à gauche. Cela reflète bien la variance hétérogène introduite par w. On voit également que la tendance centrale suit une droite approximative de pente 1.

Estimation par moindres carrés pondérés

Pour estimer la fonction de régression \(f\), on peut effectuer une régression linéaire pondérée. Comme on sait que la variance conditionnelle de Y vaut \(\mathrm{Var}(Y_i\mid x_i) = w_i^2\), un choix naturel est d’utiliser des poids égaux à \(1 / w_i^2\), c’est‑à‑dire l’inverse de la variance. De cette manière, les observations les moins bruitées (petit w) ont plus d’influence que les observations très bruitées.

# Ajustement par moindres carrés pondérés
reg1 <- lm(Y ~ x, data = tab, weights = 1/(w^2))
summary(reg1)

Call:
lm(formula = Y ~ x, data = tab, weights = 1/(w^2))

Weighted Residuals:
    Min      1Q  Median      3Q     Max 
-2.7710 -0.5259 -0.0963  0.7512  2.1938 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   1.4869     0.9604   1.548    0.127    
x             0.8320     0.1207   6.891 4.83e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.145 on 57 degrees of freedom
Multiple R-squared:  0.4545,    Adjusted R-squared:  0.4449 
F-statistic: 47.49 on 1 and 57 DF,  p-value: 4.832e-09

# Tracé de la droite ajustée
plot(tab$x, tab$Y,
     pch = 16, cex = 0.7, col = "steelblue",
     xlab = "x", ylab = "Y",
     main = "Exercice 1 – MCP pondérés : droite ajustée")
grid()

abline(reg1, col = "darkorange", lwd = 2)

Observation : l’estimation linéaire pondérée affiche des coefficients proches de 1 pour la pente et une ordonnée à l’origine faible (proche de 0). La droite ajustée suit bien la tendance centrale des données. Le fait de pondérer corrige l’influence excessive des points très bruités situés à droite et donne un ajustement plus stable que si l’on utilisait une régression ordinaire.

Estimation par polynômes locaux (loess)

Une alternative non paramétrique consiste à utiliser un estimateur local par polynômes (fonction loess), qui ajuste une courbe lisse sans supposer une forme linéaire globale. On utilise ici le lissage par défaut de loess.

# Ajustement par polynômes locaux
locpoly <- loess(Y ~ x, data = tab)
fitted_values <- fitted(locpoly)

# Tracé du nuage et de la courbe lissée
plot(tab$x, tab$Y,
     pch = 16, cex = 0.7, col = "steelblue",
     xlab = "x", ylab = "Y",
     main = "Exercice 1 – Loess : courbe ajustée")
grid()

lines(tab$x, fitted_values, col = "forestgreen", lwd = 2)

Observation : la courbe issue du lissage local suit globalement la tendance linéaire mais on peut observer de petites ondulations dues au bruit. Le lissage présente l’avantage de ne pas imposer une forme particulière pour \(f\), tout en s’adaptant localement aux données. Si la fenêtre (span) était plus petite, la courbe serait plus ondulée ; si elle était plus grande, la courbe se rapprocherait d’une droite.

Exercice 2 – Estimation non paramétrique d’une série chronologique

On considère une fonction périodique \(f(t) = 0{,}5 + 0{,}4 \sin(2\pi t)\) et une série chronologique bruitée \(Y_t = f(t) + \varepsilon_t\) où \(\varepsilon_t\sim \mathcal{N}(0, 0{,}25)\). L’objectif est d’illustrer différentes méthodes de lissage et de comparer leur comportement.

Représentation de f(t) et de Y_t sur [0, 1]

On crée une grille fine t allant de 0 à 1 par pas de 0,01, on calcule la fonction f et on simule un bruit gaussien de variance 0,25 (écart‑type 0,5). On trace ensuite les courbes.

set.seed(123)
t <- seq(0, 1, by = 0.01)
f <- 0.5 + 0.4 * sin(2 * pi * t)
epsilon <- rnorm(length(t), mean = 0, sd = 0.5)
Y <- f + epsilon

plot(t, f, type = "l", lwd = 2, col = "darkgreen",
     ylim = range(c(f, Y)),
     xlab = "t", ylab = "Valeur",
     main = "Exercice 2 – Fonction f(t) et série Y_t")
lines(t, Y, col = "steelblue", lty = 2)

legend("topright", legend = c("f(t)", "Y_t"),
       col = c("darkgreen", "steelblue"), lty = c(1, 2), lwd = 2)

Observation : la fonction \(f\) est sinusoïdale et régulière. La courbe Y_t suit cette tendance périodique mais les observations sont perturbées par un bruit gaussien, d’où une dispersion autour de f(t).

Génération d’une grille T1 et d’une réalisation Y sur T1

On crée maintenant une grille T1 sur \([0, 1]\) avec un pas de 0,05 et on simule une réalisation de Y sur ces points uniquement. Cette grille plus grossière sert de base aux estimations.

T1 <- seq(0, 1, by = 0.05)
fT1 <- 0.5 + 0.4 * sin(2 * pi * T1)
set.seed(456)
Y_T1 <- fT1 + rnorm(length(T1), mean = 0, sd = 0.5)
data_ex2 <- data.frame(T1 = T1, Y = Y_T1)
head(data_ex2)

Observation : la grille T1 comporte 21 points. La série Y générée sur ces points présente des écarts parfois importants par rapport à la fonction sous‑jacente f(T1) en raison du bruit.

Lissage par polynômes locaux (locpoly)

On utilise la fonction locpoly du package KernSmooth pour estimer la fonction de régression par polynômes locaux d’ordre 1. Le paramètre bandwidth \(h\) contrôle la largeur de la fenêtre de lissage. Pour h = 0.05, on obtient l’estimation suivante :

library(KernSmooth)
h <- 0.05
reg_locpoly <- locpoly(T1, Y_T1, bandwidth = h)

plot(T1, Y_T1, pch = 16, cex = 0.7, col = "steelblue",
     xlab = "t", ylab = "Y",
     main = paste("Exercice 2 – locpoly avec h =", h))
grid()

lines(reg_locpoly, col = "darkorange", lwd = 2)
lines(T1, fT1, col = "darkgreen", lwd = 2, lty = 2)

legend("topright", legend = c("Y_t", "locpoly", "f(t)"),
       col = c("steelblue", "darkorange", "darkgreen"),
       lty = c(NA, 1, 2), pch = c(16, NA, NA), lwd = c(NA, 2, 2))

Observation : pour h = 0.05, l’estimateur locpoly suit de près la structure sinusoïdale de la fonction réelle. Les variations rapides sont bien capturées, bien que le bruit puisse induire quelques petites oscillations. Le choix d’un h trop petit rendrait l’estimation trop « courbée », tandis qu’un h trop grand lisserait exagérément la courbe.

Effet du paramètre de lissage h

On fait varier la valeur de h pour observer l’influence du lissage. Ci‑dessous, on considère h = 0.15, soit une fenêtre trois fois plus large. Cela produit un lissage plus agressif.

h2 <- 0.15
reg_locpoly2 <- locpoly(T1, Y_T1, bandwidth = h2)

plot(T1, Y_T1, pch = 16, cex = 0.7, col = "steelblue",
     xlab = "t", ylab = "Y",
     main = paste("Exercice 2 – locpoly avec h =", h2))
grid()

lines(reg_locpoly2, col = "darkorange", lwd = 2)
lines(T1, fT1, col = "darkgreen", lwd = 2, lty = 2)

legend("topright", legend = c("Y_t", "locpoly", "f(t)"),
       col = c("steelblue", "darkorange", "darkgreen"),
       lty = c(NA, 1, 2), pch = c(16, NA, NA), lwd = c(NA, 2, 2))

Observation : avec h = 0.15, l’estimation est beaucoup plus lisse. La courbe suit toujours la tendance globale, mais les sommets et creux de la sinusoïde sont atténués. Un h trop élevé peut masquer des structures réelles des données, tandis qu’un h trop faible peut sur‑ajuster le bruit.

Lissage à noyau (ksmooth)

La fonction ksmooth fournit un estimateur de Nadaraya–Watson basé sur un noyau (par défaut gaussien). On teste deux choix de bande de lissage : bandwidth = 0.2 et bandwidth = 0.05 (note : la valeur h choisie pour locpoly précédemment). On représente ces estimations en rouge et bleu.

library(stats)

estim1 <- ksmooth(T1, Y_T1, x.points = T1, bandwidth = 0.2, kernel = "normal")
estim2 <- ksmooth(T1, Y_T1, x.points = T1, bandwidth = 0.05, kernel = "normal")

plot(T1, Y_T1, pch = 16, cex = 0.7, col = "steelblue",
     xlab = "t", ylab = "Y",
     main = "Exercice 2 – Comparaison ksmooth")
grid()

lines(estim1$x, estim1$y, col = "red", lwd = 2)
lines(estim2$x, estim2$y, col = "blue", lwd = 2)

lines(T1, fT1, col = "darkgreen", lwd = 2, lty = 2)
legend("topright", legend = c("Y_t", "ksmooth h=0.2", "ksmooth h=0.05", "f(t)"),
       col = c("steelblue", "red", "blue", "darkgreen"),
       lty = c(NA, 1, 1, 2), pch = c(16, NA, NA, NA), lwd = c(NA, 2, 2, 2))

Observation : l’estimateur à noyau avec h = 0.2 (rouge) est plus lisse que celui avec h = 0.05 (bleu). La courbe rouge suit la tendance globale en atténuant fortement les fluctuations locales, tandis que la courbe bleue réagit davantage aux variations point par point. Le choix du paramètre bandwidth est donc crucial pour ajuster le degré de lissage.

Influence de différentes valeurs de h

On peut tester des valeurs de bande beaucoup plus petites ou plus grandes. Ci‑dessous, on compare h = 0.5 (lissage très large) et h = 0.1 (lissage modéré).

estim_h05 <- ksmooth(T1, Y_T1, x.points = T1, bandwidth = 0.5, kernel = "normal")
estim_h01 <- ksmooth(T1, Y_T1, x.points = T1, bandwidth = 0.1, kernel = "normal")

plot(T1, Y_T1, pch = 16, cex = 0.7, col = "steelblue",
     xlab = "t", ylab = "Y",
     main = "Exercice 2 – Effet de h sur ksmooth")
grid()

lines(estim_h05$x, estim_h05$y, col = "purple", lwd = 2)
lines(estim_h01$x, estim_h01$y, col = "darkorange", lwd = 2)

lines(T1, fT1, col = "darkgreen", lwd = 2, lty = 2)
legend("topright", legend = c("Y_t", "ksmooth h=0.5", "ksmooth h=0.1", "f(t)"),
       col = c("steelblue", "purple", "darkorange", "darkgreen"),
       lty = c(NA, 1, 1, 2), pch = c(16, NA, NA, NA), lwd = c(NA, 2, 2, 2))

Observation : avec h = 0.5, la courbe (violette) est très lisse et ne suit plus vraiment la forme sinusoïdale : le lissage écrase la structure périodique. Avec h = 0.1 (orange), l’ajustement est intermédiaire : on voit la sinusoïde, mais l’estimation est moins oscillante qu’avec h = 0.05.

Comparaison des qualités prédictives

Pour comparer les méthodes locpoly et ksmooth, on calcule l’erreur quadratique moyenne (EQM) entre les estimations et la fonction réelle f(T1). Une EQM plus faible indique une meilleure précision prédictive. On considère locpoly avec h=0.05 et ksmooth avec les mêmes deux paramètres.

# Prévisions locpoly
pred_loc <- approx(x = reg_locpoly$x, y = reg_locpoly$y, xout = T1)$y
# Prévisions ksmooth
pred_k05 <- estim1$y
pred_k005 <- estim2$y

# Calculs de l'EQM par rapport à fT1
EQM_loc <- mean((pred_loc - fT1)^2)
EQM_k05 <- mean((pred_k05 - fT1)^2)
EQM_k005 <- mean((pred_k005 - fT1)^2)

EQM <- data.frame(Méthode = c("locpoly h=0.05", "ksmooth h=0.2", "ksmooth h=0.05"),
                  EQM = c(EQM_loc, EQM_k05, EQM_k005))
EQM

Observation : en général, l’estimateur locpoly et l’estimateur ksmooth avec un paramètre de lissage raisonnable donnent des EQM du même ordre. Une bande trop large augmente l’EQM car la courbe lissée s’éloigne de la structure réelle, alors qu’une bande trop petite peut sur‑ajuster le bruit et également augmenter l’EQM. Le choix optimal de h se situe entre ces extrêmes.

Exercice 3 – Estimation non paramétrique à partir d’un petit jeu de données

On cherche à estimer la fonction de régression \(f\) par deux estimateurs de Nadaraya–Watson : l’un utilisant un noyau uniforme (fenêtre rectangulaire) et l’autre un noyau gaussien, tous deux avec une fenêtre de lissage h = 1.5. On comparera les courbes obtenues et on analysera les résidus.

On dispose d’un petit échantillon de cinq observations :

x3 <- c(-2, -1, -1, 1, 0)
y3 <- c( 2,  1,  4, 2, -2)
data3 <- data.frame(x = x3, y = y3)
data3

Estimateur de Nadaraya–Watson avec noyau uniforme

Pour calculer l’estimateur avec noyau uniforme, on utilise la fonction ksmooth avec l’option kernel = "box". On calcule les prédictions sur une grille fine de valeurs de x allant de min(x3) à max(x3).

grid_x <- seq(min(x3) - 0.5, max(x3) + 0.5, length.out = 200)
nw_uniform <- ksmooth(x3, y3, kernel = "box", bandwidth = 1.5, x.points = grid_x)

plot(x3, y3, pch = 16, cex = 1.2, col = "steelblue",
     xlab = "x", ylab = "y",
     main = "Exercice 3 – Nadaraya–Watson noyau uniforme (h = 1,5)")
grid()

lines(nw_uniform$x, nw_uniform$y, col = "darkorange", lwd = 2)

Estimateur de Nadaraya–Watson avec noyau gaussien

On répète l’opération avec un noyau gaussien (kernel = "normal") et la même bande de lissage h = 1.5.

nw_gaussian <- ksmooth(x3, y3, kernel = "normal", bandwidth = 1.5, x.points = grid_x)

plot(x3, y3, pch = 16, cex = 1.2, col = "steelblue",
     xlab = "x", ylab = "y",
     main = "Exercice 3 – Nadaraya–Watson noyau gaussien (h = 1,5)")
grid()

lines(nw_gaussian$x, nw_gaussian$y, col = "forestgreen", lwd = 2)

Comparaison et analyse des résidus

On superpose les deux estimateurs sur le même graphique et on calcule les résidus (observations moins valeurs ajustées) en chacun des points x3. Pour le noyau uniforme, on obtient les valeurs ajustées en interpolant la courbe lissée sur les points x3.

# Superposition des deux courbes
plot(x3, y3, pch = 16, cex = 1.2, col = "steelblue",
     xlab = "x", ylab = "y",
     main = "Exercice 3 – Comparaison des estimateurs")
grid()

lines(nw_uniform$x, nw_uniform$y, col = "darkorange", lwd = 2)
lines(nw_gaussian$x, nw_gaussian$y, col = "forestgreen", lwd = 2)

legend("topright", legend = c("Observations", "NW uniforme", "NW gaussien"),
       col = c("steelblue", "darkorange", "forestgreen"),
       lty = c(NA, 1, 1), pch = c(16, NA, NA), lwd = c(NA, 2, 2))

# Interpolation pour calculer les valeurs ajustées aux points x3
fit_uniform <- approx(x = nw_uniform$x, y = nw_uniform$y, xout = x3)$y
fit_gaussian <- approx(x = nw_gaussian$x, y = nw_gaussian$y, xout = x3)$y

resid_uniform <- y3 - fit_uniform
resid_gaussian <- y3 - fit_gaussian

data.frame(x = x3, y = y3, NW_uniform = fit_uniform, NW_gaussian = fit_gaussian,
           resid_uniform = resid_uniform, resid_gaussian = resid_gaussian)

Observation : les deux estimateurs donnent des courbes relativement similaires, mais le noyau gaussien produit une courbe plus lisse et moins anguleuse que le noyau uniforme. Les résidus montrent que certains points (notamment x = -1 avec y = 4) sont mal ajustés, car leur valeur est plus éloignée des autres observations. Dans un échantillon aussi petit, le choix du noyau et de la bande de lissage influe fortement sur l’aspect de la courbe.

Exercice 4 – Analyse de la vitesse d’un camion

On dispose des mesures de la vitesse instantanée V (en km/h) d’un camion en fonction du temps t (en secondes). Les données sont présentées dans un tableau et reproduites ci‑dessous.

t4 <- c(0, 5, 10, 15, 20, 30, 40, 50, 60, 70, 80, 90, 95, 100, 105, 110, 115, 120)
V4 <- c(0, 10, 20, 30, 40, 52, 60, 63, 65, 63, 58, 50, 45, 38, 29, 20, 12, 3)
camion <- data.frame(t = t4, V = V4)
camion

Représentation graphique et ajustement linéaire

On trace d’abord la vitesse en fonction du temps, puis on ajuste une droite de tendance au moyen d’une régression linéaire lm(V ~ t). On observe si une relation linéaire est plausible.

# Ajustement linéaire
reg_camion_lin <- lm(V ~ t, data = camion)

plot(camion$t, camion$V,
     pch = 16, cex = 1.0, col = "steelblue",
     xlab = "Temps (s)", ylab = "Vitesse (km/h)",
     main = "Exercice 4 – Vitesse du camion et ajustement linéaire")
grid()

abline(reg_camion_lin, col = "darkorange", lwd = 2)

summary(reg_camion_lin)

Call:
lm(formula = V ~ t, data = camion)

Residuals:
    Min      1Q  Median      3Q     Max 
-35.797 -16.838   2.468  19.876  28.468 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept) 35.79735    9.72052   3.683  0.00201 **
t            0.01224    0.13139   0.093  0.92693   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 22.55 on 16 degrees of freedom
Multiple R-squared:  0.0005421, Adjusted R-squared:  -0.06192 
F-statistic: 0.008679 on 1 and 16 DF,  p-value: 0.9269

Observation : la vitesse augmente d’abord, atteint un maximum autour de 65 km/h vers 60 secondes, puis décroît progressivement. La droite de tendance linéaire ne reproduit pas cette forme en cloche : elle sur‑estime la vitesse au début et la sous‑estime après 60 s. La relation n’est donc pas vraiment linéaire.

Ajustement quadratique : modèle \(v(t) = θ_0 + θ_1 t + θ_2 t^2\)

Puisque la relation semble parabolique, on ajuste un modèle de régression quadratique. On calcule les coefficients \(\theta_0\), \(\theta_1\) et \(\theta_2\) via lm(V ~ t + I(t^2)). On trace ensuite la courbe quadratique obtenue.

# Ajustement quadratique
reg_camion_quad <- lm(V ~ t + I(t^2), data = camion)
summary(reg_camion_quad)

Call:
lm(formula = V ~ t + I(t^2), data = camion)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.4624 -1.0905 -0.3722  0.7969  2.8613 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.1053526  0.8143340   1.357    0.195    
t            2.1613244  0.0346062  62.455   <2e-16 ***
I(t^2)      -0.0179830  0.0002813 -63.925   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.408 on 15 degrees of freedom
Multiple R-squared:  0.9963,    Adjusted R-squared:  0.9959 
F-statistic:  2044 on 2 and 15 DF,  p-value: < 2.2e-16

# Courbe quadratique
t_seq <- seq(min(camion$t), max(camion$t), length.out = 200)
pred_quad <- predict(reg_camion_quad, newdata = data.frame(t = t_seq))

plot(camion$t, camion$V,
     pch = 16, cex = 1.0, col = "steelblue",
     xlab = "Temps (s)", ylab = "Vitesse (km/h)",
     main = "Exercice 4 – Ajustement quadratique de la vitesse")
grid()

lines(t_seq, pred_quad, col = "forestgreen", lwd = 2)

Observation : la courbe quadratique épouse bien la trajectoire en cloche de la vitesse : accélération jusqu’à environ 60 s, puis décélération. Les coefficients obtenus reflètent cette forme : un terme quadratique négatif rend la fonction concave. Le modèle quadratique améliore nettement l’ajustement par rapport à la droite.

Ajustement quadratique via poly

On peut retrouver le même résultat en utilisant la fonction poly(t, deg = 2) qui génère des polynômes orthogonaux de degré 2. L’ajustement se fait ensuite avec lm(V ~ polynome).

reg_camion_poly <- lm(V ~ poly(t, 2), data = camion)
summary(reg_camion_poly)

Call:
lm(formula = V ~ poly(t, 2), data = camion)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.4624 -1.0905 -0.3722  0.7969  2.8613 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   36.556      0.332 110.115   <2e-16 ***
poly(t, 2)1    2.101      1.409   1.492    0.157    
poly(t, 2)2  -90.035      1.409 -63.925   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.408 on 15 degrees of freedom
Multiple R-squared:  0.9963,    Adjusted R-squared:  0.9959 
F-statistic:  2044 on 2 and 15 DF,  p-value: < 2.2e-16

# Comparaison des prédictions
t_seq <- seq(min(camion$t), max(camion$t), length.out = 200)
pred_poly <- predict(reg_camion_poly, newdata = data.frame(t = t_seq))

plot(camion$t, camion$V,
     pch = 16, cex = 1.0, col = "steelblue",
     xlab = "Temps (s)", ylab = "Vitesse (km/h)",
     main = "Exercice 4 – Ajustement polynomial (degré 2)")
grid()

lines(t_seq, pred_poly, col = "darkred", lwd = 2)

Observation : l’utilisation de poly(t, deg=2) donne un ajustement identique à celui obtenu avec I(t^2). La syntaxe est plus compacte et permet de généraliser facilement à des polynômes de degré supérieur. Les courbes se superposent parfaitement.

Exercice 5 – Régression non paramétrique sur des données de verre

On s’intéresse à des mesures portant sur des fragments de verre. On veut expliquer l’indice de réfraction RI en fonction de la teneur en aluminium Al. Les données proviennent du jeu de données Glass du package mlbench.

Importation et représentation des données

On charge le jeu de données, on conserve les colonnes d’intérêt et on réalise un nuage de points. On appelle verres le data.frame contenant les variables Al et RI.

library(mlbench)
data(Glass)

# Création du data.frame verres
verres <- data.frame(Al = Glass$Al, RI = Glass$RI)
head(verres)

plot(verres$Al, verres$RI,
     pch = 16, cex = 0.7, col = "steelblue",
     xlab = "Aluminium (%)", ylab = "Indice de réfraction (RI)",
     main = "Exercice 5 – Nuage de points : RI en fonction de Al")
grid()

Observation : la relation entre Al et RI semble non linéaire : l’indice de réfraction diminue pour des teneurs en aluminium faibles, puis augmente légèrement. On explore maintenant des méthodes non paramétriques pour estimer cette relation.

Estimateur par noyau avec la fonction npreg

Le package np fournit des outils d’estimation non paramétrique multivariée. On commence par déterminer la largeur de bande optimale à l’aide de npregbw, puis on ajuste le modèle avec npreg. On peut visualiser la courbe ajustée et les intervalles de confiance asymptotiques.

library(np)
attach(verres)                # attache des variables Al et RI
h_np <- npregbw(RI ~ Al)      # estimation automatique de la bande

Multistart 1 of 1 |
Multistart 1 of 1 |
Multistart 1 of 1 |
Multistart 1 of 1 /
Multistart 1 of 1 |
Multistart 1 of 1 |
                   

summary(h_np)

Regression Data (214 observations, 1 variable(s)):

Regression Type: Local-Constant
Bandwidth Selection Method: Least Squares Cross-Validation
Formula: RI ~ Al
Bandwidth Type: Fixed
Objective Function Value: 6.854372e-06 (achieved on multistart 1)

Exp. Var. Name: Al Bandwidth: 0.132813 Scale Factor: 1.247614

Continuous Kernel Type: Second-Order Gaussian
No. Continuous Explanatory Vars.: 1
Estimation Time: 0.022 seconds

verres_np <- npreg(h_np)      # ajustement non paramétrique
plot(verres_np, plot.errors.method = "asymptotic")

detach(verres)

Observation : l’estimation indique une courbe lisse qui suit la forme générale du nuage de points. Les bandes d’erreur asymptotiques montrent où l’incertitude est plus grande, notamment aux extrémités où les données sont moins nombreuses.

Estimateur à noyau via ksmooth

Une méthode plus simple consiste à utiliser ksmooth avec un noyau gaussien et une bande de lissage fixée, ici bandwidth = 5. On trace cette courbe pour la comparer à celle obtenue avec npreg.

smooth_ks <- ksmooth(verres$Al, verres$RI, kernel = "normal", bandwidth = 5, x.points = verres$Al)

plot(verres$Al, verres$RI,
     pch = 16, cex = 0.7, col = "steelblue",
     xlab = "Aluminium (%)", ylab = "RI",
     main = "Exercice 5 – ksmooth (bandwidth = 5)")
grid()

lines(smooth_ks$x, smooth_ks$y, col = "darkorange", lwd = 2)

Observation : la courbe ksmooth est relativement lisse et suit bien la tendance moyenne. Cependant, le choix de la bande bandwidth = 5 peut être arbitraire ; une bande trop large sur‑lisse les données, tandis qu’une bande trop petite conduit à une courbe plus ondulée.

Comparaison avec un ajustement par polynômes locaux (loess)

Enfin, on ajuste la relation RI ~ Al par des polynômes locaux via loess avec un paramètre span = 0.5 et une famille gaussienne. On trace la courbe obtenue et on la compare à celle de ksmooth.

verres_loess <- loess(RI ~ Al, data = verres, span = 0.5, family = "gaussian")
pred_loess <- predict(verres_loess, newdata = data.frame(Al = verres$Al))

plot(verres$Al, verres$RI,
     pch = 16, cex = 0.7, col = "steelblue",
     xlab = "Aluminium (%)", ylab = "RI",
     main = "Exercice 5 – Comparaison des ajustements")
grid()

lines(smooth_ks$x, smooth_ks$y, col = "darkorange", lwd = 2)
lines(verres$Al, pred_loess, col = "forestgreen", lwd = 2)

legend("topright", legend = c("ksmooth", "loess"),
       col = c("darkorange", "forestgreen"), lty = 1, lwd = 2)

Observation : la courbe loess (vert) suit étroitement la forme du nuage et peut capter des variations locales grâce à un paramètre span modéré. La courbe ksmooth (orange) est légèrement plus lisse. Selon la valeur de span ou de bandwidth, l’une ou l’autre méthode peut offrir un meilleur compromis entre biais et variance.

FIN