KDE
MŁ
Estymacja gęstości
1D
Losujemy dane - dwa rozkłady normalne o średnich 1 i 5 i odchyleniach odpowiednio 1 i 2. Powinny być dwie górki w 1 i w 5, a ta w 5 powinna być bardziej spłaszczona od tej w 1.
set.seed(123)
x <- c(rnorm(100, 0, 1), rnorm(100, 5, 2))Na domyślnych argumentach:
den_x <- density(x)
bw_x <- den_x$bw
png('kde1.png')
par(mfrow=c(1,3))
#1)
hist(x)
#2)
plot(den_x)
#3)
hist(x, freq = FALSE, ylim=c(0, 0.17))
lines(den_x$x, den_x$y)
text(9.5, 0.15, paste0("bw=", round(bw_x, 4)))
dev.off()
Nie na domyślnych: density(x, kernel=..., bw=...). Tak możemy zmienić jądro i parametr wygładzający. bw można podać jako liczbę lub jako nazwę (np. domyślnie jest bw = "nrd0").
den_x_1 <- density(x, bw=0.1)
den_x_2 <- density(x, bw=3)
png('kde2.png')
par(mfrow=c(1,2))
#1)
hist(x, freq = FALSE, ylim=c(0, 0.25))
lines(den_x_1$x, den_x_1$y)
#2)
hist(x, freq = FALSE, ylim=c(0, 0.25))
lines(den_x_2$x, den_x_2$y)
dev.off()
2D
set.seed(321)
x <- c(rnorm(1000, 0, 1), rnorm(1000, 5, 2))
y <- c(rnorm(1000, 0, 1), rnorm(1000, 5, 2))
png('kde3.png')
par(mfrow=c(1,2))
f <- MASS::kde2d(x, y)
persp(f, col="orange", phi=20, theta=50, xlab="x", ylab="y", zlab="Denstiy")
h1 <- bw.SJ(x)
h2 <- bw.SJ(y)
f.SJ <- MASS::kde2d(x, y, n=50, h = c(h1, h2))
persp(f.SJ, col="orange", phi=20, theta=50, xlab="x", ylab="y", zlab="Denstiy")
dev.off()
Ćwiczenie / Praca domowa dla zainteresowanych
Dotyczy ostatnej strony wykładu (tej pisanej ręcznie, nie slajdów).
Będziemy generować obserwacje z rozkładu \(\hat{f}\), a rozkład \(\hat{f}\) otrzymamy na podstawie \(20\) obserwacji. Czyli po kolei:
- wygeneruj \(20\) obserwacji z rozkładu jednostajnego na przedziale 0-10
- dopasuj estymator jądrowy: jądro powinno być takie, żeby dało się z niego losować, parametr wygładzający dowolny, ale znany (najlepiej nie za duży, ale tylko po to, żeby dało się sprawdzić, czy procedura działa - łatwiej porównywać górki i dołki)
- wykonaj algorytm z końca wykładu tzn. 1) wylosuj indeks \(i\in\{1, 2, \ldots, 20\}\), 2) wylosuj \(\epsilon \sim K\), 3) obserwacja z rozkładu \(\hat{f}\) to \(Y = X_{i} + \epsilon h\). Powtórz tyle razy, ile nowych obserwacji chcesz otrzymać (powiedzmy 1000)
- do nowych obserwacji dopasuj estymator jądrowy \(\hat{g}\)
- porównaj \(\hat{f}\) i \(\hat{g}\)
Przykładowe wyniki poniżej, ale można eksperymentować z różnymi parametrami i licznościami prób.
Podpunkt pierwszy:
set.seed(123)
X <- runif(20, 0, 10)
png('kde4.png')
plot(X, rep(1, 20), type="h", ylab="")
dev.off()
Podpunkt 2:
Podpunkt 3 i 4 (Y to nowe obserwacje):
Podpunkt 3 i 4 dla N = 10 000:
