Некоторые хитрости в ML

\(Ŷ = \beta{}_{0} + \beta{}_{1}*X\)

Где

\(\beta{}_{1} = \frac{sd_{x}}{sd_{y}} * r_{xy}\),

\(\beta{}_{0} = \hat{Y} - \beta_{1} * X\).

library(MASS)
lm.fit <- lm(medv~lstat, data=Boston)
plot(Boston$lstat, Boston$medv)
abline(lm.fit, lwd=2, col="red")

library(MASS)
lm.fit <- lm(medv~lstat+I(lstat^2), data=Boston)
plot(Boston$lstat, Boston$medv)
curve(predict(lm.fit, newdata=data.frame(lstat=x)), add=T, col="red", lwd=2)

1. Построить scatter plot и посмотреть, что данные подчиняются линейному закону.
2. Посчитать \(R^{2}\) и выяснить подчиняются ли данные линейной корреляции. \(R^{2}\) описывает какой процент данных (от 0 до 1) описывается линейной регрессией.
3. Проверить гомоскедастичность данных (homoscedasticity vs heteroscedasticity)
4. Высчитать \(\beta_{0}\) и \(\beta_{1}\).
5. Проверить нулевую гипотезу \(H_{0}: \beta_{1} = 0\)

1. Посчитать F-статистику для данных и посмотреть вероятность нулевой гипотезы. Нулевая гипотеза для множественной регресси это: \(H_{0}: \beta_{1} = \beta_{2} = ... = \beta_{n} = 0\) против альтернативной гипотезы: \(H_{a}: \nexists{} \beta_{j} \neq{} 0\). Посчитав F-значение для наших данных, необходимо вычислить $p$-вероятность такого значения в F-распределении.

   
   lm.fit <- lm(poverty~.-state, data=df)
   summary(lm.fit)$fstatistic
   

     value   numdf   dendf
   20.5849  4.0000 46.0000
   
   

2. Отобрать только значимые переменные. В некоторых случаях число переменных (features) \(p\) может быть больше или равно числу измерений \(n\). В таком случае, линейная регрессия будет работать медленно и некорректно. Необходимо выбрать только те переменные, которые действительно описывают нашу регрессию. Есть много вариантов такого выбора: Adjusted R², Forward selection, Backward selection, Mixed selection.

3. Проверка качества нашей модели Необходимо проверить насколько хорошо наша модель описывает закон, которому подчиняются данные. Один из способов это проверить — посчитать \(R^{2}\). Так же мы можем проверить какие переменные реально влияют на распределение данных. Если последовательно считать \(R^{2}\) для \(n+1\) переменной, то можно увидить значительно или нет изменяется значение \(R^{2}\). Значение \(R^{2}\) будет расти в любом случае при добавлении новой переменной, но иногда это тысячные, значит что переменная скорее мешает нашей модели чем объясняет её.

   
   lm.fit <- lm(poverty~.-state, data=df)
   summary(lm.fit)[c("r.squared", "adj.r.squared")]
   

   $r.squared
   [1] 0.6415759
   
   $adj.r.squared
   [1] 0.6104086
   
   
   

4. Посчитать доверительный интервал для наших предсказанных значений. В таком случае, предсказанное новое значение \(\hat{Y} = \hat{f}(X)\) будет представлено интервалом значений.

   
   lm.fit <- lm(poverty~.-state, data=df)
   confint(lm.fit)
   

                     2.5 %      97.5 %
   (Intercept)  41.1343729 91.81867786
   metro_res    -0.0956751 -0.01697061
   white        -0.1146987  0.01840974
   hs_grad      -0.7658790 -0.34354762
   female_house -0.4391957  0.54028167
   
   

1. Линейная зависимость переменных. Для проверки, мы можем построить residuals vs. fitted values график.

   
   lm.fit <- lm(poverty~.-state, data=df)
   plot(predict(lm.fit), residuals(lm.fit))
   

   

2. Нормальное распределение остатков

3. Гомоскедастичность остатков

   
   fit_var <- lm(poverty~metro_res+white+hs_grad, df)
   fit_res <- fit_var$residuals
   
   par(mfrow=c(1,2))
   plot(
     fit_res,
     main = 'Анализ остатков',
     xlab = 'Предсказанные значения',
     ylab = 'Остатки'
   )
   points(fit_res, col = 'black', lwd=5, cex=.5)
   abline(h = 0, col = 'red', lwd=2, lty=2)
   qqnorm(fit_res, lwd=2)
   qqline(fit_res, col='red', lwd=2)
   par(mfrow=c(1,1))
   

   

4. Выборсы. Выбросы (измерения, имеющие значительные отличия по \(y_{i}\)) не значительно, но влияют на линейную регрессию. В частности, они уменьшают значение \(R^{2}\) Для того чтобы обнаружить выборсы, полезно воспользоваться studentized residuals графиком. Все точки, которые имеют значение studentized residual больше 3 — это выбросы.

   
   lm.fit <- lm(poverty~metro_res+white+hs_grad, df)
   plot(residuals(lm.fit), rstudent(lm.fit))
   

   

5. High Leverage Points (плечо?) В отличие от выбросов, high leverage points это точки, которые имеют значительно отличие по параметру X (в одной или нескольких плоскостях). Эти точки значительно влияют на регрессионную модель. Чтобы обнаружить такие точки, полезно построить график leverage vs studentized residuals

   
   fit_var <- lm(poverty~metro_res+white+hs_grad, df)
   par(mfrow=c(2, 2))
   plot(fit_var)
   

   

6. Проверка на мультиколлинеарность

   library("GGally")
   
   
   ggpairs(df[,-1])
   

   

   library(psych)
   
   
   pairs.panels(
     df[, -1],
     method = "pearson",
     hist.col = "cornflowerblue",
     density = T,
     ellipses = F
   )
   

   

   import pandas as pd
   import matplotlib.pyplot as plt
   
   DATA = pd.read_csv(
       'http://d396qusza40orc.cloudfront.net/statistics/lec_resources/states.csv'
   )
   
   AXES = pd.plotting.scatter_matrix(
       DATA, figsize=(6,6), diagonal='kde', grid=True
   )
   
   CORR = DATA.corr().values
   for i, j in zip(*plt.np.triu_indices_from(AXES, k=1)):
       AXES[i, j].annotate(
           '%.3f' % CORR[i, j],
           (0.8, 0.8),
           xycoords='axes fraction',
           ha='center',
           va='center'
       )
   
   figpath = 'Py.png'
   plt.savefig(figpath)
   return figpath
   

   

   Другой способ проверки на мультиколлинеарность: высчитать VIF — variance inflation factor. Самое маленькое значение VIF = 1, это означает полное отсутствие коллинеарности, значение VIF от 5 до 10 сигнализирует о проблеме с коллинеарностью данных. Преимущество VIF в отличие от матрицы коллинеарности в том, что он считает коэффециент коллинеарности не только между двумя предикторами.

   library(car)
   
   
   fit_var <- lm(poverty~metro_res+white+hs_grad, df)
   vif(fit_var)
   

   metro_res     white   hs_grad
    1.146522  1.215062  1.072928
   
   

7. Нормальное распределение переменных (желательно)

Используя R (или воспользоваться Rserve + pyRserve: https://www.rforge.net/Rserve/doc.html) можно сделать очень быстро и просто

library(MASS)
model <- lm(medv ~ ., data=Boston)
par(mfrow=c(2,2))
plot(model)

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

from sklearn.datasets import load_boston
from statsmodels.graphics.gofplots import ProbPlot

plt.style.use('seaborn')  # pretty matplotlib plots
plt.rc('font', size=14)
plt.rc('figure', titlesize=18)
plt.rc('axes', labelsize=15)
plt.rc('axes', titlesize=18)

boston = load_boston()

X = pd.DataFrame(boston.data, columns=boston.feature_names)
y = pd.DataFrame(boston.target)

# generate OLS model
model = sm.OLS(y, sm.add_constant(X))
model_fit = model.fit()

# create dataframe from X, y for easier plot handling
dataframe = pd.concat([X, y], axis=1)

# model values
model_fitted_y = model_fit.fittedvalues
# model residuals
model_residuals = model_fit.resid
# normalized residuals
model_norm_residuals = model_fit.get_influence().resid_studentized_internal
# absolute squared normalized residuals
model_norm_residuals_abs_sqrt = np.sqrt(np.abs(model_norm_residuals))
# absolute residuals
model_abs_resid = np.abs(model_residuals)
# leverage, from statsmodels internals
model_leverage = model_fit.get_influence().hat_matrix_diag
# cook's distance, from statsmodels internals
model_cooks = model_fit.get_influence().cooks_distance[0]

Ng предлагает использовать такой порядок \(\alpha\): \(0.001 \dots{} 0.003 \dots{} 0.01 \dots{} 0.03 \dots{} 0.1 \dots{} 0.3 \dots{} 1\)

\(\mathbf{n}\) = 1000 уже стоит использовать Gradient Descent.

используем слишком простую модель, в итоге получаем плохой результат для тренировочных данных и для тестовых данных.

используем слишком сложную модель, в итоге получаем идеальный результат для тренировочных данных (квадрат ошибок вплоть до 0), но на тестовых данных всё будет плохо, так как модель заточена только под конкретный набор тренировочных данных.

Добавляем слагаемое к \(RSS + \lambda{} * \sum_{j=1}^{p}(\theta{}_j^2)\) для всех \(\theta{} \in{} 1,\dots{},j\). Таким образом мы уменьшаем значения \(\theta{}\) даже для очень сложных многочленов, чтобы \(J(\theta{})\) было минимальное. Параметр \(\lambda{}\) стоит брать поменьше, но не \(0\), иначе это просто выключает регуляризацию.

Для classification-задач можно сделать такую матрицу значений

в scikit это можно сделать следующим образом:

from sklearn.metrics import confusion_matrix
from sklearn.metrics import classification_report

# Do some classifications

confusion_matrix(y_tested, y_predicted)  # -> [[52, 7], [3, 112]] for example
classification_report(y_tested, y_predicted) # -> table with columns [precision, recall, f1-score, support]