Линейная регрессия и результаты в области здравоохранения

Набор данных о диабете

Мы будем использовать набор данных, полученный от пациентов с диабетом. Данные состоят из 442 образцов и 10 переменных (все физиологические характеристики), поэтому они очень высокие и узкие. Зависимой переменной является количественное измерение прогрессирования заболевания после базового уровня за один год.

Это классический набор данных, который использовался в статье Эфрона, Хасти, Джонстона и Тибширани о регрессии по методу наименьших углов, а также включён в множество наборов данных scikit-learn.

data = datasets.load_diabetes()

feature_names=['age', 'sex', 'bmi', 'bp', 's1', 's2', 's3', 's4', 's5', 's6']

trn,test,y_trn,y_test = train_test_split(data.data, data.target, test_size=0.2)

trn.shape, test.shape

# ((353, 10), (89, 10))

Линейная регрессия в Sklearn

Рассмотрим систему Xβ=y, где строки X больше столбцов. Это происходит, когда у вас больше данных образцов, чем переменных. Мы хотим найти , чтобы минимизировать:

Давайте начнём с реализации в sklearn:

regr = linear_model.LinearRegression()
%timeit regr.fit(trn, y_trn)

# 458 µs ± 62.4 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

pred = regr.predict(test)

Есть несколько показателей, которые могут помочь нам оценить качество наших прогнозов. Мы рассмотрим среднеквадратичную ошибку (L2) и среднюю абсолютную ошибку (L1).

def regr_metrics(act, pred):
    return (math.sqrt(metrics.mean_squared_error(act, pred)), 
    metrics.mean_absolute_error(act, pred))

regr_metrics(y_test, regr.predict(test))

# (75.36166834955054, 60.629082113104403)

Полиномиальные признаки

Линейная регрессия находит оптимальные коэффициенты βi:

Добавление полиномиальных признаков всё ещё является линейной регрессией, просто требуется больше членов:

Нам нужно использовать исходные данные X для вычисления других полиномиальных признаков.

trn.shape

# (353, 10)

Теперь мы хотим попробовать улучшить производительность модели, добавив больше функций. В настоящее время наша модель линейна для каждой переменной, но мы можем добавить полиномиальные функции, чтобы изменить её.

poly = PolynomialFeatures(include_bias=False)

trn_feat = poly.fit_transform(trn)

', '.join(poly.get_feature_names(feature_names))

# 'age, sex, bmi, bp, s1, s2, s3, s4, s5, s6, age^2, age sex, age bmi, age bp, age s1, age s2, age s3, age s4, age s5, age s6, sex^2, sex bmi, sex bp, sex s1, sex s2, sex s3, sex s4, sex s5, sex s6, bmi^2, bmi bp, bmi s1, bmi s2, bmi s3, bmi s4, bmi s5, bmi s6, bp^2, bp s1, bp s2, bp s3, bp s4, bp s5, bp s6, s1^2, s1 s2, s1 s3, s1 s4, s1 s5, s1 s6, s2^2, s2 s3, s2 s4, s2 s5, s2 s6, s3^2, s3 s4, s3 s5, s3 s6, s4^2, s4 s5, s4 s6, s5^2, s5 s6, s6^2'

trn_feat.shape

# (353, 65)

regr.fit(trn_feat, y_trn)

# LinearRegression(copy_X=True, fit_intercept=True, n_jobs=1, normalize=False)

regr_metrics(y_test, regr.predict(poly.fit_transform(test)))

# (55.747345922929185, 42.836164292252235)

Время квадратично зависит от количества признаков и линейно от количества образцов, так что это будет очень медленно!

%timeit poly.fit_transform(trn)

# 635 µs ± 9.25 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

Ускорение генерации признаков

Мы хотим ускорить процесс. Мы будем использовать Numba, библиотеку Python, которая компилирует код непосредственно в C.

Numba — это компилятор.

Ресурсы

Учебник Джейка Вандер Пласа — хороший обзор. Здесь Джейк использует Numba для реализации нетривиального алгоритма (неравномерное быстрое преобразование Фурье).

Cython — ещё один вариант. Я обнаружил, что Cython в основном требует больше знаний, чем Numba (он ближе к C), но обеспечивает аналогичное ускорение.

Здесь представлен компилятор AOT (Ahead-of-Time), JIT (Just-In-Time) компилятор и интерпретатор.

Ответ на вопрос о различиях между традиционным интерпретатором, JIT-компилятором, JIT-интерпретатором.

Эксперименты с векторизацией и нативным кодом

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

%matplotlib inline

import math, numpy as np, matplotlib.pyplot as plt
from pandas_summary import DataFrameSummary
from scipy import ndimage

from numba import jit, vectorize, guvectorize, cuda, float32, void, float64

Мы продемонстрируем следующие эффекты:

• Избегание выделения памяти и копирования (медленнее, чем вычисления на CPU)
• Улучшенная локальность
• Векторизация

Если мы используем numpy для всего массива сразу, он создаст много временных значений и не сможет использовать кэш. Если мы перебираем элементы массива с помощью цикла Numba, нам не нужно выделять большой временный массив, и мы можем повторно использовать данные кэша, потому что мы выполняем многократные вычисления для каждого элемента массива. proc_python(x,y)

49.8 ms ± 1.19 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

### NumPy

Numpy позволяет нам работать с векторизацией:

```py
# С типами и векторизацией
def proc_numpy(x, y):
    z = np.zeros(nobs, dtype='float32')
    x = x * 2 - (y * 55)
    y = x + y * 2
    z = x + y + 99
    z = z * (z - .88)
    return z

np.allclose(proc_numpy(x, y), proc_python(x, y), atol=1e-4)

True

%timeit proc_numpy(x, y)  # Типизированная и векторизованная

35.9 µs ± 166 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

Numba

Numba предоставляет несколько различных декораторов. Мы попробуем два разных метода:

• @jit — очень общий;
• @vectorize — не нужно писать циклы for. Полезно при работе с векторами одинакового размера.

Сначала мы будем использовать декоратор компилятора Numba @jit, без явной векторизации. Это позволяет избежать большого объёма выделения памяти, поэтому у нас есть лучшая локальность:

@jit()
def proc_numba(xx, yy, zz):
    for j in range(nobs):
        x, y = xx[j], yy[j]
        x = x * 2 - (y * 55)
        y = x + y * 2
        z = x + y + 99
        z = z * (z - .88)
        zz[j] = z
    return zz

z = np.zeros(nobs).astype('float32')
np.allclose(proc_numpy(x, y), proc_numba(x, y, z), atol=1e-4)

True

%timeit proc_numba(x, y, z)

6.4 µs ± 17.6 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

Теперь мы будем использовать декоратор @vectorize. Компилятор Numba оптимизирует его более умным способом по сравнению с обычным Python и Numpy. Он создаёт для вас Numpy ufunc, что традиционно включает написание на C и не так просто.

@vectorize
def vec_numba(x, y):
    x = x * 2 - (y * 55)
    y = x + y * 2
    z = x + y + 99
    return z * (z - .88)

np.allclose(vec_numba(x, y), proc_numba(x, y, z), atol=1e-4)

True

%timeit vec_numba(x, y)

5.82 µs ± 14.4 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

Numba великолепна. Посмотрите, как быстро это работает!

Многочлены Numba

@jit(nopython=True)
def vec_poly(x, res):
    m, n = x.shape
    feat_idx = 0
    for i in range(n):
        v1 = x[:, i]
        for k in range(m): res[k, feat_idx] = v1[k]
        feat_idx += 1
        for j in range(i, n):
            for k in range(m): res[k, feat_index] = v1[k] * x[k, j]
            feat_index += 1

Строчный и столбцовый порядок хранения

Из статьи Эли Бендерски:

«Матрицы в строчном порядке располагаются так, что первая строка находится в непрерывной памяти, затем вторая строка и так далее. Матрицы в столбцовом порядке располагаются таким образом, что первый столбец находится в непрерывной памяти, за ним следует второй столбец и так далее».

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

«Оказывается, способ работы алгоритма с макетом данных может определять производительность приложения».

«Короче говоря, всегда проходите через данные в порядке макета».

Столбцовый порядок: Fortran, Matlab, R и Julia.

Строчный порядок: C, C++, Python, Pascal, Mathematica.

trn = np.asfortranarray(trn)
test = np.asfortranarray(test)

m, n = trn.shape
n_feat = n * (n + 1) // 2 + n
trn_feat = np.zeros((m, n_feat), order='F')
test_feat = np.zeros((len(y_test), n_feat), order='F')

vec_poly(trn, trn_feat)
vec_poly(test, test_feat)

regr.fit(trn_feat, y_trn)

LinearRegression(copy_X=True, fit_intercept=True, n_jobs=1, normalize=False)

regr_metrics(y_test, regr.predict(test_feat))

(55.74734592292935, 42.836164292252306)

%timeit vec_poly(trn, trn_feat)

7.33 µs ± 19.8 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

Вспомните, что это реализация sklearn PolynomialFeatures, созданная экспертами:

%timeit poly.fit_transform(trn)

635 µs ± 9.25 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

605/7.7

78.57142857142857

Это серьёзная проблема! Numba великолепна! Всего лишь одна строка кода, и мы получаем скорость, которая в 78 раз быстрее, чем у scikit-learn (созданного экспертами).

Регуляризация и шум

Регуляризация — это метод уменьшения переобучения и создания модели, которая лучше обобщается на новые данные.

Регуляризация

Lasso регрессия использует L1 штраф, который приводит к разреженным коэффициентам. Параметр α используется для взвешивания штрафной функции. Scikit Learn LassoCV использует множество различных значений α для перекрестной проверки.

Посмотрите видео Coursera о Lasso регрессии, чтобы узнать больше. ### Шум

Теперь мы добавим в данные некоторый шум.

idxs = np.random.randint(0, len(trn), 10)

y_trn2 = np.copy(y_trn)
y_trn2[idxs] *= 10 # label noise

regr = linear_model.LinearRegression()
regr.fit(trn, y_trn)
regr_metrics(y_test, regr.predict(test))

# (51.1766253181518, 41.415992803872754)

regr.fit(trn, y_trn2)
regr_metrics(y_test, regr.predict(test))

# (62.66110319520415, 53.21914420254862)

Huber-потери — это функция потерь, которая менее чувствительна к выбросам по сравнению с квадратичной функцией потерь. Для малых ошибок она квадратична, а для больших — линейна.

$L(x)= \begin{cases} \frac{1}{2}x^2, & \text{for } \lvert x\rvert\leq \delta \\ \delta(\lvert x \rvert - \frac{1}{2}\delta), & \text{otherwise} \end{cases}$

hregr = linear_model.HuberRegressor()
hregr.fit(trn, y_trn2)
regr_metrics(y_test, hregr.predict(test))

# (51.24055602541746, 41.670840571376822)