Add TP2

M2/Statistiques Non Paramétrique/TP2.Rmd

@@ -0,0 +1,197 @@
+# Estimateur de Nadaraya-Watson
+
+1. 
+```{r}
+n = 500
+sigma2 = 0.01
+epsilon = rnorm(n, 0, sqrt(sigma2))
+
+X = runif(n, 0, 1)
+Y = exp(abs(X - 0.5)) + epsilon
+
+plot(X, Y, pch = 3)
+```
+
+2. 
+
+```{r}
+library(np)
+
+bw = npregbw(Y ~ X)
+ghat = npreg(bw)
+x = seq(0, 1, length.out = 100)
+
+plot(ghat)
+points(x, exp(abs(x - 0.5)), type = 'l', col = 'darkgreen')
+legend(
+  'top',
+  c(expression(hat(g)[h]), 'g'),
+  col = c(1, 'darkgreen'),
+  lty = rep(1, 2)
+)
+```
+
+3.  
+
+...
+
+# Estimation par projection
+
+4. 
+
+```{r}
+mmax = 100
+Dmax = 2 * mmax + 1
+phiX = matrix(1, n, Dmax)
+for (j in 1:mmax) {
+  phiX[, 2 * j] = sqrt(2) * cos(2 * pi * j * X)
+  phiX[, 2 * j + 1] = sqrt(2) * sin(2 * pi * j * X)
+}
+G = crossprod(phiX) / n
+```
+
+5. 
+
+```{r}
+hatb = crossprod(phiX, Y) / n
+hata = solve(G, hatb)
+```
+
+6. 
+
+```{r}
+p = 500
+t = seq(0, 1, length.out = p)
+phi = matrix(1, Dmax, p)
+
+for (j in 1:mmax) {
+  phi[2 * j, ] = sqrt(2) * cos(2 * pi * j * t)
+  phi[2 * j + 1, ] = sqrt(2) * sin(2 * pi * j * t)
+}
+
+hatg = crossprod(phi, hata)
+
+plot(t, hatg, type = 'l', col = 'red')
+points(t, exp(abs(t - 0.5)), type = 'l', col = 'darkgreen')
+legend(
+  'topright',
+  c('g', expression(hat(g)[501])),
+  lty = c(1, 1),
+  col = c('darkgreen', 'red')
+)
+```
+
+7. 
+
+
+```{r}
+hata1 = solve(G[1, 1], hatb[1])
+hatg1 = hata1 * phi[1, ]
+
+hata3 = solve(G[1:3, 1:3], hatb[1:3])
+hatg3 = crossprod(phi[1:3, ], hata3)
+
+hata5 = solve(G[1:5, 1:5], hatb[1:5])
+hatg5 = crossprod(phi[1:5, ], hata5)
+
+plot(t, hatg, type = 'l', col = 'red', main = 'Projection estimators')
+points(t, exp(abs(t - 0.5)), type = 'l', col = 'darkgreen')
+points(t, hatg1, type = 'l', col = 'blue')
+points(t, hatg3, type = 'l', col = 'green')
+points(t, hatg5, type = 'l', col = 'purple')
+legend(
+  'topright',
+  c(
+    'g',
+    expression(hat(g)[1]),
+    expression(hat(g)[3]),
+    expression(hat(g)[5]),
+    expression(hat(g)[501])
+  ),
+  lty = rep(1, 5),
+  col = c('darkgreen', 'blue', 'green', 'purple', 'red')
+)
+```
+
+```{r}
+kappa = 2.5
+mmax = 10
+gamman = rep(NA, mmax)
+penm = sigma2 * (2 * (1:mmax) + 1) / n
+for (m in 1:mmax) {
+  hatam = solve(G[1:(2 * m + 1), 1:(2 * m + 1)], hatb[1:(2 * m + 1)])
+  hatgmX = crossprod(t(phiX[, 1:(2 * m + 1)]), hatam)
+  gamman[m] = mean((Y - hatgmX)^2)
+}
+crit = gamman + kappa * penm
+plot(
+  (2 * (1:mmax) + 1),
+  gamman,
+  type = 'l',
+  ylim = range(c(gamman, crit, kappa * penm)),
+  xlab = 'm',
+  ylab = 'crit(m)',
+  col = 'blue'
+)
+points((2 * (1:mmax) + 1), kappa * penm, type = 'l', col = 'red')
+points((2 * (1:mmax) + 1), crit, type = 'l')
+legend(
+  'left',
+  c('contraste', 'penalité', 'critère'),
+  col = c('blue', 'red', 1),
+  lty = rep(1, 3)
+)
+hatm = which.min(crit)
+points(2 * hatm + 1, crit[hatm], pch = 16)
+```
+
+```{r}
+hatahatm = solve(G[1:(2 * hatm + 1), 1:(2 * hatm + 1)], hatb[1:(2 * hatm + 1)])
+hatghatm = crossprod(phi[1:(2 * hatm + 1), ], hatahatm)
+plot(
+  t,
+  hatghatm,
+  type = 'l',
+  col = 'green',
+  main = 'Selected projection estimators',
+  ylab = expression(hat(g)[hat(m)](t))
+)
+points(t, exp(abs(t - 0.5)), type = 'l', col = 'darkgreen')
+```
+
+# Application aux données réelles
+
+9. 
+
+
+```{r}
+data(cps71)
+attach(cps71)
+plot(age, logwage, xlab = "Age", ylab = "log(wage)")
+```
+
+
+```{r}
+bw = npregbw(logwage ~ age)
+ghat = npreg(bw)
+
+plot(ghat, xlab = "Age", ylab = "log(wage)")
+points(age, logwage)
+```
+
+10. 
+
+```{r}
+data(Italy)
+attach(Italy)
+
+plot(year, gdp, xlab = "Année", ylab = "GDP")
+```
+
+
+```{r}
+bw = npregbw(gdp ~ year)
+ghat = npreg(bw)
+
+plot(ghat, xlab="Année", ylab="gdp",type='l')
+```