План исследования
Более детальное описание тут.
import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
import seaborn as sns
import scipy
from statsmodels.stats.weightstats import *
from sklearn.linear_model import RidgeCV, Ridge, Lasso, LassoCV
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import train_test_split, KFold
from sklearn.preprocessing import StandardScaler, LabelBinarizer, PolynomialFeatures
from sklearn.metrics import mean_absolute_error, make_scorer
from xgboost import XGBClassifier
from hyperopt import fmin,tpe, hp, STATUS_OK, Trials
import xgboost as xgb
from sklearn.model_selection import learning_curve, validation_curve
%matplotlib inline
Датасет содержит информацию о 53940 бриллиантах. По некоторым характеристикам (об этом позже) будем предсказывать стоимость. Данные можно скачать здесь.
С точки зрения бизнеса, ценность задачи понятна - по характеристикам бриллианта предсказать, сколько долларов за него можно получить. От бизнеса я далёк, поэтому интерес чисто спортивный: разобраться, какие характеристики и как влияют на стоимость этих камешков =)
Признаки
Целевой признак: price - стоимость бриллианта в долларах
#загружаем dataset
diamonds_df = pd.read_csv('../../data/diamonds.csv')
diamonds_df.head()
diamonds_df.describe()
Видно, что масштабы признаков отличаются. В дальнейшем нужно будет применить StandartScaler
diamonds_df.info()
В данных отсутствуют пропуски. Итого, имеется 6 вещественных, 1 целочисленный (unnamed: 0 не считаем) и 3 категориальных признака.
real_features = ['carat', 'depth', 'table', 'x', 'y', 'z','price']
# Изучим корреляцию вещественных признаков и целевой переменной
sns.heatmap(diamonds_df[real_features].corr(method='spearman'));
Признаки carat, x,y,z имеют большую корреляцию, как между собой, так и с целевой переменной, что не удивительно. При этом, корреляция целевой переменной и признаков depth, table почти отсутствует
cat_features = ['cut','color','clarity']
fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(16, 10))
for idx, feature in enumerate(cat_features):
sns.countplot(diamonds_df[feature], ax=axes[idx % 3], label=feature)
Реальные значения категориальных признаков не отличаются от тех, что заявлены в описании. Кроме того, видно, что уникальных значений не много, так что One Hot encoding должен отлично сработать.
sns.distplot(diamonds_df['price'])
Распределение имеет тяжелый правый хвост. Применим логарифмирование.
sns.distplot(diamonds_df['price'].map(np.log1p))
Помогло это не сильно: получилось бимодальное распределение. Зато хвост исчез =) Для наглядности, построим QQ график
stats.probplot(diamonds_df['price'], dist="norm", plot=plt);
# Начнем с построения гистограмм вещественных признаков
fig, axes = plt.subplots(nrows=2, ncols=3, figsize=(16, 10))
for idx, feature in enumerate(real_features[:-1]): #price рисовать не будем
sns.distplot(diamonds_df[feature], ax=axes[idx // 3, idx % 3], label=feature)
Распределение признаков depth, table, y, z отдаленно, но напоминает колокол. У depth хвосты тяжеловаты для нормального распределения; carat и table скорее бимодальные. Кроме того, у них тяжелые правые хвосты, так что np.log1p не помешает. По графикам выше не видно выбросов. Проверим, что это действительно так, с помощью boxplot
fig, axes = plt.subplots(nrows=2, ncols=3, figsize=(16, 10))
for idx, feature in enumerate(real_features[:-1]): #price рисовать не будем
sns.boxplot(diamonds_df[feature], ax=axes[idx // 3, idx % 3], orient='v')
Каких-либо серьезных аномалий в рассматриваемых данных нет. На всякий случай посмотрим бриллиант с y=60, z = 32 и carat > 4. Если он стоит дорого, то за выброс его считать не будем.
diamonds_df[diamonds_df['y'] > 55].head()
diamonds_df[diamonds_df['z'] > 30].head()
diamonds_df[diamonds_df['carat'] > 4].head()
Видно, что это просто очень дорогие камни. Посмотрим, как рассматриваемые признаки взаимосвязаны с целевой перменной
sns.pairplot(diamonds_df[real_features], diag_kind="kde")
Посмотрим, как целевая переменная зависит от категориальных признаков
# цвет бриллианта
fig, axes = plt.subplots(nrows=3, ncols=3, figsize=(16, 10))
for idx, (color, sub_df) in enumerate(pd.groupby(diamonds_df, 'color')):
ax = sns.distplot(sub_df['price'], kde=False, ax=axes[idx // 3, idx % 3])
ax.set(xlabel=color)
Распределения для всех значений цветов имеют тяжелый правый хвост и не сильно отличаются друг от друга.
# чистота бриллианта
fig, axes = plt.subplots(nrows=3, ncols=3, figsize=(16, 10))
for idx, (clarity, sub_df) in enumerate(pd.groupby(diamonds_df, 'clarity')):
ax = sns.distplot(sub_df['price'], kde=False, ax=axes[idx // 3, idx % 3])
ax.set(xlabel=clarity)
Хвосты у всех тяжелые, но у SI1,SI2 присутствуют дополнительные пики в районе 5000.
# качество огранки
fig, axes = plt.subplots(nrows=2, ncols=3, figsize=(16, 10))
for idx, (cut, sub_df) in enumerate(pd.groupby(diamonds_df, 'cut')):
ax = sns.distplot(sub_df['price'], kde=False, ax=axes[idx // 3, idx % 3])
ax.set(xlabel=cut)
И снова пики в районе 5000 (у Good и Premium). А в целом графики похожи.
Нарисуем boxplot для каждого значения
fig, axes = plt.subplots(nrows=3, ncols=1, figsize=(16, 10))
# Отобразим строки в числа в порядке от худшего к лучшему. Так удобнее на графике смотреть
df = diamonds_df.copy()
df['color'] = df['color'].map({'J': 0, 'I': 1, 'H': 2, 'G': 3, 'F': 4, 'E': 5, 'D': 6})
df['clarity'] = df['clarity'].map({'I1': 0, 'SI2': 1, 'SI1': 2, 'VS2': 3, 'VS1': 4, 'VVS2': 5, 'VVS1': 6, 'IF': 7 })
df['cut'] = df['cut'].map({'Fair': 0, 'Good': 1, 'Very Good': 2, 'Premium': 3, 'Ideal': 4})
for idx, feature in enumerate(cat_features):
sns.boxplot(x=feature, y='price',data=df,hue=feature, ax=axes[idx])
Тут уже интереснее. Начнем с огранки. Видно, что медиана максимальна для Very Good и Premium. Для ideal медианное значение цены гораздо меньше. Аналогичные наблюдения можно сделать для цвета и чистоты. Возможно, бриллианты с наилучшими свойствами на очень большие, и, соответсвенно, их цена ниже. Проверим это.
fig, axes = plt.subplots(nrows=3, ncols=1, figsize=(16, 10))
for idx, feature in enumerate(cat_features):
sns.boxplot(x=feature, y='carat',data=df,hue=feature, ax=axes[idx])
Действительно, медианное значение веса для бриллиантов с очень хорошими характеристиками меньше, чем для бриллиантов с плохими харакетристиками. Напоследок, посмотрим сколько бриллиантов с той или иной харакеристикой присутствует в данных.
fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(16, 10))
for idx, feature in enumerate(cat_features):
sns.countplot(df[feature], ax=axes[idx % 3], label=feature)
Видно, что очень мало камней с плохой огранкой. Также мало камней с плохими цветовыми харакеристиками. Но и не очень много с идеальными. Распределение чистоты камня напоминает лапласовское распределение.
# Для начала, выделим выборку для тестирования
X = diamonds_df.drop(['price'], axis=1).values[:,1:] # отсекаем индекс
y = diamonds_df['price']
X_train, X_test, y_train, y_test = train_test_split(X,y,test_size=0.3,random_state=4444, shuffle=True)
# признаки с индексами 1,2,3 категориальные. Применим к ним ohe
label_bin = LabelBinarizer()
X_train_cut_ohe = label_bin.fit_transform(X_train[:,1])
X_test_cut_ohe = label_bin.transform(X_test[:,1])
X_train_color_ohe = label_bin.fit_transform(X_train[:,2])
X_test_color_ohe = label_bin.transform(X_test[:,2])
X_train_clarity_ohe = label_bin.fit_transform(X_train[:,3])
X_test_clarity_ohe = label_bin.transform(X_test[:,3])
# carat, x и целевую переменную логарифмируем
log_vect = np.vectorize(np.log1p)
X_train_сarat_log = log_vect(X_train[:,0]).reshape(-1,1)
X_test_сarat_log = log_vect(X_test[:,0]).reshape(-1,1)
X_train_x_log = log_vect(X_train[:,6]).reshape(-1,1)
X_test_x_log = log_vect(X_test[:,6]).reshape(-1,1)
y_train_log = log_vect(y_train)
y_test_log = log_vect(y_test)
# масштабириуем вещественные признаки
scaler = StandardScaler()
X_train_real = np.hstack((X_train_сarat_log, X_train_x_log, X_train[:,[7,8,4,5]]))
X_test_real = np.hstack((X_test_сarat_log, X_test_x_log, X_test[:,[7,8,4,5]]))
X_train_real_scaled = scaler.fit_transform(X_train_real)
X_test_real_scaled = scaler.transform(X_test_real)
# В качестве дополнительных признаков будем рассматривать полиномиальные признаки
#Данные признаки должны улучшить качество линейной модели.
X_train_additional = PolynomialFeatures().fit_transform(X_train_real)
X_test_additional = PolynomialFeatures().fit_transform(X_test_real)
X_train_additional_scaled = scaler.fit_transform(X_train_additional)
X_test_additional_scaled = scaler.transform(X_test_additional)
# Объединяем все преобразованные признаки
X_train_transformed = np.hstack((X_train_real_scaled,X_train_cut_ohe, X_train_color_ohe, X_train_clarity_ohe))
X_test_transformed = np.hstack((X_test_real_scaled,X_test_cut_ohe, X_test_color_ohe, X_test_clarity_ohe))
Смотри предыдущий пункт
Рассмотрим сначала линейную модель. Данные разделим на 5 фолдов. С помощью RidgeCV и LassoCV будем оптимизировать силу регуляризации.
# функция потерь для рассматриваемой задачи. Ошибку смотрим на исходных данных
def mean_absolute_exp_error(model, X,y):
return -mean_absolute_error(np.expm1(model.predict(X)), np.expm1(y))
cv = KFold(n_splits=5, shuffle=True, random_state=4444)
alphas = np.logspace(-5,2,100)
ridge_cv = RidgeCV(alphas=alphas, scoring=mean_absolute_exp_error, cv=cv)
lasso_cv = LassoCV(alphas=alphas, cv=cv, random_state=4444)
ridge_cv.fit(X_train_transformed, y_train_log)
lasso_cv.fit(X_train_transformed, y_train_log)
print('Optimized alpha: Ridge = %f, Lasso = %f' % (ridge_cv.alpha_, lasso_cv.alpha_))
score_ridge = mean_absolute_error(y_test, np.expm1(ridge_cv.predict(X_test_transformed)))
score_lasso = mean_absolute_error(y_test, np.expm1(lasso_cv.predict(X_test_transformed)))
print('Ridge regression score = %f' % score_ridge)
print('Lasso regression score = %f' % score_lasso)
Оба метода показали схожий результат. Что будет, если мы добавим новые признаки?
X_train_transformed_add = np.hstack((X_train_transformed, X_train_additional_scaled))
X_test_transformed_add = np.hstack((X_test_transformed, X_test_additional_scaled))
ridge_cv.fit(X_train_transformed_add, y_train_log)
lasso_cv.fit(X_train_transformed_add, y_train_log)
print('Optimized alpha: Ridge = %f, Lasso = %f' % (ridge_cv.alpha_, lasso_cv.alpha_))
score_ridge = mean_absolute_error(y_test, np.expm1(ridge_cv.predict(X_test_transformed_add)))
score_lasso = mean_absolute_error(y_test, np.expm1(lasso_cv.predict(X_test_transformed_add)))
print('Ridge regression score = %f' % score_ridge)
print('Lasso regression score = %f' % score_lasso)
Ошибка значительно уменьшилась. Построим кривые валидации и обучения
%%time
# код из статьи на хабре
def plot_with_err(x, data, **kwargs):
mu, std = data.mean(1), data.std(1)
lines = plt.plot(x, mu, '-', **kwargs)
plt.fill_between(x, mu - std, mu + std, edgecolor='none',
facecolor=lines[0].get_color(), alpha=0.2)
model = Ridge(random_state=4444)
alphas = np.logspace(1,2,10) + 10 # Если коэффициент регуляризации мал, то значения получаются заоблочными
val_train, val_test = validation_curve(model, X_train_transformed_add, y_train_log,'alpha', alphas, cv=cv,scoring=mean_absolute_exp_error)
plot_with_err(alphas, -val_train, label='training scores')
plot_with_err(alphas, -val_test, label='validation scores')
plt.xlabel(r'$\alpha$'); plt.ylabel('MAE')
plt.legend();
Судя по кривым валидации, модель недообучилась: ошибки лежат близко друг к другу.
# код из статьи на хабре
def plot_learning_curve(model, X,y):
train_sizes = np.linspace(0.05, 1, 20)
N_train, val_train, val_test = learning_curve(model,X, y, train_sizes=train_sizes, cv=5,scoring=mean_absolute_exp_error, random_state=4444)
plot_with_err(N_train, -val_train, label='training scores')
plot_with_err(N_train, -val_test, label='validation scores')
plt.xlabel('Training Set Size'); plt.ylabel('MAE')
plt.legend()
model = Ridge(alpha=52.140083,random_state=4444)
plot_learning_curve(model, X_train_transformed_add, y_train_log)
Кривые лежат близко друг к другу почти с самого начала. Вывод: наблюдений у нас достаточно, нужно двигаться в сторону усложнения модели
Случайный лес должен хорошо работать "из коробки". Поэтому будем оптимизировать только число деревьев.
%%time
model = RandomForestRegressor(n_estimators=100, random_state=4444)
n_estimators = [10,25,50,100,250,500,1000]
val_train, val_test = validation_curve(model, X_train_transformed, y_train_log,'n_estimators', n_estimators, cv=cv,scoring=mean_absolute_exp_error)
plot_with_err(n_estimators, -val_train, label='training scores')
plot_with_err(n_estimators, -val_test, label='validation scores')
plt.xlabel('n_estimators'); plt.ylabel('MAE')
plt.legend();
Видно, что начиная с 200 деревьев качество практически не изменяется. Поэтому в качестве еще одной модели будем рассматривать случайный лес именно с таким количеством деревьев.
forest_model = RandomForestRegressor(n_estimators=200, random_state=4444)
forest_model.fit(X_train_transformed, y_train_log)
forest_prediction = np.expm1(forest_model.predict(X_test_transformed))
score = mean_absolute_error(y_test, forest_prediction)
print('Random forest score: %f' % score)
# посмотрим на важность признаков
np.argsort(forest_model.feature_importances_)
Первые четыре столбца обучающей выборки соответствуют признакам carat, x,y,z. Как и предполагалось в начале, 3 из 4 этих признаков имеют наибольшую важность для модели
%%time
# Построим, также, кривую обучения
plot_learning_curve(model, X_train_transformed, y_train_log)
График выходит на полку, так что больше данных нам не нужно
X_train_boosting, X_valid_boosting, y_train_boosting, y_valid_boosting = train_test_split(
X_train_transformed, y_train_log, test_size=0.3, random_state=4444)
def score(params):
from sklearn.metrics import log_loss
print("Training with params:")
print(params)
params['max_depth'] = int(params['max_depth'])
dtrain = xgb.DMatrix(X_train_boosting, label=y_train_boosting)
dvalid = xgb.DMatrix(X_valid_boosting, label=y_valid_boosting)
model = xgb.train(params, dtrain, params['num_round'])
predictions = model.predict(dvalid).reshape((X_valid_boosting.shape[0], 1))
score = mean_absolute_error(np.expm1(y_valid_boosting), np.expm1(predictions))
# score = mean_absolute_error(y_valid_boosting, predictions)
print("\tScore {0}\n\n".format(score))
return {'loss': score, 'status': STATUS_OK}
def optimize(trials):
space = {
'num_round': 200,
'learning_rate': hp.quniform('eta', 0.05, 0.5, 0.005),
'max_depth': hp.quniform('max_depth', 3, 14, 1),
'min_child_weight': hp.quniform('min_child_weight', 1, 10, 1),
'subsample': hp.quniform('subsample', 0.5, 1, 0.05),
'gamma': hp.quniform('gamma', 0.5, 1, 0.01),
'colsample_bytree': hp.quniform('colsample_bytree', 0.4, 1, 0.05),
'eval_metric': 'mae',
'objective': 'reg:linear',
'nthread' : 4,
'silent' : 1,
'seed': 4444
}
best = fmin(score, space, algo=tpe.suggest,trials=trials, max_evals=100)
return best
%%time
# Оптимизация параметров
trials = Trials()
best_params = optimize(trials)
best_params
params = {
'num_round': 200,
'colsample_bytree': 0.65,
'eta': 0.145,
'gamma': 0.55,
'max_depth': 10,
'min_child_weight': 4.0,
'subsample': 1.0,
'eval_metric': 'mae',
'objective': 'reg:linear',
'nthread' : 4,
'silent' : 1,
'seed': 4444}
dtrain = xgb.DMatrix(X_train_transformed, label=y_train_log)
dvalid = xgb.DMatrix(X_test_transformed, label=y_test_log)
boosting_model = xgb.train(params, dtrain, params['num_round'])
predictions = boosting_model.predict(dvalid).reshape((X_test_transformed.shape[0], 1))
score = mean_absolute_error(y_test, np.expm1(predictions))
print('Boosting score: %f' % score)
В большом количестве в предыдущем пункте
В большом количестве в части 7
Приведем результаты различных моделей на тестовой выборке. Как уже оговаривалось ранее, в качестве метрики используем MAE
pure_ridge = Ridge(random_state=4444, alpha=0.00001) # гребневая регрессия на исходных данных
pure_ridge.fit(X_train_transformed, y_train_log)
pure_ridge_score = mean_absolute_error(y_test, np.expm1(pure_ridge.predict(X_test_transformed)))
print('Ridge regression score: %f' % pure_ridge_score)
poly_ridge = Ridge(random_state=4444, alpha=52.140083) # гребневая регрессия с полиномиальными признаками
poly_ridge.fit(X_train_transformed_add, y_train_log)
poly_ridge_score = mean_absolute_error(y_test, np.expm1(poly_ridge.predict(X_test_transformed_add)))
print('Ridge regression score with poly features: %f' % poly_ridge_score)
forest_score = mean_absolute_error(y_test, np.expm1(forest_model.predict(X_test_transformed)))
print('Random forest score: %f' % forest_score)
boosting_score = mean_absolute_error(y_test, np.expm1(boosting_model.predict(dvalid)))
print('XGBoost score: %f' % boosting_score)
Результаты близки к тем, что получались на кросс-валидации. Так что всё хорошо =)
В данном проекте рассматривались достаточно "простые" данные, поэтому основной упор был сделан на применение различных моделей для их анализа. С одной стороны, случайный лес без какой-либо настройки гиперпараметров показал лучший результат. С другой стороны, если потратить больше времени на оптимизацию градиентного бустинга, возможно, он сможет показать результат лучше, чем у случайного леса. Стоит, также, отметить линейную модель: после добавления полиномиальных признаков она показала очень неплохой результат (если сравнивать с моделью без дополнительных признаков =)). Зато сложность гораздо меньше. Если вдруг кому-то по жизни придется оценивать бриллианты, можете смело использовать предложенную модель случайного леса. В среднем будете терять по 275 $ с одного камушка :p
Спасибо за внимание!