---
title: "Przewidywanie cen mieszkań"
output: 
  html_notebook: 
    theme: cerulean
    highlight: textmate
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(warning = FALSE, message = FALSE)
```


Dwa poprzednie przykłady były problemami klasyfikacji, w których celem było przewidzenie etykiety opisującej dane wejściowe. Innym spotykanym często problemem uczenia maszynowego jest regresja, która polega na przewidywaniu wartości o charakterze ciągłym, a nie dyskretnej etykiety. Przykładem regresji jest przewidywanie jutrzejszej temperatury na podstawie zgromadzonych danych meteorologicznych lub przewidywanie czasu potrzebnego na skończenie oprogramowania na podstawie jego specyfikacji..

Nie myl regresji z algorytmem regresji logistycznej. Regresja logistyczna wbrew swej nazwie nie jest algorytmem regresji. Jest algorytmem klasyfikacji.

## Zbiór cen mieszkań w Bostonie


Będziemy starali się przewidzieć medianę cen mieszkań w podmiejskich dzielnicach Bostonu w połowie lat 70. XX w. na podstawie danych określających konkretną dzielnicę, takich jak współczynnik przestępczości czy lokalny podatek od nieruchomości. Tym razem będziemy korzystali ze zbioru danych, który różni się od wcześniejszych zbiorów dwiema rzeczami. Ma dość mało elementów: tylko 506 (podzielono je na zbiór treningowy zawierający 404 próbki i zbiór testowy zawierający 102 próbki). Każda cecha danych wejściowych (przykładem cechy jest współczynnik przestępczości) jest wyrażona w innej skali. Niektóre wartości są ułamkami przyjmującymi wartości od 0 do 1, inne przyjmują wartości z zakresu od 1 do 12, a jeszcze inne — od 0 do 100 itd.

Przyjrzyjmy się danym:

```{r}
library(keras)

dataset <- dataset_boston_housing()
c(c(train_data, train_targets), c(test_data, test_targets)) %<-% dataset
```

```{r}
str(train_data)
```

```{r}
str(test_data)
```

Jak widać, zbiór treningowy składa się z 404 próbek, a zbiór testowy ze 102 próbek. Każda próbka jest opisana za pomocą 13 cech numerycznych:

 1.    Współczynnik przestępczości
 2.    Część działek o powierzchni przekraczającej 2300 metrów kwadratowych.
 3.    Część powierzchni działek działalności gospodarczych nie zajmujących się sprzedażą.
 4.    Sztuczna zmienna rzeki Charles River (= 1 jeżeli działka znajduje się przy rzece; 0 w pozostałych przypadkach).
 5.    Stopień koncentracji tlenków azotu(cząstek na 10 milionów).
 6.    Średnia liczba pokoi w budynku.
 7.    Część mieszkań zajętych przez właścicieli, które zostały wybudowane przed rokiem 1940.
 8.    Średnia ważona odległości od pięciu stref, w których znajduje się najwięcej zakładów pracy.
 9.    Indeks dostępności dróg szybkiego ruchu
10.    Pełny podatek od nieruchomości za 10000 $.
11.    Stosunek liczby uczniów do liczby nauczycieli.
12.    1000 * (Bk - 0.63) ** 2, gdzie Bk ułamkiem czarnoskórej ludności miasta.
13.    % niższego statusu populacji.


Celem jest określenie median wartości domów zamieszkanych przez właścicieli — wyrażonych w tysiącach dolarów:

```{r}
str(train_targets)
```

Ceny zwykle wahają się od 10 000 dolarów do 50 000 dolarów. Tanio? Pamiętaj o tym, że to wartości z połowy lat 70. Ceny w tym zbiorze nie uwzględniają inflacji.

## Przygotowywanie danych


Ładowanie do sieci neuronowej wartości należących do kilku różnych zakresów może sprawić problem. Niektóre sieci są w stanie automatycznie dopasować do siebie tak różne dane, ale z pewnością utrudni to proces uczenia. Najlepszą praktyką podczas pracy z takimi danymi jest przeprowadzenie normalizacji poszczególnych cech: w przypadku każdej cechy danych wejściowych (kolumny macierzy danych wejściowych) należy przeprowadzić operację odejmowania od wartości średniej i dzielenia przez odchylenie standardowe — wówczas wartości cech zostaną wyśrodkowane wokół zera i będą charakteryzowały się jednostkowym odchyleniem standardowym. Operację tę można z łatwością przeprowadzić dzięki funkcji scale().

```{r}
mean <- apply(train_data, 2, mean)
std <- apply(train_data, 2, sd)
train_data <- scale(train_data, center = mean, scale = std)
test_data <- scale(test_data, center = mean, scale = std)
```

Zwróć uwagę na to, że wielkości używane podczas normalizacji testowego zbioru danych są obliczane na podstawie treningowego zbioru danych. Nigdy nie powinniśmy korzystać z wartości obliczonych na podstawie treningowego zbioru danych. Dotyczy to nawet tak prostych zadań jak normalizacja danych.

## Budowanie sieci

Dysponujemy małą liczbą próbek, a więc zbudujemy bardzo małą sieć zawierającą dwie warstwy ukryte, składające się z 64 jednostek każda. Ogólnie rzecz biorąc, im mniejszą ilością danych treningowych dysponujemy, tym bardziej jesteśmy narażeni na nadmierne dopasowanie sieci. W celu zminimalizowania efektu nadmiernego dopasowania można między innymi korzystać z małej sieci.

```{r}
# Będziemy tworzyć wiele instancji tego samego modelu, 
# a więc konstruując je, będziemy korzystać z funkcji.
build_model <- function() {
  model <- keras_model_sequential() %>% 
    layer_dense(units = 64, activation = "relu", 
                input_shape = dim(train_data)[[2]]) %>% 
    layer_dense(units = 64, activation = "relu") %>% 
    layer_dense(units = 1) 
    
  model %>% compile(
    optimizer = "rmsprop", 
    loss = "mse", 
    metrics = c("mae")
  )
}
```

Sieć kończy się pojedynczą jednostką bez funkcji aktywacji (jest to warstwa liniowa). To typowe rozwiązanie stosowane w regresji skalarnej (regresji, w której próbuje się przewidzieć jedną wartość o charakterze ciągłym). Zastosowanie funkcji aktywacji ograniczyłoby zakres wartości wyjściowych możliwych do wygenerowania. Gdybyśmy zastosowali w ostatniej warstwie tej sieci funkcję aktywacji sigmoid, to sieć mogłaby generować tylko wartości znajdujące się w zakresie od 0 do 1. W praktyce zastosowaliśmy ostatnią warstwę o charakterze liniowym, a więc możemy przewidywać dowolne wartości.

Zwróć uwagę na to, że sieć jest kompilowana z funkcją straty mse (średniego błędu kwadratowego). Funkcja ta oblicza kwadrat różnicy między wartościami przewidywanymi przez sieć i wartościami docelowymi. Ta funkcja straty jest często używana w czasie rozwiązywania problemów regresji.

Podczas trenowania monitorowana jest nowa metryka: średni błąd bezwzględny (mae). Jest to bezwzględna wartość różnicy między wartościami przewidywanymi przez sieć a wartościami docelowymi. Średnia wartość błędu bezwzględnego o wartości równej np. 0,5 w przypadku tego problemu oznacza, że przewidywane ceny średnio odbiegają od wartości docelowych o 500 dolarów.

## K-składowa walidacja krzyżowa

W celu oceny sprawności działania sieci podczas dostrajania jej parametrów, takich jak liczba epok trenowania, możemy tak jak wcześniej podzielić dane treningowe na podzbiór treningowy i podzbiór walidacyjny, ale nasz zbiór jest na tyle mały, że podzbiór walidacyjny utworzony w ten sposób byłby bardzo mały (zawierałby np. tylko 100 elementów). W związku z tym wynik walidacji mógłby ulegać dużej zmianie w zależności od tego, które elementy treningowego zbioru danych byłyby używane podczas walidacji, a które podczas trenowania. Wyniki walidacji mogłyby charakteryzować się dużą wariancją zależną od podziału zbioru testowego na podzbiór testowy i walidacyjny. W takiej sytuacji nie można przeprowadzić wiarygodnej walidacji.

W takim przypadku najlepiej jest skorzystać z walidacji krzyżowej k-składowych (patrz rysunek 3.11). Polega ona na podziale dostępnych danych na k części (zwykle 4 lub 5), utworzeniu k identycznych modeli i trenowaniu każdego z nich na k – 1 częściach zbioru oraz przeprowadzaniu ewaluacji na pozostałej, nieużytej wcześniej części zbioru dostępnych danych. Wynik walidacji modelu jest średnią wyników walidacji wszystkich składowych modeli. Kod implementujący to rozwiązanie jest dość prosty:


```{r, echo=TRUE, results='hide'}
k <- 4
indices <- sample(1:nrow(train_data))
folds <- cut(1:length(indices), breaks = k, labels = FALSE) 

num_epochs <- 100
all_scores <- c()
for (i in 1:k) {
  cat("processing fold #", i, "\n")
  # Przygotuj dane walidacyjne: dane z k-tej składowej.
  val_indices <- which(folds == i, arr.ind = TRUE) 
  val_data <- train_data[val_indices,]
  val_targets <- train_targets[val_indices]
  
  # Przygotuj dane treningowe: dane z pozostałych składowych.
  partial_train_data <- train_data[-val_indices,]
  partial_train_targets <- train_targets[-val_indices]
  
  # Zbuduj model Keras (model został skompilowany wcześniej).
  model <- build_model()
  
  # Trenuj model w trybie cichym (parametr verbose = 0).
  model %>% fit(partial_train_data, partial_train_targets,
                epochs = num_epochs, batch_size = 1, verbose = 0)
                
  # Przeprowadź ewaluację modelu przy użyciu danych walidacyjnych.
  results <- model %>% evaluate(val_data, val_targets, verbose = 0)
  all_scores <- c(all_scores, results$mean_absolute_error)
}  
```

```{r}
all_scores
```

```{r}
mean(all_scores)
```

Podczas poszczególnych iteracji uzyskujemy dość zróżnicowane wartości walidacji (od 2,1 do 2,8). Średnia wartość (2,4) jest wartością, na której można o wiele bardziej polegać niż na poszczególnych wynikach walidacji składowych — właśnie to chcieliśmy uzyskać, stosując k-składową walidację krzyżową. W tym przypadku odchodzimy od wartości docelowych średnio o 2400 dolarów, co jest znaczącą kwotą przy cenach znajdujących się w zakresie od 10 000 do 50 000 dolarów.
Spróbujmy wydłużyć proces trenowania do 500 epok. W celu obserwacji wydajności modelu w każdej epoce zmodyfikujemy pętlę treningową tak, aby zapisywała wynik walidacji poszczególnych epok w dzienniku pracy:

```{r}
# Some memory clean-up
k_clear_session()
```

```{r, echo=TRUE, results='hide'}
num_epochs <- 500
all_mae_histories <- NULL
for (i in 1:k) {
  cat("processing fold #", i, "\n")
  
  # Przygotowuje dane walidacyjne: dane z k-tej składowej.
  val_indices <- which(folds == i, arr.ind = TRUE)
  val_data <- train_data[val_indices,]
  val_targets <- train_targets[val_indices]
  
  # Przygotowuje dane treningowe: dane z pozostałych składowych.
  partial_train_data <- train_data[-val_indices,]
  partial_train_targets <- train_targets[-val_indices]
  
  # Buduje model Keras (model został skompilowany wcześniej).
  model <- build_model()
  
  # Przeprowadza ewaluację modelu przy użyciu danych walidacyjnych.
  history <- model %>% fit(
    partial_train_data, partial_train_targets,
    validation_data = list(val_data, val_targets),
    epochs = num_epochs, batch_size = 1, verbose = 0
  )
  mae_history <- history$metrics$val_mean_absolute_error
  all_mae_histories <- rbind(all_mae_histories, mae_history)
}
```

Teraz możemy obliczyć średni wynik walidacji wszystkich składowych poszczególnych epok.

```{r}
average_mae_history <- data.frame(
  epoch = seq(1:ncol(all_mae_histories)),
  validation_mae = apply(all_mae_histories, 2, mean)
)
```

Przedstawmy te dane na wykresie:

```{r}
library(ggplot2)
ggplot(average_mae_history, aes(x = epoch, y = validation_mae)) + geom_line()
```

Wykres może okazać się dość trudny do przeanalizowania z powodu skali i dużej wariancji. Spróbujmy uzyskać bardziej czytelny obraz za pomocą funkcji geom_smooth():

```{r}
ggplot(average_mae_history, aes(x = epoch, y = validation_mae)) + geom_smooth()
```

Z ostatniego wykresu wynika, że średni błąd bezwzględny przestaje ulegać poprawie po 125 epokach. Po przekroczeniu tego punktu model zaczyna ulegać przeuczeniu.

Po zakończeniu dostrajania pozostałych parametrów modelu (poza liczbą epok możemy zmienić również rozmiar warstw ukrytych) przeprowadzamy trenowanie ostatecznej wersji modelu na całym zbiorze danych treningowych (w procesie tym korzystamy z optymalnych parametrów), a następnie sprawdzamy jego wydajność na zbiorze testowym.

```{r, echo=FALSE, results='hide'}
# Utwórz nową, skompilowaną wersję modelu.
model <- build_model()

# Trenuj model na całym zbiorze danych treningowych.
model %>% fit(train_data, train_targets,
          epochs = 80, batch_size = 16, verbose = 0)

result <- model %>% evaluate(test_data, test_targets)
```

```{r}
result
```

We are still off by about \$2,680.

## Wnioski


Oto wnioski, które należy wynieść z tego przykładu:

* Regresję przeprowadza się przy użyciu innej funkcji straty od tej, z której korzystaliśmy podczas klasyfikacji. W czasie rozwiązywania problemów dotyczących regresji często stosuje się funkcję straty mse (średniego błędu kwadratowego).
* Podczas pracy nad problemem regresji stosuje się również inne metryki ewaluacyjne (nie używa się tych samych metryk, które są stosowane w problemach klasyfikacji). Oczywiście wynika to z tego, że koncepcja dokładności nie odnosi się do regresji. Metryką stosowaną w regresji jest średni błąd bezwzględny realizowany przez funkcję mae.
* Jeżeli cechy danych wejściowych przyjmują wartości o różnych zakresach, to każda cecha powinna zostać niezależnie przeskalowana podczas wstępnej obróbki danych.
* W pracy z małym zbiorem danych ewaluację modelu warto jest oceniać przy użyciu techniki k-składowej walidacji krzyżowej.
* W pracy z małym zbiorem danych lepiej jest korzystać z małej sieci z kilkoma warstwami ukrytymi (zwykle stosuje się od 1 do 2 warstw), co pozwoli uniknąć przeuczenia modelu.
