#!/usr/bin/env python
# coding: utf-8

# In[1]:


import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

get_ipython().run_line_magic('matplotlib', 'inline')


# In[2]:


data = pd.read_csv('student/student-mat.csv', sep=';')
print(data.shape)
data.head()


# In[3]:


data.columns


#  - school - тип школы ("GP" - Gabriel Pereira или "MS" - Mousinho da Silveira)
#  - sex - пол ("F" - female или "M" - male)
#  - age - возраст (от 15 до 22)
#  - address - откуда студент ("U" - urban или "R" - rural)
#  - famsize - размер семьи ("LE3" - меньше или равно 3 или "GT3" - больше 3)
#  - Pstatus - в каких отношениях родители ("T" - живут вместе "A" - раздельно)
#  - Medu - образование матери (0 - никакого,  1 - начальное образование (4 класса), 2 – от 5 до 9 классов, 3 – среднеспециальное или 4 – высшее)
#  - Fedu - образование отца (0 - никакого,  1 - начальное образование (4 класса), 2 – от 5 до 9 классов, 3 – среднеспециальное или 4 – высшее)
#  - Mjob - работа матери ("teacher", "health" care related, civil "services" (e.g. administrative or police), "at_home" or "other")
#  - Fjob - работа отца ("teacher", "health" care related, civil "services" (e.g. administrative or police), "at_home" or "other")
#  - reason - причина выбора школы (близко к дому — "home", репутация школы — "reputation", "course" предпочтение некоторым предметам или "other")
#  - guardian - опекун ("mother", "father" или "other")
#  - traveltime - время от дома до школы (1 - меньше 15 мин., 2 - 15 до 30 мин., 3 - 30 мин. до 1 часа, или 4 - больше 1 часа)
#  - studytime - количество часов обучения в неделю (1 - меньше 2 часов, 2 - от 2 до 5 часов, 3 - от 5 до 10 часов, или 4 - больше 10 часов)
#  - failures - колисечтво ранее не сданных предметов (n if 1 <= n < 3, else 4)
#  - schoolsup - дополнительные занятия (yes or no)
#  - famsup - помощь от семьи при выполнении заданий (yes or no)
#  - paid - дополнительные платные занятия (yes or no)
#  - activities - внеклассная деятельность (yes or no)
#  - nursery - посещал детский сад (yes or no)
#  - higher - желание высшего образования (yes or no)
#  - internet - домашний интернет (yes or no)
#  - romantic - состоит в романтических отношениях (yes or no)
#  - famrel - насколько хорошо отношения в семье (от 1 - очень плохие до 5 - превосходные)
#  - freetime - наличие свободного времени после школы (от 1 - очень мало до 5 - очень много)
#  - goout - гуляет с друзьями (от 1 - редко до 5 - очень часто)
#  - Dalc - употребление алкоголя в будние дни (от 1 - очень редко до 5 - очень часто)
#  - Walc - употребление алкоголя в выходные (от 1 - очень редко до 5 - очень часто)
#  - health - текущее состояние здоровья (от 1 - очень плохое до 5 - очень хорошее)
#  - absences - количество школьных пропусков (от 0 до 93)
#  - G1 - оценка за первый семестр (от 0 до 20)
#  - G2 - оценка за второй семестр (от 0 до 20)
#  - G3 - итоговая оценка (от 0 до 20)

# Добавим в наши данные целевую переменную из задания и посмотрим на ее распределение:

# In[4]:


data['Talc'] = (data['Dalc'] * 5 + data['Walc'] * 2) / 7.0


# In[5]:


data.hist(['Dalc', 'Walc', 'Talc'], layout=(1, 3), figsize=(15, 5), bins=30);


# **Наблюдение:** распределение Dalc и Walc не очень похожи..
# 
# **Вывод:** может быть стоит настраивать две регрессии: одну для предсказания Dalc, другую для Walc и пробовать комбинировать результаты?

# В данном датасете присутствуют признаки разных типов: вещественные, категориальные, порядковые и бинарные. Не стоит путать порядковые и категориальные признаки! Для первых мы, например, можем, осмысленно считать корреляцию, а для вторых нет.

# In[6]:


binary_features = []
categorical_features = []
real_features = []
targets = ['Walc', 'Dalc', 'Talc']
for feature_name in data.drop(targets, axis=1):
    num_values = len(data[feature_name].unique())
    if num_values == 2:
        binary_features.append(feature_name)
        continue
    if num_values < 10:
        categorical_features.append(feature_name)
        continue
    real_features.append(feature_name)


# In[7]:


print('Binary: {}\nCategorical: {}\nReal: {}'.format(
        ", ".join(binary_features), 
        ", ".join(categorical_features), 
        ", ".join(real_features)))


# ### Посмотрим на вещественные переменные:

# In[8]:


sns.pairplot(data[['Talc'] + real_features])


# Видно несколько особенностей. Во-первых, графики зависимости переменных G{1,2,3} друг от друга выглядят линейно. Это означает, что эти признаки сильно коррелированы (они почти линейно зависимы с небольшими случайными отклонениями). Как мы знаем, это в среднем не очень хорошо, потому что мультиколлинеарность ведет к переобучению. Давайте посмотрим на коэффициенты корреляции между этими переменными.

# In[9]:


np.round(data.corr(), 2)


# Глядя на эту табличку стоит обратить внимание на несколько вещей. Во-первых, видно, что корреляция посчиталась не для всего
# множества переменных, а только для тех признаков, для которых осмысленно можно говорить о линейной зависимости (если у вас
# признаки -- это индикаторы некоторых неупорядоченных множеств, то нельзя сказать, что увеличение одного из них ведет к
# увеличению другого). Во-вторых, мы можем еще раз подтвердить наше наблюдение, что признаки G1, G2 и G3 сильно коррелированы.
# Обратите внимание, что сами по себе значения корреляции ничего не показывают! Они могут быть высокими из-за того, что в данных
# присутствуют выбросы, и также могут быть близкими к нулю, если зависимость между переменнами очень сильная, но не линейная.
# Также не стоит по большим значениям корреляции делать выводы о причинно следственных связях, они чаще всего будут неверными
# (если, например, на обе наших переменных оказывает влияние некоторая третья). Вот несколько примеров, с которыми стоит
# внимательно разобраться:
# 
# https://en.wikipedia.org/wiki/Anscombe%27s_quartet
# 
# http://www.tylervigen.com/
# 
# http://guessthecorrelation.com/
# 
# http://www.machinelearning.ru/wiki/images/e/e7/Psad_corr.pdf (слайды 1-6, особенно 5-ый!)
# 

# In[10]:


dc = data.corr()
idx = np.where(np.abs(dc) > 0.5)
indices = [(dc.index[x], dc.columns[y]) for x, y in zip(*idx) if x < y]
indices


# **Наблюдение:** есть сильная корреляция между некоторыми признаками.
# 
# **Вывод:** можно попробовать(!) преобразовать коррелированные признаки (выбросить часть? взять среднее?)

# **Наблюдение:** есть сильная корреляция между некоторыми признаками и целевой переменной
# 
# **Вывод:** ?

# **Наблюдение:** некоторые переменные очень слабо коррелированы с целевой переменной
# 
# **Вывод:** ?

# Посмотрим теперь на признак absences. По графикам видно, что есть небольшое количество точек, которые стоят обособленно -- это
# скорее всего выбросы (очень большие, нехарактерные значения признака). С ними нужно что-то сделать. Возможно стоит просто 
# выбросить эти наблюдения. Но также можно подумать в сторону того, чтобы разбить признак absences на группы (много прогуливал, 
# мало прогуливал, средне и т.д.) и тогда эти выбросы просто попадут в последнюю группу и могут оказаться полезными (а могут и 
# не оказаться).

# In[11]:


sns.jointplot('Talc', 'absences', data)
data.hist('absences', bins=50)


# **Наблюдение:** признак absences содержит выбросы
# 
# **Вывод:** выбросим эти объекты? добавим новый признак??
# 
# *feature engineering (добавление новых признаков) — это то, что обычно дает наибольший прирост в качестве*

# Попробуем теперь разобраться со странной полоской, возникающей в признаках G{1,2,3} (эта полоска соответствует G2 = 0 или G3 = 
# 0). То есть, она соответствует тем студентам, которым поставили 0 за какой-то экзамен. Странно, что они выбиваются из общей 
# зависимости и, скорее всего, для нас также являются выбросами. Попробуем посмотреть на данные и понять, что повляило на эту 
# оценку: может быть они много болели или может быть много прогуливали?

# In[12]:


pd.options.display.max_columns = 100


# In[13]:


data[data['G1'] == 0]


# In[14]:


data[data['G2'] == 0]


# In[15]:


data[data['G3'] == 0]


# Глядя на данные объяснить эту зависимость довольно сложно. Поэтому, скорее всего она вызвана какими-то причинами нам неизвестными (может быть этих студентов на самом деле отчислили, а потом автоматом ставили 0, или может быть у них что-то случилось в семье, или что-то другое). Так как эту зависимость мы не знаем, то для нас это чистые выбросы, которые опять же стоит удалить.

# **Наблюдение:** отклонение от зависимости на графиках G1, G2, G3. Явного объяснения из данных найти не удается.
# 
# **Вывод:** добавить признак? удалить выбросы?

# Стоит ли добавлять квадраты признаков? Глядя на графики, кажется, что наша целевая переменная не зависит квадратично от вещественных признаков. Но не стоит торопиться с выводами! Так как мы рассматриваем не только вещественные переменные, но и категориальные и бинарные, то реально, у нас может образоваться какое-то подмножество объектов (например, только девочки) от которых уже будет квадратичная зависимость. Поэтому, подобные выводы всегда стоит проверять на практике -- обучать модели и смотреть, насколько хорошо они работают.

# ### Посмотрим на категориальные переменные:

# In[16]:


figure, axes = plt.subplots(4, 4, figsize=(15, 15))
c = 0
for i in range(4):
    for j in range(4):
        sns.swarmplot(categorical_features[c], 'Talc', data=data, ax=axes[i, j])
        c += 1
        if c == len(categorical_features):
            break


# In[17]:


figure, axes = plt.subplots(4, 4, figsize=(15, 15))
c = 0
for i in range(4):
    for j in range(4):
        sns.violinplot(categorical_features[c], 'Talc', data=data, ax=axes[i, j])
        c += 1
        if c == len(categorical_features):
            break


# Как интерпретировать полученные графики можно посмотреть, например, вот здесь: https://stanford.edu/~mwaskom/software/seaborn/tutorial/categorical.html.
# 
# Глядя на графики видно, что в некоторых категориях у нас очень мало объектов (например, для больших значений возраста). Это означает, что мы, скорее всего, не сможем адекватно настроить веса для этих категорий и с ними что-то нужно сделать. Можно ли просто выбросить эти объекты? Правильный ответ -- нет, потому что объекты с такими значениями категорий могут оказаться в тестовой выборке. Но что-то нужно сделать, например, можно объединить несколько категорий в одну.

# **Наблюдение:** есть значения категорий для которых очень мало данных
# 
# **Вывод:** можно попробовать объединить категории?

# Можем попробовать проверить предоположение о том, что оценка за экзамен зависит от количества неудач в прошлом. Действительно, кажется, что некоторая зависимость наблюдается:

# In[18]:


sns.violinplot('failures', 'G3', data=data)


# Мы не будем делать никаких содержательных выводов из этого наблюдения -- это просто иллюстрация того, насколько удобно использовать визуализацию для проверки ваших гипотез о данных.

# ### Посмотрим на бинарные переменные:

# In[19]:


figure, axes = plt.subplots(4, 4, figsize=(15, 15))
c = 0
for i in range(4):
    for j in range(4):
        sns.swarmplot(binary_features[c], 'Talc', data=data, ax=axes[i, j])
        c += 1
        if c == len(binary_features):
            break


# In[20]:


figure, axes = plt.subplots(4, 4, figsize=(15, 15))
c = 0
for i in range(4):
    for j in range(4):
        sns.violinplot(binary_features[c], 'Talc', data=data, ax=axes[i, j])
        c += 1
        if c == len(binary_features):
            break


# **Наблюдение:** распределения для некоторых бинарных признаков выглядят очень похожи
# 
# **Вывод:** удалим эти признаки? Нет!

# Глядя, например, на график признака activities можно заметить, что у нас практически одинаковое распределение на целевую переменную для обоих возможных значений признака. Значит ли это, что мы можем его удалить? Нет! Он может входить в содержательные зависимости вместе с другими признаками (например, если мы посмотрим на наличие активностей только у девочек, то зависимость целевой переменной теоретически может появиться). Поэтому перед тем как делать выводы стоит попробовать убедиться, что таких зависимостей нет:

# In[21]:


figure, axes = plt.subplots(4, 4, figsize=(15, 15))
c = 0
for i in range(4):
    for j in range(4):
        sns.swarmplot(binary_features[c], 'Talc', hue='activities', data=data, ax=axes[i, j])
        c += 1
        if c == len(binary_features):
            break


# In[22]:


figure, axes = plt.subplots(4, 4, figsize=(15, 15))
c = 0
for i in range(4):
    for j in range(4):
        sns.violinplot(binary_features[c], 'Talc', hue='activities', data=data, ax=axes[i, j])
        c += 1
        if c == len(binary_features):
            break


# In[23]:


figure, axes = plt.subplots(4, 4, figsize=(15, 15))
c = 0
for i in range(4):
    for j in range(4):
        sns.violinplot(categorical_features[c], 'Talc', hue='activities', split=True, data=data, ax=axes[i, j])
        c += 1
        if c == len(categorical_features):
            break


# Вроде бы подобных зависимостей не наблюдается (за исключением нескольких случаев, когда у нас просто очень мало данных). Но все равно выбрасывать этот признак стоит только после экспериментальной проверка, что он бесполезен, потому что мы не учитывали, например, тройные взаимодействия.

# ### Преобразуем данные

# In[24]:


pd.get_dummies(data).head()


# In[25]:


for feature in categorical_features + binary_features:
    data[feature] = data[feature].astype('category')


# In[26]:


pd.get_dummies(data).head()


# In[27]:


pd.get_dummies(data).shape


# Обратите внимание, что правильный способ преобразовать данные из категориальных в вещественные -- это не просто применить функцию get_dummies, но и предварительно сказать pandas, какие из наших переменных реально категориальные (иначе признаки, кодируемые цифрами так ими и останутся).

# In[28]:


data_transformed = pd.get_dummies(data)
# преобразуем данные после анализа
# data_transformed = ...


# In[29]:


num_data = data_transformed.drop(targets, axis=1).as_matrix()
num_targets = data_transformed['Talc'].as_matrix()


# In[30]:


np.min(num_data, axis=0), np.max(num_data, axis=0)


# **Наблюдение:** диапазоны признаков различаются
# 
# **Вывод:** признаки стоит нормировать. Как нормируем второй признак?

# Признаки практически всегда стоит отнормировать. Для обычной линейной регрессии это совершенно не обязательно, а вот если вы используете регуляризаторы, то это становится просто необходимо (почему?) Но будьте внимательны -- если у нас минимальное значение признаков на обучающей выборке не равно реальному минимальному значению (второй признак), то не стоит его бездумно вычитать, ведь мы можем получить меньшие значения на объектах тестовой выборки и все будет не очень хорошо!

# In[31]:


num_data = num_data / np.max(num_data, axis=0)


# In[32]:


from sklearn.cross_validation import train_test_split


# In[33]:


X_train, X_test, y_train, y_test = train_test_split(num_data, num_targets, test_size=0.2, random_state=42)


# In[34]:


from sklearn.linear_model import LogisticRegression
from sklearn.linear_model import LinearRegression


# Что используем?

# In[35]:


from sklearn.linear_model import LogisticRegression
from sklearn.linear_model import LinearRegression
from sklearn.linear_model import Ridge
from sklearn.linear_model import Lasso
from sklearn.ensemble import RandomForestClassifier
from sklearn.naive_bayes import GaussianNB
from sklearn.neighbors import KNeighborsClassifier
from sklearn.ensemble import GradientBoostingClassifier
# ...