# coding: utf-8

# # Statystyka praktyczna w data science (Python)
# # Rozdział 1. Badania eksploracyjne
# > (c) 2019 Peter C. Bruce, Andrew Bruce, Peter Gedeck

# Importuje wymagane pakiety Pythona.

from pathlib import Path
import pandas as pd
import numpy as np
from scipy.stats import trim_mean
from statsmodels import robust
import wquantiles
import seaborn as sns
import matplotlib.pylab as plt
from matplotlib.collections import EllipseCollection
from matplotlib.colors import Normalize


# Definiuje ścieżki do zestawów danych. Jeżeli przechowujesz dane w innym katalogu, wprowadź poniżej stosowne zmiany.



DATA = Path('.').resolve().parents[1] / 'dane'

AIRLINE_STATS_CSV = DATA / 'airline_stats.csv'
KC_TAX_CSV = DATA / 'kc_tax.csv.gz'
LC_LOANS_CSV = DATA / 'lc_loans.csv'
AIRPORT_DELAYS_CSV = DATA / 'dfw_airline.csv'
SP500_DATA_CSV = DATA / 'sp500_data.csv.gz'
SP500_SECTORS_CSV = DATA / 'sp500_sectors.csv'
STATE_CSV = DATA / 'state.csv'


# # Miary położenia
# ## Przykład: miara położenia dla wielkości populacji i wskaźnika morderstw



# Tabela 1.2
state = pd.read_csv(STATE_CSV)
print(state.head(8))


# Oblicza średnią, średnią ucinaną i medianę kolumny Population. W przypadku `średniej` i `mediany` możemy skorzystać z metod ramki danych w pakiecie _pandas_. Średnia ucinana wymaga funkcji `trim_mean` z pakietu _scipy.stats_.



state = pd.read_csv(STATE_CSV)
print(state['Population'].mean())




print(trim_mean(state['Population'], 0.1))




print(state['Population'].median())


# Średnia ważona jest dostępna w pakiecie numpy. Medianę ważoną obliczymy za pomocą wyspecjalizowanego pakietu `wquantiles` (https://pypi.org/project/wquantiles/).



print(state['Murder.Rate'].mean())




print(np.average(state['Murder.Rate'], weights=state['Population']))




print(wquantiles.median(state['Murder.Rate'], weights=state['Population']))


# # Miary rozproszenia



# Tabela 1.2
print(state.head(8))


# Odchylenie standardowe:



print(state['Population'].std())


# Przedział międzykwartylowy obliczamy jako różnicę pomiędzy 75. a 25. percentylem.



print(state['Population'].quantile(0.75) - state['Population'].quantile(0.25))


# Medianę odchylenia bezwględnego z mediany można obliczyć za pomocą metody dostępne w pakiecie _statsmodels_



print(robust.scale.mad(state['Population']))
print(abs(state['Population'] - state['Population'].median()).median() / 0.6744897501960817)


# ## Percentyle i boxploty
# Pakiet _Pandas_ zawiera metodę `quantile` dostosowaną do ramek danych.



print(state['Murder.Rate'].quantile([0.05, 0.25, 0.5, 0.75, 0.95]))




# Tabela 1.4
percentages = [0.05, 0.25, 0.5, 0.75, 0.95]
df = pd.DataFrame(state['Murder.Rate'].quantile(percentages))
df.index = [f'{p * 100}%' for p in percentages]
print(df.transpose())


# Pakiet _Pandas_ zawiera wiele podstawowych wykresów eksploracyjnych; do tej kategorii należą boxploty.



ax = (state['Population']/1_000_000).plot.box(figsize=(3, 4))
ax.set_ylabel('Populacja (w milionach)')

plt.tight_layout()
plt.show()


# ## Tablica częstości i histogramy
# Metoda `cut` z pakietu _pandas_ rozdziela dane z zestawu danych na przedziały. Metoda ta zawiera wiele argumentów. W poniższym listingu tworzymy przedziały o takich samych rozmiarach. Metoda `value_counts` zwraca tablicę częstości.



binnedPopulation = pd.cut(state['Population'], 10)
print(binnedPopulation.value_counts())




# Tabela 1.5
binnedPopulation.name = 'binnedPopulation'
df = pd.concat([state, binnedPopulation], axis=1)
df = df.sort_values(by='Population')

groups = []
for group, subset in df.groupby(by='binnedPopulation'):
    groups.append({
        'Zakres': group,
        'Liczba': len(subset),
        'Stany': ','.join(subset.Abbreviation)
    })
print(pd.DataFrame(groups))


# Pakiet _Pandas_ obsługuje także histogramy w analizie eksploracyjnej.



ax = (state['Population'] / 1_000_000).plot.hist(figsize=(4, 4))
ax.set_xlabel('Populacja (w milionach)')

plt.tight_layout()
plt.show()


# ## Szacowanie gęstości
# Gęstość stanowi rozwiązanie alternatywne dla histogramów; dzięki niej widać wyraźniej rozkład punktów danych. Gładkość krzywej gęstości określamy za pomocą argumentu `bw_method`.



ax = state['Murder.Rate'].plot.hist(density=True, xlim=[0, 12], 
                                    bins=range(1,12), figsize=(4, 4))
state['Murder.Rate'].plot.density(ax=ax)
ax.set_xlabel('Wskaźnik morderstw (na 100 000 osób)')

plt.tight_layout()
plt.show()


# # Badanie danych binarnych i skategoryzowanych



# Tabela 1.6
dfw = pd.read_csv(AIRPORT_DELAYS_CSV)
print(100 * dfw / dfw.values.sum())


# Pakiet _Pandas_ obsługuje rónież wykresy słupkowe służące do wyświetlania poszczególnych zmiennych skategoryzowanych.



ax = dfw.transpose().plot.bar(figsize=(4, 4), legend=False)
ax.set_xlabel('Powód opóźnienia')
ax.set_ylabel('Liczba')

plt.tight_layout()
plt.show()


# # Korelacja
# Najpierw wczytujemy wymagane zestawy danych:



sp500_sym = pd.read_csv(SP500_SECTORS_CSV)
sp500_px = pd.read_csv(SP500_DATA_CSV, index_col=0)




# Tabela 1.7
# Określa symbole telekomunikacyjne
telecomSymbols = sp500_sym[sp500_sym['sector'] == 'telecommunications_services']['symbol']

# Filtruje dane dla zakresu dat od lipca 2012 r. do czerwca 2015 r.
telecom = sp500_px.loc[sp500_px.index >= '2012-07-01', telecomSymbols]
telecom.corr()
print(telecom)


# Następnie skoncentrujmy się na zwrotach dla głównych funduszy (sector == 'etf'). 



etfs = sp500_px.loc[sp500_px.index > '2012-07-01', 
                    sp500_sym[sp500_sym['sector'] == 'etf']['symbol']]
print(etfs.head())


# Z powodu znacznej liczby kolumn w tej tabeli analiza macierzy korelacji okazuje się bardzo żmudna i o wiele wygodniej jest zilustrować korelację w postaci mapy cieplnej. Pakiet _seaborn_ zawiera wygodną implementację map cieplnych.



fig, ax = plt.subplots(figsize=(5, 4))
ax = sns.heatmap(etfs.corr(), vmin=-1, vmax=1, 
                 cmap=sns.diverging_palette(20, 220, as_cmap=True),
                 ax=ax)

plt.tight_layout()
plt.show()


# Powyższa mapa cieplna jest czytelna w przypadku, gdy mamy możliwość oglądania jej w kolorze. W przypadku wykresów czarno-białych musimy także zwizualizować kierunek. Poniższy listing wizualizuje siłę korelację za pomocą elips.




def plot_corr_ellipses(data, figsize=None, **kwargs):
    ''' https://stackoverflow.com/a/34558488 '''
    M = np.array(data)
    if not M.ndim == 2:
        raise ValueError('Dane muszą być tablicą dwuwymiarową.')
    fig, ax = plt.subplots(1, 1, figsize=figsize, subplot_kw={'aspect':'equal'})
    ax.set_xlim(-0.5, M.shape[1] - 0.5)
    ax.set_ylim(-0.5, M.shape[0] - 0.5)
    ax.invert_yaxis()

    # współrzędne xy środka każdej elipsy
    xy = np.indices(M.shape)[::-1].reshape(2, -1).T

    # wyznacza względne rozmiary osi dłuższej/krótszej zgodnie z siłą korelacji dodatniej/ujemnej.
    w = np.ones_like(M).ravel() + 0.01
    h = 1 - np.abs(M).ravel() - 0.01
    a = 45 * np.sign(M).ravel()

    ec = EllipseCollection(widths=w, heights=h, angles=a, units='x', offsets=xy,
                           norm=Normalize(vmin=-1, vmax=1),
                           transOffset=ax.transData, array=M.ravel(), **kwargs)
    ax.add_collection(ec)

    # jeżeli dane mieszczą się w obiekcie DataFrame, nazwy rzędu/kolumny zostają użyte jako etykiety jednostek na osiach
    if isinstance(data, pd.DataFrame):
        ax.set_xticks(np.arange(M.shape[1]))
        ax.set_xticklabels(data.columns, rotation=90)
        ax.set_yticks(np.arange(M.shape[0]))
        ax.set_yticklabels(data.index)

    return ec

m = plot_corr_ellipses(etfs.corr(), figsize=(5, 4), cmap='bwr_r')
cb = fig.colorbar(m)
cb.set_label('Współczynnik korelacji')

plt.tight_layout()
plt.show()


# ## Wykres punktowy
# Pakiet _pandas_ obsługuje proste wykresy punktowe. Wartość znacznika `$\u25EF$` wyznacza okrąg dla każdego punktu.



ax = telecom.plot.scatter(x='T', y='VZ', figsize=(4, 4), marker='$\u25EF$')
ax.set_xlabel('ATT (T)')
ax.set_ylabel('Verizon (VZ)')
ax.axhline(0, color='grey', lw=1)
ax.axvline(0, color='grey', lw=1)

plt.tight_layout()
plt.show()




ax = telecom.plot.scatter(x='T', y='VZ', figsize=(4, 4), marker='$\u25EF$', alpha=0.5)
ax.set_xlabel('ATT (T)')
ax.set_ylabel('Verizon (VZ)')
ax.axhline(0, color='grey', lw=1)
print(ax.axvline(0, color='grey', lw=1))


# # Badanie dwóch lub więcej zmiennych
# Wczytujemy zestaw danych kc_tax i filtrujemy go na podstawie zróżnicowania kryteriów:



kc_tax = pd.read_csv(KC_TAX_CSV)
kc_tax0 = kc_tax.loc[(kc_tax.TaxAssessedValue < 750000) & 
                     (kc_tax.SqFtTotLiving > 100) &
                     (kc_tax.SqFtTotLiving < 3500), :]
print(kc_tax0.shape)


# ## Wykres hexagonal binning  i kontury 
# ### Przedstawianie danych numerycznych względem danych numerycznych

# W przypadku dużej liczby punktów danych wykres punktowy przestaje być czytelny. Przydatniejsze stają się metody wizualizujące gęstość. Jednym z lepszych rozwiązań jest metoda `hexbin` z ramek danych _pandas_.



ax = kc_tax0.plot.hexbin(x='SqFtTotLiving', y='TaxAssessedValue',
                         gridsize=30, sharex=False, figsize=(5, 4))
ax.set_xlabel('Powierzchnia domu')
ax.set_ylabel('Wartość podatku')

plt.tight_layout()
plt.show()


# Metoda kdeplot z pakietu _seaborn_ stanowi dwuwymiarowe rozszerzenie wykresu gęstości. 



fig, ax = plt.subplots(figsize=(4, 4))
ax = sns.kdeplot(kc_tax0.SqFtTotLiving, kc_tax0.TaxAssessedValue, ax=ax)
ax.set_xlabel('Powierzchnia domu')
ax.set_ylabel('Wartość podatku')

plt.tight_layout()
plt.show()


# ## Dwie zmienne skategoryzowane
# Wczytujemy zestaw danych `lc_loans`:



lc_loans = pd.read_csv(LC_LOANS_CSV)




# Tabela 1.8(1)
crosstab = lc_loans.pivot_table(index='grade', columns='status', 
                                aggfunc=lambda x: len(x), margins=True)
print(crosstab)




# Tabela 1.8(2)
df = crosstab.copy().loc['A':'G',:]
df.loc[:,'Charged Off':'Late'] = df.loc[:,'Charged Off':'Late'].div(df['All'], axis=0)
df['All'] = df['All'] / sum(df['All'])
perc_crosstab = df
print(perc_crosstab)


# ## Dane kategoryzowane i numeryczne
# Boxploty kolumny w pakiecie _Pandas_ można pogrupować ze względu na różne kolumny.



airline_stats = pd.read_csv(AIRLINE_STATS_CSV)
airline_stats.head()
ax = airline_stats.boxplot(by='airline', column='pct_carrier_delay',
                           figsize=(5, 5))
ax.set_xlabel('')
ax.set_ylabel('Procent opóźnionych lotów')
plt.suptitle('')

plt.tight_layout()
plt.show()


# Pakiet _Pandas_ obsługuje także odmianę boxplotów o zwaną _wykresami skrzypcowymi_. 



fig, ax = plt.subplots(figsize=(5, 5))
sns.violinplot(airline_stats.airline, airline_stats.pct_carrier_delay,
               ax=ax, inner='quartile', color='white')
ax.set_xlabel('')
ax.set_ylabel('Procent opóźnionych lotów')

plt.tight_layout()
plt.show()


# ## Wizualizacja wielu zmiennych



zip_codes = [98188, 98105, 98108, 98126]
kc_tax_zip = kc_tax0.loc[kc_tax0.ZipCode.isin(zip_codes),:]
kc_tax_zip

def hexbin(x, y, color, **kwargs):
    cmap = sns.light_palette(color, as_cmap=True)
    plt.hexbin(x, y, gridsize=25, cmap=cmap, **kwargs)

g = sns.FacetGrid(kc_tax_zip, col='ZipCode', col_wrap=2)
g.map(hexbin, 'SqFtTotLiving', 'TaxAssessedValue', 
      extent=[0, 3500, 0, 700000])
g.set_axis_labels('Powierzchnia domu', 'Wartość podatku')
g.set_titles('Kod pocztowy: {col_name:.0f}')

plt.tight_layout()
plt.show()










