• Jak oszacować \(f\)?
• Możemy oszacować \(f(x) = \beta_0 + \beta_1 x + \beta_2 x^2 + \beta_3 x^3\)
• Lub \(f(x) = \beta_0 + \beta_1x^{\beta_2}\)
• OK jeśli znasz odpowiednią postać.
• Ale lepiej często założyć, że \(f\) jest ciągła i wygładzona.

• \(K\) jest funkcją typu kernel, gdzie: \(\int K = 1\), \(K(a)=K(-a)\), oraz \(K(0)\ge K(a)\), dla wszystkich \(a\).
• \(\hat{f}\) jest ważoną średnią ruchomą
• Musimy dobrać \(K\) oraz \(h\).

\[K(x)= \begin{cases} \frac{1}{2} & -1 < x < 1 \\ 0 & \text{dla pozostałych}. \end{cases}\]

\[K(x) = \begin{cases} \frac{3}{4}(1-x^2) & -1 < x < 1 \\ 0 & \text{dla pozostałych}. \end{cases}\]

\[K(x) = \begin{cases} c(1-|x|^3)^3 & -1 < x < 1 \\ 0 & \text{dla pozostałych}. \end{cases}\]

\[K(x) = \frac{1}{\sqrt{2\pi}} e^{-\frac{1}{2}x^2}\]

• Gładkie jądro jest lepsze, ale poza tym wybór jądra nie robi wielkiej różnicy.
• Optymalnym jądrem (minimalizującym MSE) jest Epanechnikov. Jest ono również szybkie.
• Wybór wartości \(h\) jest kluczowy.

\[CV(h) = \frac{1}{n} \sum_{j=1}^n (y_j - \hat{f}^{(-j)}_h(x_j))^2\]

• \((-j)\) wskazuje \(j\)-ty punkt pominięty w oszacowaniu
• Wybierz \(h\), które minimalizuje CV.
• Działa ok pod warunkiem, że nie ma zduplikowanych par \((x,y)\). Ale sporadycznie pojawiają się dziwne wyniki.

Niech \(y=f(x) + \varepsilon\) gdzie \(\varepsilon \sim IID(0,\sigma^2)\). Wtedy

• \(r_K = \displaystyle\int K^2(x) dx\), \(v_K = \displaystyle\int x^2K(x)dx\)

• Estymator spójny, jeśli \(nh\rightarrow\infty\) oraz \(h\rightarrow0\) jako \(n\rightarrow\infty\).

• \(h_{\text{optimal}} \approx\displaystyle \left(\frac{r_K}{n[f''(x)]^2v_K^2}\right)^{1/5}\)

Dostępnych jest wiele pakietów. Jeden z lepszych to KernSmooth:

fit <- locpoly(faithful$eruptions, faithful$waiting,
        degree=0, bandwidth=0.3) %>% as.tibble
ggplot(faithful) +
  geom_point(aes(x=eruptions,y=waiting)) +
  geom_line(data=fit, aes(x=x,y=y), col='blue')

Wygładzanie kernelem jest metodą o stałej lokalnej:

\(\text{WLS}(x) = \sum_{j=1}^n w_j(x) (y_j - a_0)^2\)

gdzie

\(w_j(x) = \frac{K\left(\frac{x-x_j}{h}\right)}{\sum_{i=1}^n K\left(\frac{x-x_i}{h}\right)}\)

jest minimalizowane przez

\(\hat{f}(x) = \hat{a}_0 = \sum_{j=1}^n w_j(x)y_j\)

Zamiast tego możemy obliczyć lokalny liniowy:

\[\hat{f}(x) = \hat{a}_0\]

\[\hat{f}(x) = \hat{a}_0\]

• Lokalne liniowe i lokalne kwadratowe są powszechnie stosowane.
• Odporna regresja może być użyta jako zamiennik!
• Mniej stronniczy na granicach niż wygładzanie jądra!
• Lokalna kwadratura mniej nieobiektywna przy szczytach i dołkach niż lokalna liniowa lub kernelowa

• Bardzo wygodną funkcją jestKernSmooth::locpoly(x, y, degree, bandwidth)
• Wykorzystuje się tutaj kernel Gaussowski.
• dpill może być użyte do wyboru bandwidth \(h\) jako degree=1.
• Wiele osób wydaje się używać metody prób i błędów — znalezienie największego \(h\), który uchwyci to, co myślą, że widzą na oko.

Najbardziej znaną implementacją jest loess (lokalnie kwadratowy)

fit <- loess(y ~ x, span=0.75, degree=2,
  family="gaussian", data)

• Używa się tutaj jądra tri-cube i zmiennej szerokości pasma.
• span kontroluje szerokość pasma. Określone w procentach pokrytych danych.
• degree jest rzędem wielomianu
• Użyj family="symmetric" dla solidnego dopasowania

smr <- loess(waiting ~ eruptions, data=faithful)
ggplot(faithful) +
    geom_point(aes(x=eruptions,y=waiting)) +
    ggtitle("Old Faithful (Loess, span=0.75)") +
    geom_line(aes(x=eruptions, y=fitted(smr)),
      col='blue')

ggplot(exa) +
    geom_point(aes(x=x,y=y)) +
    geom_smooth(aes(x=x,y=y), method='loess',
      span=0.22)

• Ponieważ lokalne wielomiany wykorzystują lokalne modele liniowe, możemy łatwo znaleźć błędy standardowe dla dopasowanych wartości.

• Połączone razem, tworzą one punktowy przedział ufności.

• Automatycznie tworzone przy użyciu geom_smooth.

• Parametry ograniczone tak, aby \(f(x)\) było ciągłe.
• Dalsze ograniczenia nałożone, aby dać ciągłe pochodne.
• Kwadratowe sploty są najczęściej występujące, dla ciągłych \(f'\), \(f''\)

• \(\lambda\) jest parametrem wygładzania, który należy wybrać
• \(\int [f''(x)]^2 dx\) to miara chropowatości.
• Rozwiązanie: \(\hat{f}\) jest splajnem sześciennym z węzłami \(\kappa_i=x_i\), \(i=1,\dots,n\) (ignorując duplikaty).
• Inne ograniczenia prowadzą do splotów wyższego rzędu
• Walidacja krzyżowa może być użyta do wyboru \(\lambda\).

smr <- smooth.spline(
  faithful$eruptions,
  faithful$waiting,
  cv=TRUE)
smr <- data.frame(x=smr$x,y=smr$y)
ggplot(faithful) +
  geom_point(aes(x=eruptions,y=waiting)) +
  ggtitle("Old Faithful (Sploty interpolujące, lambda wybrana przez CV)") +
  geom_line(data=smr, aes(x=x, y=y), col='blue')

• Mniej węzłów niż splajnów wygładzających
• Potrzeba wyboru węzłów, a nie parametru wygładzania.
• Można oszacować jako model liniowy po wybraniu węzłów.

• Niech \(\kappa_1<\kappa_2<\cdots<\kappa_K\) będzie ``knots’’ w przedziale \((a,b)\).
• Niech \(x_1=x\), \(x_2 = x^2\), \(x_3=x^3\), \(x_j = (x-\kappa_{j-3})_+^3\) dla \(j=4,\dots,K+3\).
• Wtedy regresja \(y\) wobec \(x_1,\dots,x_{K+3}\) jest sześcienna, ale gładka w węzłach.
• Wybór węzłów może być trudny i arbitralny.
• Algorytmy automatycznego wyboru węzłów są bardzo wolne.
• Często używa się równo rozłożonych węzłów. Wtedy wystarczy wybrać tylko \(K\).

• B-splajny dostarczają równoważnego zestawu funkcji bazowych.
• Naturalne splajny sześcienne są wariacją na temat B-splajnów z liniowymi warunkami brzegowymi.
• Są one zazwyczaj bardziej stabilne
• Zaimplementowane w funkcji splines::ns w R
• Można określić węzły jawnie, lub df. Wtedy df-2 węzły są wybierane na kwantylach \(x\).

fit <- lm(waiting ~ ns(eruptions, df=6), faithful)
ggplot(faithful) +
    geom_point(aes(x=eruptions,y=waiting)) +
    ggtitle("Old Faithful (Naturalne sploty, 6 df)") +
    geom_line(aes(x=eruptions, y=fitted(fit)), col='blue')

ggplot(exa) +
    geom_point(aes(x=x,y=y)) +
    geom_smooth(aes(x=x,y=y), method='gam',
      formula = y ~ s(x,k=12))

• Ponieważ splajny regresyjne używają lokalnych modeli liniowych, możemy łatwo znaleźć błędy standardowe dla dopasowanych wartości.

• Połączone razem, tworzą one punktowy przedział ufności.

• Automatycznie tworzone przy użyciu geom_smooth.

Większość metod w naturalny sposób rozszerza się na więcej wymiarów.

• Metody jądra wielowymiarowego
• Wielowymiarowe lokalne powierzchnie kwadratowe
• Cienkie splajny (2-d wersja splajnów wygładzających)

Większość danych leży w pobliżu granicy.

x <- matrix(runif(1e6,-1,1), ncol=100)
boundary <- function(z) { any(abs(z) > 0.95) }

mean(apply(x[,1,drop=FALSE], 1, boundary))

## [1] 0.0532

mean(apply(x[,1:2], 1, boundary))

## [1] 0.0985

mean(apply(x[,1:5], 1, boundary))

## [1] 0.2285

Większość danych leży w pobliżu granicy.

x <- matrix(runif(1e6,-1,1), ncol=100)
boundary <- function(z) { any(abs(z) > 0.95) }

mean(apply(x[,1:10], 1, boundary))

## [1] 0.4031

mean(apply(x[,1:50], 1, boundary))

## [1] 0.9221

mean(apply(x[,1:100], 1, boundary))

## [1] 0.9933

Dane są niekompletne!

x <- matrix(runif(1e6,-1,1), ncol=100)
nearby <- function(z) { all(abs(z) < 0.5) }
mean(apply(x[,1,drop=FALSE], 1, nearby))

## [1] 0.4966

mean(apply(x[,1:2], 1, nearby))

## [1] 0.2397

mean(apply(x[,1:5], 1, nearby))

## [1] 0.0286

Dane są niekompletne!

x <- matrix(runif(1e6,-1,1), ncol=100)
nearby <- function(z) { all(abs(z) < 0.5) }
mean(apply(x[,1:10], 1, nearby))

## [1] 7e-04

mean(apply(x[,1:50], 1, nearby))

## [1] 0

mean(apply(x[,1:100], 1, nearby))

## [1] 0

• Ilość danych dostępnych w oknie jest proporcjonalna do \(n^{-d}\).
• Niech \(h= 0.5\), i załóżmy, że potrzebujemy 100 obserwacji, aby oszacować nasz model lokalnie.

lomod <- loess(sr ~ pop15 + ddpi, data=savings)

library(mgcv)
smod <- gam(sr ~ s(pop15, ddpi), data=savings)

Zestaw danych “Prestige” (z pakietu “car”) zawiera dane nt. prestiżu n=102 Kanadyjskich zawodów z 1971 roku, a także średni dochód w danym zawodzie. Do zbadania zależności między prestiżem a dochodem wykorzystaj metody regresji nieparametrycznej.

##                     education income women prestige census type
## gov.administrators      13.11  12351 11.16     68.8   1113 prof
## general.managers        12.26  25879  4.02     69.1   1130 prof
## accountants             12.77   9271 15.70     63.4   1171 prof
## purchasing.officers     11.42   8865  9.11     56.8   1175 prof
## chemists                14.62   8403 11.68     73.5   2111 prof
## physicists              15.64  11030  5.13     77.6   2113 prof

Najpierw załaduj dane i zwizualizuj relację pomiędzy dochodem (X) a prestiżem (Y).

Uwaga: związek wygląda na nieliniowy. Dla zawodów, które zarabiają mniej niż $10K, istnieje silna (pozytywna) liniowa zależność pomiędzy dochodem a prestiżem. Jednak w przypadku zawodów, które zarabiają od 10 do 25 tysięcy dolarów, związek ten ma znacznie inne (osłabione) nachylenie.

Zbiór danych “mcycle” (z pakietu MASS) zawiera n=133 pary punktów czasowych (w ms) i obserwowanych przyspieszeń głowy (w g), które zostały zarejestrowane w symulowanym wypadku motocyklowym.

Do zbadania zależności między czasem a przyspieszeniem wykorzystaj metody regresji nieparametrycznej.

##   times accel
## 1   2.4   0.0
## 2   2.6  -1.3
## 3   3.2  -2.7
## 4   3.6   0.0
## 5   4.0  -2.7
## 6   6.2  -2.7

Najpierw wczytaj dane i zwizualizuj zależność między czasem (X) a przyspieszeniem (Y).

Uwaga: zależność wygląda na nieliniową.

Przyspieszenie jest stabilne od 0-15 ms, spada od ok. 15-20 ms, rośnie od 20-30 ms, spada od 30-40 ms, a następnie zaczyna się stabilizować.