Ml3a linear practice2

Линейная регрессия Практика 2¶

Сравнение подходов¶

• Сравним матричные операции с градиентном спуске
• Посмотрим как используя линейные модели моделировать нелинейную взаимосвязь

Импортируем модули¶

In [24]:

import numpy as np
import pandas as pd
from sklearn.datasets import load_boston
from sklearn.metrics import mean_squared_error, f1_score, accuracy_score, roc_curve, roc_auc_score
from sklearn.model_selection import train_test_split

from matplotlib import pyplot as plt

import numpy as np import pandas as pd from sklearn.datasets import load_boston from sklearn.metrics import mean_squared_error, f1_score, accuracy_score, roc_curve, roc_auc_score from sklearn.model_selection import train_test_split from matplotlib import pyplot as plt

Когда использовать матричные операции¶

вместо градиентного впуска

In [25]:

def print_regression_metrics(y_true, y_pred):
    mse = mean_squared_error(y_true, y_pred)
    rmse = np.sqrt(mse)
    print(f'MSE = {mse:.2f}, RMSE = {rmse:.2f}')
    
def prepare_boston_data():
    data = load_boston()
    X, y = data['data'], data['target']
    # Нормализовать даннные с помощью стандартной нормализации
    X = (X - X.mean(axis=0)) / X.std(axis=0)
    # Добавить фиктивный столбец единиц (bias линейной модели)
    X = np.hstack([np.ones(X.shape[0])[:, np.newaxis], X])
    
    return X, y

def print_regression_metrics(y_true, y_pred): mse = mean_squared_error(y_true, y_pred) rmse = np.sqrt(mse) print(f'MSE = {mse:.2f}, RMSE = {rmse:.2f}') def prepare_boston_data(): data = load_boston() X, y = data['data'], data['target'] # Нормализовать даннные с помощью стандартной нормализации X = (X - X.mean(axis=0)) / X.std(axis=0) # Добавить фиктивный столбец единиц (bias линейной модели) X = np.hstack([np.ones(X.shape[0])[:, np.newaxis], X]) return X, y

Прежде чем начать, обернем написанную нами

• линейную регрессию методом матричны операций в класс:

In [26]:

class LinRegAlgebra():
    def __init__(self):
        self.theta = None
    
    def fit(self, X, y):
        self.theta = np.linalg.inv(X.T.dot(X)).dot(X.transpose()).dot(y)
    
    def predict(self, X):
        return X.dot(self.theta)

class LinRegAlgebra(): def __init__(self): self.theta = None def fit(self, X, y): self.theta = np.linalg.inv(X.T.dot(X)).dot(X.transpose()).dot(y) def predict(self, X): return X.dot(self.theta)

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

In [27]:

X, y = prepare_boston_data()

X, y = prepare_boston_data()

In [28]:

linreg_alg = LinRegAlgebra()
linreg_alg.fit(X, y)
y_pred = linreg_alg.predict(X)

# Посчитать значение ошибок MSE и RMSE для тренировочных данных
print_regression_metrics(y, y_pred)

linreg_alg = LinRegAlgebra() linreg_alg.fit(X, y) y_pred = linreg_alg.predict(X) # Посчитать значение ошибок MSE и RMSE для тренировочных данных print_regression_metrics(y, y_pred)

MSE = 21.89, RMSE = 4.68

In [29]:

class RegOptimizer():
    def __init__(self, alpha, n_iters):
        self.theta = None
        self._alpha = alpha
        self._n_iters = n_iters
    
    def gradient_step(self, theta, theta_grad):
        return theta - self._alpha * theta_grad
    
    def grad_func(self, X, y, theta):
        raise NotImplementedError()

    def optimize(self, X, y, start_theta, n_iters):
        theta = start_theta.copy()

        for i in range(n_iters):
            theta_grad = self.grad_func(X, y, theta)
            theta = self.gradient_step(theta, theta_grad)

        return theta
    
    def fit(self, X, y):
        m = X.shape[1]
        start_theta = np.ones(m)
        self.theta = self.optimize(X, y, start_theta, self._n_iters)
        
    def predict(self, X):
        raise NotImplementedError()

class RegOptimizer(): def __init__(self, alpha, n_iters): self.theta = None self._alpha = alpha self._n_iters = n_iters def gradient_step(self, theta, theta_grad): return theta - self._alpha * theta_grad def grad_func(self, X, y, theta): raise NotImplementedError() def optimize(self, X, y, start_theta, n_iters): theta = start_theta.copy() for i in range(n_iters): theta_grad = self.grad_func(X, y, theta) theta = self.gradient_step(theta, theta_grad) return theta def fit(self, X, y): m = X.shape[1] start_theta = np.ones(m) self.theta = self.optimize(X, y, start_theta, self._n_iters) def predict(self, X): raise NotImplementedError()

In [30]:

class LinReg(RegOptimizer):

    def grad_func(self, X, y, theta):
        n = X.shape[0]
        grad = 1. / n * X.T.dot(X.dot(theta) - y)
        return grad
    
    def predict(self, X):
        if self.theta is None:
            raise Exception('You should train the model first')
        
        y_pred = X.dot(self.theta)
        
        return y_pred

class LinReg(RegOptimizer): def grad_func(self, X, y, theta): n = X.shape[0] grad = 1. / n * X.T.dot(X.dot(theta) - y) return grad def predict(self, X): if self.theta is None: raise Exception('You should train the model first') y_pred = X.dot(self.theta) return y_pred

In [31]:

linreg_crit = LinReg(0.2,1000)
linreg_crit.fit(X, y)
y_pred = linreg_crit.predict(X)

# Посчитать значение ошибок MSE и RMSE для тренировочных данных
print_regression_metrics(y, y_pred)

linreg_crit = LinReg(0.2,1000) linreg_crit.fit(X, y) y_pred = linreg_crit.predict(X) # Посчитать значение ошибок MSE и RMSE для тренировочных данных print_regression_metrics(y, y_pred)

MSE = 21.89, RMSE = 4.68

Теперь измерим скорость выполнения

In [32]:

%timeit linreg_alg.fit(X, y)

%timeit linreg_alg.fit(X, y)

165 µs ± 1.29 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

In [33]:

%timeit linreg_crit.fit(X, y)

%timeit linreg_crit.fit(X, y)

15.4 ms ± 79.8 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

In [34]:

linreg_crit.fit(X, y)

linreg_crit.fit(X, y)

Когда какой метод использовать¶

Реализация на матричных операциях опережает реализацию на градиентном спуске в сотни раз

• Но всегда ли это так и какие подводные камни могут быть?
• Ниже приведен набор случаев, при которых версия с градентным спуском предпочтительнее:

1. Градиентный спуск работает быстрее в случае матриц с большим количеством признаков. Основная по сложности операция — нахождение обратной матрицы $(X^T X)^{-1}$.
2. Нахождение обратной матрицы может также потребовать больше оперативной памяти
3. Матричные операции могут также проигрывать и в случае небольших объемов данных, но при плохой параллельной реализации или недостаточных ресурсах.
4. Градиентный спуск может быть усовершенствован до так называемого стохастического градиентного спуска, (данные для оптимизации подгружаются небольшими наборами), что уменьшает требования по памяти.
5. В некоторых случаях (например, в случае линейно-зависимых строк) алгебраический способ решения не будет работать совсем в виду невозможности найти обратную матрицу.

Превращение линейной модели в нелинейную¶

• Нелинейные зависимости в данных встречаются намного чаще линейных
• На самом деле простейшая линейная регрессия способна обнаруживать нелинейные зависимости
• Для этого необходимо рассмотреть дополнительные признаки, полученные из исходных применением различных нелинейных функций
• Возьмем уже знакомый датасет с ценами на квартиры в Бостоне и последовательно станем применять различные функции к исходным признакам:

In [35]:

# Boston Data. Attribute Information (in order):
#     - CRIM     per capita crime rate by town
#     - ZN       proportion of residential land zoned for lots over 25,000 sq.ft.
#     - INDUS    proportion of non-retail business acres per town
#     - CHAS     Charles River dummy variable (= 1 if tract bounds river; 0 otherwise)
#     - NOX      nitric oxides concentration (parts per 10 million)
#     - RM       average number of rooms per dwelling
#     - AGE      proportion of owner-occupied units built prior to 1940
#     - DIS      weighted distances to five Boston employment centres
#     - RAD      index of accessibility to radial highways
#     - TAX      full-value property-tax rate per `$10000`
#     - PTRATIO  pupil-teacher ratio by town
#     - B        1000(Bk - 0.63)^2 where Bk is the proportion of blacks by town
#     - LSTAT    % lower status of the population
#     - MEDV     Median value of owner-occupied homes in $1000's

# Boston Data. Attribute Information (in order): # - CRIM per capita crime rate by town # - ZN proportion of residential land zoned for lots over 25,000 sq.ft. # - INDUS proportion of non-retail business acres per town # - CHAS Charles River dummy variable (= 1 if tract bounds river; 0 otherwise) # - NOX nitric oxides concentration (parts per 10 million) # - RM average number of rooms per dwelling # - AGE proportion of owner-occupied units built prior to 1940 # - DIS weighted distances to five Boston employment centres # - RAD index of accessibility to radial highways # - TAX full-value property-tax rate per `$10000` # - PTRATIO pupil-teacher ratio by town # - B 1000(Bk - 0.63)^2 where Bk is the proportion of blacks by town # - LSTAT % lower status of the population # - MEDV Median value of owner-occupied homes in $1000's

In [36]:

def prepare_boston_data_new():
    data = load_boston()
    X, y = data['data'], data['target']
    
    X = np.hstack([X, np.sqrt(X[:, 5:6]), X[:, 6:7] ** 3])
    
    # Нормализовать даннные с помощью стандартной нормализации
    X = (X - X.mean(axis=0)) / X.std(axis=0)
    # Добавить фиктивный столбец единиц (bias линейной модели)
    X = np.hstack([np.ones(X.shape[0])[:, np.newaxis], X])
    
    return X, y

def prepare_boston_data_new(): data = load_boston() X, y = data['data'], data['target'] X = np.hstack([X, np.sqrt(X[:, 5:6]), X[:, 6:7] ** 3]) # Нормализовать даннные с помощью стандартной нормализации X = (X - X.mean(axis=0)) / X.std(axis=0) # Добавить фиктивный столбец единиц (bias линейной модели) X = np.hstack([np.ones(X.shape[0])[:, np.newaxis], X]) return X, y

Создадим несколько нелинейных признаков¶

Мы добавили два новых признака:

1. Взяли корень из признака RM (среднее число комнат на сожителя)
2. Возвели в куб значения признака AGE

Это только два примера. Всевозможных комбинаций признаков и примененных к ним функций неограниченное количество.

In [37]:

def train_validate(X, y):
    # Разбить данные на train/valid
    X_train, X_valid, y_train, y_valid = train_test_split(X, y, test_size=0.2, shuffle=True, random_state=1)

    # Создать и обучить линейную регрессию
    linreg_alg = LinRegAlgebra()
    linreg_alg.fit(X_train, y_train)

    # Сделать предсказания по валидционной выборке
    y_pred = linreg_alg.predict(X_valid)

    # Посчитать значение ошибок MSE и RMSE для валидационных данных
    print_regression_metrics(y_valid, y_pred)

def train_validate(X, y): # Разбить данные на train/valid X_train, X_valid, y_train, y_valid = train_test_split(X, y, test_size=0.2, shuffle=True, random_state=1) # Создать и обучить линейную регрессию linreg_alg = LinRegAlgebra() linreg_alg.fit(X_train, y_train) # Сделать предсказания по валидционной выборке y_pred = linreg_alg.predict(X_valid) # Посчитать значение ошибок MSE и RMSE для валидационных данных print_regression_metrics(y_valid, y_pred)

In [38]:

# Подготовить данные без модификации признаков
X, y = prepare_boston_data()
# Провести эксперимент
train_validate(X, y)

# Подготовить данные без модификации признаков X, y = prepare_boston_data() # Провести эксперимент train_validate(X, y)

MSE = 23.38, RMSE = 4.84

In [39]:

# Подготовить данные без модификации признаков
X, y = prepare_boston_data_new()
# Провести эксперимент
train_validate(X, y)

# Подготовить данные без модификации признаков X, y = prepare_boston_data_new() # Провести эксперимент train_validate(X, y)

MSE = 14.28, RMSE = 3.78

Как видно из результатов, мы добились улучшения точности предсказаний на 40%, всего лишь добавив пару нелинейных признаков в имеющимся. Можете поиграть с признаками и еще улучшить результат.

Задание¶

Задание 3.7.1¶

• Сделайте для градиентного спуска остановку алгоритма, сли максимальное из абсолютных значений компонент градиента становится меньше 0.01.
• Сравните скорость обучения градиентным спуском и матричными операциями
• Для градиентного спуска установите alpha = 0.2.
• На какой итерации останавливается градиентный спуск?

In [40]:

# Функция для гдариентного спуска
class RegOptimizer():
    
    def __init__(self, alpha, n_iters, limiter=None):
        self.theta = None
        self._alpha = alpha
        self._n_iters = n_iters
        self._limiter = limiter
    
    def gradient_step(self, theta, theta_grad):
        return theta - self._alpha * theta_grad
    
    def grad_func(self, X, y, theta):
        raise NotImplementedError()

    def optimize(self, X, y, start_theta, n_iters):

        theta = start_theta.copy()
        for i in range(n_iters):
            theta_grad = self.grad_func(X, y, theta)
            theta = self.gradient_step(theta, theta_grad)
            if(self._limiter != None):
                if(max(theta_grad)<self._limiter):
                    print(f'max theta_grad reached: {i}')
                    break
        
        return theta
    
    def fit(self, X, y):
        m = X.shape[1]
        start_theta = np.ones(m)
        self.theta = self.optimize(X, y, start_theta, self._n_iters)
        
    def predict(self, X):
        raise NotImplementedError()

# Линейная регрессия используя градиентный спуск
class LinReg(RegOptimizer):
    
    def grad_func(self, X, y, theta):
        n = X.shape[0]
        grad = 1. / n * X.transpose().dot(X.dot(theta) - y)
        return grad
    
    def predict(self, X):
        if self.theta is None:
            raise Exception('You should train the model first')
        
        y_pred = X.dot(self.theta)
        
        return y_pred

# Функция для гдариентного спуска class RegOptimizer(): def __init__(self, alpha, n_iters, limiter=None): self.theta = None self._alpha = alpha self._n_iters = n_iters self._limiter = limiter def gradient_step(self, theta, theta_grad): return theta - self._alpha * theta_grad def grad_func(self, X, y, theta): raise NotImplementedError() def optimize(self, X, y, start_theta, n_iters): theta = start_theta.copy() for i in range(n_iters): theta_grad = self.grad_func(X, y, theta) theta = self.gradient_step(theta, theta_grad) if(self._limiter != None): if(max(theta_grad)

In [41]:

# Для задания 3.7.1 
def train_validate_limiter(X, y, limiter=0.01):
    
    # Разбить данные на train/valid
    X_train, X_valid, y_train, y_valid = train_test_split(X, y, 
                                                          test_size=0.2, 
                                                          shuffle=True, 
                                                          random_state=1)

    # Создать и обучить линейную регрессию
    linreg_alg = LinReg(0.2,1000,limiter)
    linreg_alg.fit(X_train, y_train)

    # Сделать предсказания по валидционной выборке
    y_pred = linreg_alg.predict(X_valid)

    # Посчитать значение ошибок MSE и RMSE для валидационных данных
    print_regression_metrics(y_valid, y_pred)

X1,y1 = prepare_boston_data()
train_validate_limiter(X1,y1,0.01)

# Для задания 3.7.1 def train_validate_limiter(X, y, limiter=0.01): # Разбить данные на train/valid X_train, X_valid, y_train, y_valid = train_test_split(X, y, test_size=0.2, shuffle=True, random_state=1) # Создать и обучить линейную регрессию linreg_alg = LinReg(0.2,1000,limiter) linreg_alg.fit(X_train, y_train) # Сделать предсказания по валидционной выборке y_pred = linreg_alg.predict(X_valid) # Посчитать значение ошибок MSE и RMSE для валидационных данных print_regression_metrics(y_valid, y_pred) X1,y1 = prepare_boston_data() train_validate_limiter(X1,y1,0.01)

max theta_grad reached: 217
MSE = 23.40, RMSE = 4.84

In [42]:

data = load_boston()
data.feature_names
df = pd.DataFrame(data.data,columns=data.feature_names)
df.head()

data = load_boston() data.feature_names df = pd.DataFrame(data.data,columns=data.feature_names) df.head()

Out[42]:

	CRIM	ZN	INDUS	CHAS	NOX	RM	AGE	DIS	RAD	TAX	PTRATIO	B	LSTAT
0	0.00632	18.0	2.31	0.0	0.538	6.575	65.2	4.0900	1.0	296.0	15.3	396.90	4.98
1	0.02731	0.0	7.07	0.0	0.469	6.421	78.9	4.9671	2.0	242.0	17.8	396.90	9.14
2	0.02729	0.0	7.07	0.0	0.469	7.185	61.1	4.9671	2.0	242.0	17.8	392.83	4.03
3	0.03237	0.0	2.18	0.0	0.458	6.998	45.8	6.0622	3.0	222.0	18.7	394.63	2.94
4	0.06905	0.0	2.18	0.0	0.458	7.147	54.2	6.0622	3.0	222.0	18.7	396.90	5.33

In [43]:

list(df.columns).index('DIS')

list(df.columns).index('DIS')

Out[43]:

7

Задание 3.7.2¶

• Добавьте к признакам нелинейной модели квадрат признака DIS и переобучите модель
• Какой получился RMSE?
• Подсказка: используйте написанную нами линейную регрессию методом матричных операций.

In [44]:

# Предобработка данных 

def prepare_boston_data_new():
    
    data = load_boston()
    X, y = data['data'], data['target']
    
    # Добавяем нелинейные признаки
    X = np.hstack([X, 
                   np.sqrt(X[:, 5:6]), 
                   X[:, 6:7] ** 3])
    
    # Нормализовать даннные с помощью 
    # стандартной нормализации
    model = StandardScaler()
    X = model.fit_transform(X)
    
    # Добавить фиктивный столбец единиц 
    # (bias линейной модели)
    X = np.hstack([np.ones(X.shape[0])[:, None], X])
    
    return X, y

# Разбтваем выборку на train,test
# Тренируем линейную регрессию (на основе матричных умнажении)

def train_validate_algebra(X, y):
    
    # Разбить данные на train/valid
    X_train, X_valid, y_train, y_valid = train_test_split(X, y, 
                                                          test_size=0.2, 
                                                          shuffle=True, 
                                                          random_state=1)

    # Создать и обучить линейную регрессию
    linreg_alg = LinRegAlgebra()
    linreg_alg.fit(X_train, y_train)

    # Сделать предсказания по валидционной выборке
    y_pred = linreg_alg.predict(X_valid)

    # Посчитать значение ошибок MSE и RMSE для валидационных данных
    print_regression_metrics(y_valid, y_pred)
    
# Разбтваем выборку на train,test
# Тренируем линейную регрессию (градиентный спуск)
    
def train_validate(X, y):
    
    # Разбить данные на train/valid
    X_train, X_valid, y_train, y_valid = train_test_split(X, y, test_size=0.2,shuffle=True,random_state=1)

    # Создать и обучить линейную регрессию
    linreg_alg = LinReg(0.2,1000)
    linreg_alg.fit(X_train, y_train)

    # Сделать предсказания по валидционной выборке
    y_pred = linreg_alg.predict(X_valid)

# Посчитать значение ошибок MSE и RMSE для валидационных данных
    print_regression_metrics(y_valid, y_pred)

# Предобработка данных def prepare_boston_data_new(): data = load_boston() X, y = data['data'], data['target'] # Добавяем нелинейные признаки X = np.hstack([X, np.sqrt(X[:, 5:6]), X[:, 6:7] ** 3]) # Нормализовать даннные с помощью # стандартной нормализации model = StandardScaler() X = model.fit_transform(X) # Добавить фиктивный столбец единиц # (bias линейной модели) X = np.hstack([np.ones(X.shape[0])[:, None], X]) return X, y # Разбтваем выборку на train,test # Тренируем линейную регрессию (на основе матричных умнажении) def train_validate_algebra(X, y): # Разбить данные на train/valid X_train, X_valid, y_train, y_valid = train_test_split(X, y, test_size=0.2, shuffle=True, random_state=1) # Создать и обучить линейную регрессию linreg_alg = LinRegAlgebra() linreg_alg.fit(X_train, y_train) # Сделать предсказания по валидционной выборке y_pred = linreg_alg.predict(X_valid) # Посчитать значение ошибок MSE и RMSE для валидационных данных print_regression_metrics(y_valid, y_pred) # Разбтваем выборку на train,test # Тренируем линейную регрессию (градиентный спуск) def train_validate(X, y): # Разбить данные на train/valid X_train, X_valid, y_train, y_valid = train_test_split(X, y, test_size=0.2,shuffle=True,random_state=1) # Создать и обучить линейную регрессию linreg_alg = LinReg(0.2,1000) linreg_alg.fit(X_train, y_train) # Сделать предсказания по валидционной выборке y_pred = linreg_alg.predict(X_valid) # Посчитать значение ошибок MSE и RMSE для валидационных данных print_regression_metrics(y_valid, y_pred)

In [45]:

from sklearn.preprocessing import StandardScaler

# Предобработка данных 
def prepare_boston_data_new_mod():
    
    data = load_boston()
    X, y = data['data'], data['target']
    
    # Добавяем нелинейные признаки
    X = np.hstack([X, 
                   np.sqrt(X[:, 5:6]), 
                   X[:, 6:7] ** 3,
                   X[:, 7:8] ** 2])
                   
# Нормализовать даннные с помощью 
    # стандартной нормализации
    model = StandardScaler()
    X = model.fit_transform(X)
    
    # Добавить фиктивный столбец единиц 
    # (bias линейной модели)
    X = np.hstack([np.ones(X.shape[0])[:, None], X])
    
    return X, y

X,y = prepare_boston_data_new_mod()
train_validate_algebra(X,y)

# 3.69

from sklearn.preprocessing import StandardScaler # Предобработка данных def prepare_boston_data_new_mod(): data = load_boston() X, y = data['data'], data['target'] # Добавяем нелинейные признаки X = np.hstack([X, np.sqrt(X[:, 5:6]), X[:, 6:7] ** 3, X[:, 7:8] ** 2]) # Нормализовать даннные с помощью # стандартной нормализации model = StandardScaler() X = model.fit_transform(X) # Добавить фиктивный столбец единиц # (bias линейной модели) X = np.hstack([np.ones(X.shape[0])[:, None], X]) return X, y X,y = prepare_boston_data_new_mod() train_validate_algebra(X,y) # 3.69

MSE = 13.59, RMSE = 3.69

Задание 3.7.3¶

Уберите нормализацию и оставьте добавленные признаки на основе RM и AGE, Какой получился RMSE?

In [46]:

# Предобработка данных 
def prepare_boston_data_new_mod():
    
    data = load_boston()
    X, y = data['data'], data['target']
    
    # Добавяем нелинейные признаки
    X = np.hstack([X, 
                   np.sqrt(X[:, 5:6]), 
                   X[:, 6:7] ** 3])
                   
    # Нормализовать даннные с помощью 
    # стандартной нормализации
#     model = StandardScaler()
#     X = model.fit_transform(X)
    
    # Добавить фиктивный столбец единиц 
    # (bias линейной модели)
    X = np.hstack([np.ones(X.shape[0])[:, None], X])
    
    return X, y

X,y = prepare_boston_data_new_mod()
train_validate_algebra(X,y)

# Предобработка данных def prepare_boston_data_new_mod(): data = load_boston() X, y = data['data'], data['target'] # Добавяем нелинейные признаки X = np.hstack([X, np.sqrt(X[:, 5:6]), X[:, 6:7] ** 3]) # Нормализовать даннные с помощью # стандартной нормализации # model = StandardScaler() # X = model.fit_transform(X) # Добавить фиктивный столбец единиц # (bias линейной модели) X = np.hstack([np.ones(X.shape[0])[:, None], X]) return X, y X,y = prepare_boston_data_new_mod() train_validate_algebra(X,y)

MSE = 14.28, RMSE = 3.78