# File-Name:       chapter02.R           
# Date:            2011-11-11                                
# Author:          John Myles White
# Email:           jmw@johnmyleswhite.com
# Purpose:         Kod do rozdziału 2. Prezentuje narzędzia do eksploracyjnej analizy danych.
# Data Used:       data/01_heights_weights_genders.csv
# Packages Used:   ggplot2
# Machine:         John Myles White's MacBook

# All source code is copyright (c) 2011, under the Simplified BSD License.  
# For more information on FreeBSD see: http://www.opensource.org/licenses/bsd-license.php

# All images and materials produced by this code are licensed under the Creative Commons 
# Attribution-Share Alike 3.0 United States License: http://creativecommons.org/licenses/by-sa/3.0/us/

# All rights reserved.

#
# Fragment 1
#

# Wczytywanie zbioru danych z dysku.
#options(device="png")
#png(filename = "Rplot%03d.png",
#    width = 8, height = 8, units = "in", pointsize = 12,
#     bg = "white",  res = 300)
#png(filename = "Rplot%03d.png", width = 800, height = 800, units = "px")


data.file <- file.path('data', '01_heights_weights_genders.csv')
heights.weights <- read.csv(data.file, header = TRUE, sep = ',')

# Przeliczamy jednostki i zamieniamy symbole płci
heights.weights$Height <- heights.weights$Height * 2.54
heights.weights$Weight <- heights.weights$Weight * 0.454
heights.weights$Gender <- gsub("Male", "Mężczyzna", heights.weights$Gender)
heights.weights$Gender <- gsub("Female", "Kobieta", heights.weights$Gender)

# Tworzymy wektor liczbowy zawierający same dane o wzroście.
heights <- with(heights.weights, Height)
summary(heights)

# Oczekiwany wynik:
# Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#137.8   161.3   168.4   168.6   175.7   200.7 

#
# Fragment 2
#

# Definiujemy własne funkcje średniej i mediany.
my.mean <- function(x)
{
  return(sum(x) / length(x))
}

my.median <- function(x)
{
  sorted.x <- sort(x)
  if (length(x) %% 2 == 0)
  {
    indices <- c(length(x) / 2, length(x) / 2 + 1)
    return(mean(sorted.x[indices]))
  }
  else
  {
    index <- ceiling(length(x) / 2)
    return(sorted.x[index])
  }
}

#
# Fragment 3
#

# Porównujemy średnie i mediany na prostych przykładach.
my.vector <- c(0, 100)

my.vector
# [1]	0 100

mean(my.vector)
#[1] 50

median(my.vector)
#[1] 50

my.vector <- c(0, 0, 100)

mean(my.vector)
#[1] 33.33333

median(my.vector)
#[1] 0

#
# Fragment 4
#

# Sprawdzimy, czy nasze własne funkcje mediany i średniej zwracają poprawne wyniki.
my.mean(heights)
#[1] 66.36756

my.median(heights)
#[1] 66.31807

mean(heights) - my.mean(heights)
#[1] 0

median(heights) - my.median(heights)
#[1] 0

#
# Fragment 5
#

# Eksperymenty z funkcjami pozwalającymi na ocenę zakresu danych w zbiorze.
min(heights)
#[1] 54.26313

#
# Fragment 6
#

max(heights)
#[1] 78.99874

#
# Fragment 7
#

c(min(heights), max(heights))
#[1] 54.26313 78.99874

range(heights)
#[1] 54.26313 78.99874

#
# Fragment 8
#

# Wypróbujmy funkcję 'quantile' do wyznaczenia dowolnie zdefiniowanych kwantyli.
quantile(heights)
#      0%      25%      50%      75%     100% 
#54.26313 63.50562 66.31807 69.17426 78.99874 

#
# Fragment 9
#

quantile(heights, probs = seq(0, 1, by = 0.20))
#      0%      20%      40%      60%      80%     100% 
#54.26313 62.85901 65.19422 67.43537 69.81162 78.99874 

#
# Fragment 10
#

seq(0, 1, by = 0.20)
#[1] 0.0 0.2 0.4 0.6 0.8 1.0

#
# Fragment 11
#

# Definiujemy własną funkcję wariancji, do szacowania rozrzutu danych.
my.var <- function(x)
{
  m <- mean(x)
  return(sum((x - m) ^ 2) / length(x))
}

#
# Fragment 12
#

# Testowanie własnej funkcji wariancji pod kątem poprawności obliczeń.
my.var(heights) - var(heights)

#
# Fragment 13
#

# Poprawienie funkcji wariancji w celu wyeliminowania skrzywienia wyników.
my.var <- function(x)
{
  m <- mean(x)
  return(sum((x - m) ^ 2) / (length(x) - 1))
}

# Ponowne testy funkcji wariancji pod kątem poprawności.
my.var(heights) - var(heights)

#
# Fragment 14
#

# Sprawdzenie zakresu przewidzianego przez funkcję wariancji
c(mean(heights) - var(heights), mean(heights) + var(heights))
#[1] 51.56409 81.17103

#
# Fragment 15
#

c(mean(heights) - var(heights), mean(heights) + var(heights))
#[1] 51.56409 81.17103
range(heights)
#[1] 54.26313 78.99874

#
# Fragment 16
#

# Własna funkcja odchylenie standardowego do oceny zakresów danych.
my.sd <- function(x)
{
  return(sqrt(my.var(x)))
}

#
# Fragment 17
#

# Sprawdzenie funkcji pod kątem poprawności.
my.sd(heights) - sd(heights)

#
# Fragment 18
#

c(mean(heights) - sd(heights), mean(heights) + sd(heights))
# [1] 62.52003 70.21509

range(heights)
#[1] 54.26313 78.99874

#
# Fragment 19
#

c(mean(heights) - sd(heights), mean(heights) + sd(heights))
# [1] 62.52003 70.21509

c(quantile(heights, probs = 0.25), quantile(heights, probs = 0.75))
#     25%      75% 
#63.50562 69.17426 

#
# Fragment 20
#

# Zacznijmy wizualizować dane za pomocą pakietu ggplot2.
library('ggplot2')

# Wczytanie danych od nowa, dla pewności.
data.file <- file.path('data', '01_heights_weights_genders.csv')
heights.weights <- read.csv(data.file, header = TRUE, sep = ',')

# Przeliczenie jednostek i zamiana symboli płci
heights.weights$Height <- heights.weights$Height * 2.54
heights.weights$Weight <- heights.weights$Weight * 0.454
heights.weights$Gender <- gsub("Male", "Mężczyzna", heights.weights$Gender)
heights.weights$Gender <- gsub("Female", "Kobieta", heights.weights$Gender)


# Eksperymenty z rozdzielczością histogramów
ggplot(heights.weights, aes(x = Height)) +
  geom_histogram(binwidth = 1) + labs(x = "Wzrost", y = "liczba")

#
# Fragment 21
#

ggplot(heights.weights, aes(x = Height)) +
  geom_histogram(binwidth = 5) + labs(x = "Wzrost", y = "liczba")

#
# Fragment 22
#

ggplot(heights.weights, aes(x = Height)) +
  geom_histogram(binwidth = 0.01) + labs(x = "Wzrost", y = "liczba")

#
# Fragment 23
#

# Eksperymenty z estymata gęstości jądrowej
ggplot(heights.weights, aes(x = Height)) +
  geom_density() + labs(x = "Wzrost", y = "gęstość")

#
# Fragment 24
#

# Rozdzielenie obrazu rozkładu wzrostu i wagi na podstawie płci
ggplot(heights.weights, aes(x = Height, fill = Gender)) +
  geom_density() + labs(x = "Wzrost", y = "gęstość", fill = "Płeć")

#
# Fragment 25
#

ggplot(heights.weights, aes(x = Weight, fill = Gender)) +
  geom_density() + labs(x = "Waga", y = "gęstość", fill = "Płeć")

#
# Fragment 26
#

# Pojedynczy wykres dwuaspektowy, lepiej obrazujący ukrytą strukturę danych.
ggplot(heights.weights, aes(x = Weight, fill = Gender)) +
  geom_density() + labs(x = "Wzrost", y = "gęstość", fill = "Płeć") +
  facet_grid(Gender ~ .)

#
# Fragment 27
#

# Eksperymenty z liczbami losowymi z rozkładu normalnego.
m <- 0
s <- 1
ggplot(data.frame(X = rnorm(100000, m, s)), aes(x = X)) +
  geom_density() + ylab("gęstość")

#
# Fragment 28
#

# Porównanie rozkładu normalnego z rozkładem Cauchy'ego.
set.seed(1)
normal.values <- rnorm(250, 0, 1)
cauchy.values <- rcauchy(250, 0, 1)
range(normal.values)
range(cauchy.values)

#
# Fragment 29
#

ggplot(data.frame(X = normal.values), aes(x = X)) +
  geom_density() + ylab("gęstość")
ggplot(data.frame(X = cauchy.values), aes(x = X)) +
  geom_density() + ylab("gęstość")

#
# Fragment 30
#

# Eksperymenty z liczbami losowymi z rozkładu gamma.
gamma.values <- rgamma(100000, 1, 0.001)
ggplot(data.frame(X = gamma.values), aes(x = X)) +
  geom_density() + ylab("gęstość")

#
# Fragment 31
#

# Generowanie wykresów relacji wzrostu do wagi, w celu uwidocznienia zależności.
ggplot(heights.weights, aes(x = Height, y = Weight)) +
  geom_point() + labs(x = "Wzrost", y = "Waga")

#
# Fragment 32
#

# Dodanie wygładzonej linii relacji.
ggplot(heights.weights, aes(x = Height, y = Weight)) +
  geom_point() + labs(x = "Wzrost", y = "Waga") +
  geom_smooth()

#
# Fragment 33
#

# Zauważmy, że wygładzona linia relacji poprawia się wraz ze wzrostem ilości danych.
ggplot(heights.weights[1:20, ], aes(x = Height, y = Weight)) +
  geom_point() + labs(x = "Wzrost", y = "Waga") +
  geom_smooth()
ggplot(heights.weights[1:200, ], aes(x = Height, y = Weight)) +
  geom_point() + labs(x = "Wzrost", y = "Waga") +
  geom_smooth()
ggplot(heights.weights[1:2000, ], aes(x = Height, y = Weight)) +
  geom_point() + labs(x = "Wzrost", y = "Waga") +
  geom_smooth()

#
# Fragment 34
#

# Wizualizacja zależności płci od wzrostu i wagi.
ggplot(heights.weights, aes(x = Height, y = Weight)) +
  geom_point(aes(color = Gender, alpha = 0.25)) +
  labs(x = "Wzrost", y = "Waga", color = "Płeć") +
  scale_alpha(guide = "none") + 
  scale_color_manual(values = c("Mężczyzna" = "black", "Kobieta" = "gray")) +
  theme_bw()

# Wersja barwna.
ggplot(heights.weights, aes(x = Height, y = Weight, color = Gender)) +
  geom_point() + labs(x = "Wzrost", y = "Waga", color = "Płeć")

#
# Fragment 35
#

heights.weights <- transform(heights.weights,
                             Male = ifelse(Gender == 'Mężczyzna', 1, 0))

logit.model <- glm(Male ~ Weight + Height,
                   data = heights.weights,
                   family = binomial(link = 'logit'))

ggplot(heights.weights, aes(x = Height, y = Weight)) +
  geom_point(aes(color = Gender, alpha = 0.25)) +
  labs(x = "Wzrost", y = "Waga", color = "Płeć") +
  scale_alpha(guide = "none") + 
  scale_color_manual(values = c("Mężczyzna" = "black", "Kobieta" = "gray")) +
  theme_bw() +
  stat_abline(intercept = -coef(logit.model)[1] / coef(logit.model)[2],
              slope = - coef(logit.model)[3] / coef(logit.model)[2],
              geom = 'abline',
              color = 'black')
