mirror of
https://github.com/ArthurDanjou/ArtStudies.git
synced 2026-02-02 21:31:32 +01:00
199 lines
3.4 KiB
Plaintext
199 lines
3.4 KiB
Plaintext
# 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)
|
|
|
|
print(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')
|
|
``` |