Этот набор данных содержит цены продажи жилья для округа Кинг, с административным центром в г. Сиэтл за период с мая 2014 года по май 2015 года.
Краткая информация из Википедии об округе Кинг (англ. King County)
Крупные города: Сиэтл (окружной центр и самый большой город), Белвью, Такома
Население округа: 2 188 649 (на 2017 год)
Территория:
Сиэтл (англ. Seattle) — крупнейший город на северо-западе США и в штате Вашингтон, крупный морской порт. Расположен между системой заливов Пьюджет и озером Вашингтон.
Адресс нахождения данных: https://www.kaggle.com/harlfoxem/housesalesprediction/data
id - Номер обеъкта
date - Дата продажи объекта неджвижимости
price- Цена продажи
bedrooms - Кол-во спален
bathrooms - Кол-во ванных+сан.узлов
sqft_living- Жилая площадь дома
sqft_lot - Общая площадь участка
'floors' - Кол-во этажей в доме
"waterfront" - Вид на водоем (Да/Нет)
view - Был просмотрен (не очень понял что тут имелось ввиду)
condition - Насколько хорошо состояние (в целом)
grade - Общий класс(грейд), присвоенный обеъкту неджвижимости, основанный на системе классификации Кинг Кантри
sqft_above - Площадь объекта недвижимости, не включая подвал(если есть)
sqft_basement - Площадь подвала (если есть)
yr_built - Год постройки
yr_renovated - Год реновации (если была)
zipcode - Почтовый индекс
lat - Широта
long - Долгота
sqft_living15 - Жилая площадь в 2015 году (подразумевает некоторые ремонтные работы/ перепланировку). Это может повлиять или не повлиять на общую площадь территории
sqft_lot15 - Площадь общей территории в 2015 году(подразумевает некоторые ремонтные работы/ перепланировку)
https://moto.data.socrata.com/dataset/King-County-Sheriff-s-Office/4h35-4mtu - Информация о преступлениях в округе Кинг
https://data.kingcounty.gov - Информационный портал округа Кинг(несколько разных отчетов по тематикам, Школы,Клиники,Магазины и общепит из которых я вытаскивал местонаходжение объектов и сами объекты)
Данные в агрегированном виде и основной датасет: https://cloud.mail.ru/public/DTDE/zax4jHwSL
Исследовать данные по продажам недвижимости округа Кинг.
Построить модель наиболее точно предсказывающую стоимость объектов недвижимости на основании имеющихся данных и за счет сбора дополнительных.
# Импорт всех нужных библиотек
import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
%matplotlib inline
import seaborn as sns
from sklearn.model_selection import train_test_split, GridSearchCV, cross_val_predict, cross_val_score
from sklearn.base import BaseEstimator, ClassifierMixin
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score
from sklearn.linear_model import Lasso, Ridge
from sklearn.neighbors import KNeighborsRegressor
from sklearn.pipeline import Pipeline
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.cluster import KMeans
import xgboost as xgb
import lightgbm as lgb
from scipy import stats
from scipy.special import inv_boxcox, boxcox1p
from scipy.stats.mstats import gmean
from scipy.stats import norm, skew , boxcox
from scipy.sparse import csr_matrix
pd.set_option('display.max_columns', 100)
import warnings
warnings.filterwarnings('ignore')
data = pd.read_csv('kc_house_data.csv', parse_dates =['date'])
data.head()
data.info()
Данные не имеют пропусков, что уже хорошо.
data.describe()
ind = len(data)*0.8
work_data = data.ix[:ind]
validation_data = data.ix[ind:]
data.shape, work_data.shape, validation_data.shape
Перед анализом данных, поскольку сам дата сет без пропусков и ошибок, я построю простую модель RandomForest и использую его метрику качества score - по умолчанию R2 - коэф. детерминации, которую он считает на отложенной выборке. Дополнительно буду смотреть на среднеквадратичную логарифмическую ошибку (достаточно частая метрика в пердсказании стоимости недвижимости). Также посмотрю как они (R2 и RSMLE) зависимы друг от друга практически. RandomForest также выбираю из практических соображений, чтобы сразу посмотреть на значимость признаков
#Функция среднеквадратичной логарифмической ошибки
def rmsle(y_true,y_pred):
assert len(y_true) == len(y_pred)
return np.square(np.log(y_pred + 1) - np.log(y_true + 1)).mean() ** 0.5
df_train = work_data.drop([ 'id','price', 'date'], axis=1)
target = work_data['price']
X_train, X_test, y_train, y_test = train_test_split(df_train, target, test_size=0.25, random_state=17)
clf = RandomForestRegressor(random_state=17)
clf.fit(X_train.values, y_train)
y_pred = clf.predict(X_test.values)
clf.score(X_test.values, y_test), rmsle(y_test, y_pred)
#Сразу посмотрим какие признаки алгоритм оранжировал по степени важности
feature_importance = pd.DataFrame()
feature_importance['feature'] = df_train.columns
feature_importance['importances'] = clf.feature_importances_
feature_importance.sort_values('importances', ascending=False).head(20)
Посмотрим на график зависимости цены от первых 5 самых важных признаков из базовго алгоритма, а также их зависимость между собой
sns.set()
cols = feature_importance['feature'][feature_importance['importances']> 0.031]
sns.pairplot(work_data[cols], size = 2.5)
plt.show();
Из графиков видно, что на высокий грейд линейно влияет общая жилая площадь, а также немного жилая площадь по оценке на 2015 год (учитывающая реновации).
# Посмотрим как распределены продажи по zipcode
sort_data = work_data.groupby('zipcode', sort = 'id').count().reset_index()
plt.figure(figsize=(15,10))
ax= sns.barplot(x=sort_data['zipcode'], y=sort_data['bedrooms'])
plt.xticks(rotation= 90)
plt.xlabel('Zipcode')
plt.ylabel('Кол-во объектов')
plt.title("Распределение кол-ва проданных объектов относительно их районов, выделенных по zipcode")
В некоторых районах недвижимости продается в разы больше чем в других. Невысокиие продажи могут говорить как правило о двух противоположных тенеденциях. Дорогой район и там состав жильцов не меняется. Бедный район и там мало кто покупает. Посмотрим на других срезах, подтвердится ли гипотеза.
# Посмотрим как распределена средняя стоимость проданных объектов относительно их районов, выделенных по zipcode
sort_data_pz = work_data.groupby('zipcode', sort = 'id').mean().reset_index()
plt.figure(figsize=(15,10))
ax= sns.barplot(x=sort_data_pz['zipcode'], y=sort_data_pz['price'])
plt.xticks(rotation= 90)
plt.xlabel('Zipcode')
plt.ylabel('Средняя стоимость')
plt.title("Распределение средней стоимости проданных объектов относительно их районов, выделенных по zipcode")
Видно что есть несколько индексов(4), среденяя стоимость жилья возле которых превышает 1 млн. долларов. По всей видимости это какие-то элитные районы. Можно создать дополнительно 3 признака класса жилья по стоимости. Частичное подтверждение гипотезы по индексу 98039 - продается мало объектом и средняя стоимостьсвыше 2млн. долларов. - богатый район.
98039 - Медина, Вашингтон, США - входит в топ 10 самых дорогих районов в США, ср. стоимость недвижимости (по данным из интернета за 2015 год) около 2 998 000 долларов, в нашем графике около 2 100 000, но у нас все что реально продалось, а аналитические журналы часто аппелируют средней заявленной стартовой ценой продавца.
sort_data_yr = work_data.groupby('yr_built').mean().reset_index()
plt.figure(figsize=(15,10))
ax= sns.barplot(x=sort_data_yr['yr_built'], y=sort_data_yr['price'])
plt.xticks(rotation= 90)
plt.xlabel('Год постройки')
plt.ylabel('Средняя стоимость')
plt.title("Распределение средней стоимости относительно года постройки")
Интересно, что на графике видно, что средняя стоимость постоек до 1941 года и после 1994 выше, чем в период между этими датами. Это может свидетельствовать о таких факторах, как:
Т.е. возможно, постройки до 41 и после 94 - это строительство примерно в одних и тех же районах. Нужно построить график распределения по годам построек по индексам.
plt.figure(figsize=(15,10))
sort_data_zy = work_data.groupby(['zipcode','yr_built']).count().reset_index()
sns.boxplot(x="zipcode", y="yr_built", data=sort_data_zy)
plt.xticks(rotation= 90)
plt.xlabel('Индекс')
plt.ylabel('Год постройки')
plt.title("Распределение возраста объектов недвижимости в засисимости от индекса")
Очевидно, что есть более "молодые районы" и более "старинные", также можно предположить, что индексы с "широким ящиком" это районы близкие к историческому центру городов или они и являются центрами, также из теории урбанизации можно сделать вывод что такие районы ближе к воде (к историческому порту), если город находится у судоходного водоема.
Проверил индекс 98102 - это почти центр г. Сиэтл и близкий к воде - т.е. предположение подтвердилось
#Выделим из даты месяц продажи и посмотрим на распределение продаж по месяцам
work_data['month_sold'] = work_data['date'].apply(lambda ts: ts.month)
validation_data['month_sold'] = validation_data['date'].apply(lambda ts: ts.month)
plt.figure(figsize=(10,7))
sns.countplot(x=work_data.month_sold)
plt.xlabel('Месяц продажи')
plt.ylabel('Кол-во объектов')
plt.title('График кол-ва продаж по месяцам',color = 'blue',fontsize=15)
Видна определенная сезоннасть весной продают больше, чем зимой и т.д.
# Посмотрим на расперделение по грейду
plt.figure(figsize=(10,7))
sns.countplot(x=work_data.grade)
plt.xlabel('Грейд')
plt.ylabel('Кол-во объектов')
plt.title('График распределения грейда среди всех проданных объектов',color = 'blue',fontsize=15)
# Построим график зависимости цены на недвижимость в зависимости от самых значимых признаков (значимость больше 0.01)
for i in feature_importance['feature'][feature_importance['importances']> 0.01]:
fig, ax = plt.subplots(figsize=(10, 10))
ax.scatter(x = work_data[i], y = work_data['price'])
plt.ylabel('Price', fontsize=13)
plt.xlabel(i, fontsize=13)
plt.show()
Из полезного на графике №2 наблюдается 1 выброс, который стоит удалить, дабы он не повлиял на модели, чувствительные к выбросам, если мы в итоге такие будем использовать
По графикам зависимости от широты и долготы видно, что в дата сете присуствует порядка 10 объектов, у которых местоположение сильно влисяет на стоимость, возможно это условная Luxury недвижимость в элитных районах, где продается не так много объектов ежегодно.
# Удалим выброс и построим график для проверки снова
work_data = work_data[work_data['sqft_living'] < 13000]
fig, ax = plt.subplots(figsize=(10, 10))
ax.scatter(x = work_data['sqft_living'], y = work_data['price'])
plt.ylabel('Price', fontsize=13)
plt.xlabel('Sqft living Area', fontsize=13)
plt.show()
#Посмотрим на количество уникальных значений признаков
for col in work_data.select_dtypes(exclude=['datetime64']).columns:
print([col], len(work_data[col].unique()))
Из интересного здесь.
Категориальные признаки:
Ранговые:
эту информацию можно будет использовать для генерации новых признаков
corrmat = work_data.corr()
f, ax = plt.subplots(figsize=(15, 10))
sns.heatmap(corrmat, square=True);
Из интересного здесь можно увидеть несколько интересных корреляционных групп.
Достаточно сильная корреляция между собой и ценой у:
- bedroom, bathroom, sqft_living
- waterfront, view
- grade, sqft_above
fig, ax = plt.subplots()
fig.set_size_inches(8, 8)
sns.distplot(work_data['price'] , fit=norm);
(mu, sigma) = norm.fit(data['price'])
print( '\n mu = {:.2f} and sigma = {:.2f}\n'.format(mu, sigma))
#Построим график распределения
plt.legend(['Нормальное распределение ($\mu=$ {:.2f} and $\sigma=$ {:.2f} )'.format(mu, sigma)],
loc='best')
plt.ylabel('Частота цены')
plt.title('Распределение цены')
#Дополнительно построим QQ-plot
fig = plt.figure(figsize=(8, 8))
res = stats.probplot(work_data['price'], plot=plt)
plt.show()
На 1-ом графике видно, что распределение цены имеет правый уклон, что на самомм деле достаточно типично, по моему опытуЭ для рынка недвижимости США. По крайней мере 5-6 датасетов, что мне встречались имели похожий тренд.
Уже из этих грфиков, в принципе можно предположить,что линейные модели, которые заточены под нормальное распределение, будут иметь, скорее всего не очень высокий результат без преобразования данных.
df_train = work_data.drop([ 'id','price', 'date'], axis=1)
target = work_data['price']
X_train, X_test, y_train, y_test = train_test_split(df_train, target, test_size=0.25, random_state=17)
scalar = StandardScaler()
clf = Lasso(random_state=17)
lasso_base = Pipeline([('scalar', scalar), ('Lasso', clf)])
lasso_base.set_params(Lasso__alpha=1).fit(X_train, y_train)
y_pred = lasso_base.predict(X_test)
print(" R2_score: {}, RMSLE: {}".format(lasso_base.score(X_test, y_test),rmsle(y_test, y_pred)))
Унылость результата подтверждает гипотезу.
sk_date = work_data['price'] - min(work_data['price']) + 1e-5
_, lmbda_price = boxcox(sk_date)
target = boxcox1p(work_data['price'], lmbda_price)
#Также проведем туже операцию для данных на валидации
target_val = boxcox1p(validation_data['price'], lmbda_price)
print (lmbda_price)
fig, ax = plt.subplots()
fig.set_size_inches(8, 8)
sns.distplot(target , fit=norm);
(mu, sigma) = norm.fit(work_data['price'])
print( '\n mu = {:.2f} and sigma = {:.2f}\n'.format(mu, sigma))
#Построим график распределения
plt.legend(['Нормальное распределение ($\mu=$ {:.2f} and $\sigma=$ {:.2f} )'.format(mu, sigma)],
loc='best')
plt.ylabel('Частота цены')
plt.title('Распределение цены')
#Дополнительно построим QQ-plot
fig = plt.figure(figsize=(8, 8))
res = stats.probplot(target, plot=plt)
plt.show()
Почти прекрасно
#Попробуем нормализавать распределение через геометрическое среднее
fig, ax = plt.subplots()
fig.set_size_inches(8, 8)
geom = (work_data['price']**0.104 - 1 )/ 0.104 * gmean(work_data['price'], axis=0)**(0.104 - 1)
sns.distplot(geom , fit=norm);
plt.ylabel('Частота цены')
plt.title('Распределение цены')
В данном случае, просто проверил формулу, которую видел в учебниках по статистике
# Посмотрим на распределение других количественных признаков с широким диапазоном значений
columns_disc = ['sqft_living', 'sqft_above', 'sqft_basement', 'sqft_living15', 'sqft_lot15', 'sqft_lot']
for col in columns_disc:
fig, ax = plt.subplots()
fig.set_size_inches(8, 8)
sns.distplot(work_data[col] , fit=norm);
(mu, sigma) = norm.fit(work_data[col])
#Построим график распределения
plt.legend(['Нормальное распределение ($\mu=$ {:.2f} and $\sigma=$ {:.2f} )'.format(mu, sigma)],
loc='best')
plt.ylabel('Частота')
plt.title('Распределение')
Очевидно что все вышеуказанные колличественные признаки также как и цена имеют правый уклон.
_lot15, _lot, _basement - имеют много нулевых значений. Возможн стоит дополнительно добавить преобразованные дискретные признаки "Да/Нет" от них в набор данных
skewed_feats = work_data[columns_disc].apply(lambda x: skew(x)).sort_values(ascending=False)
#col = list(data.drop(['date', 'id', 'zipcode','price', 'lat', 'long'], axis=1).columns)
#skewed_feats = data[col].apply(lambda x: skew(x)).sort_values(ascending=False)
print("\nУклон количественных признаков: \n")
skewness = pd.DataFrame({'Skew' :skewed_feats})
skewness.head(6)
Для нормально распределенных признаков значение skew должно быть равно 0 или быть близко к нему.
Для того, чтобы избежать во входных данных появления отрицательных или равных нулю значений, всегда будем находить минимальное значение входной последовательности и вычитать его из каждого ее элемента, дополнительно осуществляя сдвиг на небольшую величину, равную 1e-5.
#Вычислим оптимальное lambda для каждого
skewed_features = skewness.index
lambd = []
for feat in skewed_features:
work_data[feat] = work_data[feat] - min(work_data[feat]) + 1e-5
_, lmbda = boxcox(work_data[feat])
lambd.append(lmbda)
print(feat, lmbda)
sum([abs(x) for x in lambd])/6
lam = 0.12
for feat in skewed_features:
work_data[feat] = work_data[feat] - min(work_data[feat]) + 1e-5
work_data[feat] = boxcox1p(work_data[feat], lam)
validation_data[feat] = validation_data[feat] - min(validation_data[feat]) + 1e-5
validation_data[feat] = boxcox1p(validation_data[feat], lam)
for col in columns_disc:
fig, ax = plt.subplots()
fig.set_size_inches(8, 8)
sns.distplot(work_data[col] , fit=norm);
(mu, sigma) = norm.fit(work_data[col])
#Построим график распределения
plt.legend(['Нормальное распределение ($\mu=$ {:.2f} and $\sigma=$ {:.2f} )'.format(mu, sigma)],
loc='best')
plt.ylabel('Частота')
plt.title('Распределение')
Кроме площади цокольного этажа почти все признаки стали более нормально распределенными.
df_train = work_data.drop(['id','price', 'date'], axis=1)
X_train, X_test, y_train, y_test = train_test_split(df_train, target, test_size=0.25, random_state=17)
clf = RandomForestRegressor(random_state=17)
rnd_base = Pipeline([('scalar', scalar), ('rnd', clf)])
rnd_base.fit(X_train, y_train)
y_pred = rnd_base.predict(X_test)
print("R2_score: {}, RMSLE: {}".format(rnd_base.score(X_test, y_test),
rmsle(inv_boxcox(y_test, lmbda_price), inv_boxcox(y_pred, lmbda_price))))
Есть улучшение относительно исходных данных (0.8504728683309748, 0.19446848336861014) - улучшились оба показателя.
Посмотрим что произошло с важностью признаков.
feature_importance = pd.DataFrame()
feature_importance['feature'] = df_train.columns
feature_importance['importances'] = clf.feature_importances_
feature_importance.sort_values('importances', ascending=False).head(20)
Интересно. Широта обогнала жилую площадь.
# Посмотрим как улучшилась регрессионная модель.
df_train = work_data.drop(['id','price', 'date'], axis=1)
X_train, X_test, y_train, y_test = train_test_split(df_train, target, test_size=0.25, random_state=17)
scalar = StandardScaler()
clf = Lasso(random_state=17)
lasso_base = Pipeline([('scalar', scalar), ('Lasso', clf)])
lasso_base.set_params(Lasso__alpha=1).fit(X_train.values, y_train)
y_pred = lasso_base.predict(X_test.values)
print("R2_score: {}, RMSLE: {}".format(lasso_base.score(X_test, y_test),
rmsle(inv_boxcox(y_test, lmbda_price), inv_boxcox(y_pred, lmbda_price))))
Показатели регрессионной модели ухудшились.
Было: R2_score: 0.711325562859866, RMSLE: 0.39307021980748347
Вообще это странно - такое ухудшение. Возможно, модель просто недообучается с таким высоким значением alpha
# Проверим коэффициенты
print("Коэффициенты: {}\n всего коэффициентов: {}\n всего коэффициентов отличныхот 0: {}".format(clf.coef_, len(clf.coef_),
len(clf.coef_[clf.coef_ != 0])))
Что и следовало доказать, высокий коэффициент регуляризации занижает много весов, оставляя только 2
#Уменьшим регуляризацию
lasso_base = Pipeline([('scalar', scalar), ('Lasso', clf)])
lasso_base.set_params(Lasso__alpha=0.1).fit(X_train.values, y_train)
y_pred = lasso_base.predict(X_test.values)
print("R2_score: {}, RMSLE: {}".format(lasso_base.score(X_test, y_test),
rmsle(inv_boxcox(y_test, lmbda_price), inv_boxcox(y_pred, lmbda_price))))
print("Коэффициенты: {}\n всего коэффициентов: {}\n всего коэффициентов отличныхот 0: {}".format(clf.coef_, len(clf.coef_),
len(clf.coef_[clf.coef_ != 0])))
# Построим кривую кросс-валидации
cv_scores, holdout_scores = [], []
alphas = np.linspace(0.001, 1.0, 15)
for i in alphas:
lasso = Lasso(alpha=i)
cv_scores.append(np.mean(cross_val_score(lasso, X_train, y_train, cv=3)))
lasso.fit(X_train, y_train)
holdout_scores.append(r2_score(y_test, lasso.predict(X_test)))
plt.figure(figsize=(8, 8))
plt.plot(alphas, cv_scores, label='CV')
plt.plot(alphas, holdout_scores, label='holdout')
plt.title('Easy task. Lasso fails')
plt.legend();
Лучше оставить такой alpha = 0.1, так гораздо больше признаков работает.
Займемся поиском точек роста качества модели за счет инжиниринга признаков. В части 3, где мы проводили визуальный анализ было обнаружено несколько интересныз закономерностей, которые мы инкрементируем в модель данных и посмотрим, как они повлияют на качество обучения.
# Добавим дополнительные дискретные признаки для колличественных, где большая часть имеет нулевые значения
work_data['_lot'] = 0
work_data['_lot15'] = 0
work_data['_basement'] = 0
work_data['_lot'][work_data['sqft_lot']>0] = 1
work_data['_lot15'][work_data['sqft_lot15']>0] = 1
work_data['_basement'][work_data['sqft_basement']>0] = 1
# Выделю 3 группы объектов исходя из их стоимости
work_data['chip'] = 0
work_data['mid'] = 0
work_data['exp'] = 0
work_data['chip'][work_data['price']<= 375000] = 1
work_data['mid'][(work_data['price']> 375000)&(work_data['price']<= 1000000)] = 1
work_data['exp'][work_data['price']> 1000000] = 1
#Выделю сезоны продажи недвижимости
work_data['sum'] = 0
work_data['spr'] = 0
work_data['aut'] = 0
work_data['win'] = 0
work_data['sum'][(work_data['month_sold']>= 6) & (work_data['month_sold']<= 8)] = 1
work_data['spr'][(work_data['month_sold']>= 3) & (work_data['month_sold']<= 5)] = 1
work_data['aut'][(work_data['month_sold']>= 9) & (work_data['month_sold']<= 11)] = 1
work_data['win'][(work_data['month_sold']>= 1) & (work_data['month_sold']<= 2) | (work_data['month_sold']== 12)] = 1
# Выделю клас жилья через грейд
work_data['low'] = 0
work_data['average'] = 0
work_data['high'] = 0
work_data['low'][work_data['grade']<= 5] = 1
work_data['average'][(work_data['price']> 5)&(work_data['price']<= 9)] = 1
work_data['high'][work_data['price']> 9] = 1
#Выделю 3 группы исходя из года постройки
work_data['old'] = 0
work_data['old_to_hihg'] = 0
work_data['young'] = 0
work_data['old'][work_data['yr_built']<= 1942] = 1
work_data['old_to_hihg'][(work_data['yr_built']> 1942)&(work_data['yr_built']<= 1994)] = 1
work_data['young'][work_data['yr_built']> 1994] = 1
#Добавим 4 новых признака, 2 характеризующих "возраст" дома, наличие подвала "Да\Нет", была ли реноваци "Да\нет"
work_data['years_old'] = 2015 - work_data['yr_built']
work_data['second_young'] = 2015 - work_data['yr_renovated']
work_data['second_young'][work_data['second_young'] == 2015] = 0
work_data['basement'] = 0
work_data['basement'][work_data['sqft_basement'] >0] = 1
work_data['was_renov'] = 0
work_data['was_renov'][work_data['yr_renovated'] >0] = 1
#Напишем простую функцию для лучшего разделения признаков, характеризующих размеры
def quadratic(feature):
work_data[feature+'2'] = work_data[feature]**2
quadratic('sqft_living')
quadratic('sqft_living15')
quadratic('sqft_lot')
quadratic('sqft_above')
quadratic('sqft_basement')
# Проделю все тоже самое и для контрольной выборки
validation_data['_lot'] = 0
validation_data['_lot15'] = 0
validation_data['_basement'] = 0
validation_data['_lot'][validation_data['sqft_lot']>0] = 1
validation_data['_lot15'][validation_data['sqft_lot15']>0] = 1
validation_data['_basement'][validation_data['sqft_basement']>0] = 1
validation_data['chip'] = 0
validation_data['mid'] = 0
validation_data['exp'] = 0
validation_data['chip'][validation_data['price']<= 375000] = 1
validation_data['mid'][(validation_data['price']> 375000)&(validation_data['price']<= 1000000)] = 1
validation_data['exp'][validation_data['price']> 1000000] = 1
validation_data['sum'] = 0
validation_data['spr'] = 0
validation_data['aut'] = 0
validation_data['win'] = 0
validation_data['sum'][(validation_data['month_sold']>= 6) & (validation_data['month_sold']<= 8)] = 1
validation_data['spr'][(validation_data['month_sold']>= 3) & (validation_data['month_sold']<= 5)] = 1
validation_data['aut'][(validation_data['month_sold']>= 9) & (validation_data['month_sold']<= 11)] = 1
validation_data['win'][(validation_data['month_sold']>= 1) & \
(validation_data['month_sold']<= 2) | (validation_data['month_sold']== 12)] = 1
validation_data['low'] = 0
validation_data['average'] = 0
validation_data['high'] = 0
validation_data['low'][validation_data['grade']<= 5] = 1
validation_data['average'][(validation_data['price']> 5)&(validation_data['price']<= 9)] = 1
validation_data['high'][validation_data['price']> 9] = 1
validation_data['old'] = 0
validation_data['old_to_hihg'] = 0
validation_data['young'] = 0
validation_data['old'][validation_data['yr_built']<= 1942] = 1
validation_data['old_to_hihg'][(validation_data['yr_built']> 1942)&(validation_data['yr_built']<= 1994)] = 1
validation_data['young'][validation_data['yr_built']> 1994] = 1
# Тоже самое для валидации
validation_data['years_old'] = 2015 - validation_data['yr_built']
validation_data['second_young'] = 2015 - validation_data['yr_renovated']
validation_data['second_young'][validation_data['second_young'] == 2015] = 0
validation_data['basement'] = 0
validation_data['basement'][validation_data['sqft_basement'] >0] = 1
validation_data['was_renov'] = 0
validation_data['was_renov'][validation_data['yr_renovated'] >0] = 1
def quadratic(feature):
validation_data[feature+'2'] = validation_data[feature]**2
quadratic('sqft_living')
quadratic('sqft_living15')
quadratic('sqft_lot')
quadratic('sqft_above')
quadratic('sqft_basement')
work_data.shape, validation_data.shape
df_train = work_data.drop(['id','price', 'date'], axis=1)
X_train, X_test, y_train, y_test = train_test_split(df_train, target, test_size=0.25, random_state=17)
clf = RandomForestRegressor(random_state=17)
rnd_base = Pipeline([('scalar', scalar), ('rnd', clf)])
rnd_base.fit(X_train, y_train)
y_pred = rnd_base.predict(X_test)
print("R2_score: {}, RMSLE: {}".format(rnd_base.score(X_test, y_test),
rmsle(inv_boxcox(y_test, lmbda_price), inv_boxcox(y_pred, lmbda_price))))
Есть неплохой прирост (R2_score: 0.8757049690628793, RMSLE: 0.18806603652440654)
clf = Lasso(random_state=17)
lasso_base = Pipeline([('scalar', scalar), ('Lasso', clf)])
lasso_base.set_params(Lasso__alpha=0.1).fit(X_train.values, y_train)
y_pred = lasso_base.predict(X_test.values)
print("R2_score: {}, RMSLE: {}".format(lasso_base.score(X_test, y_test),
rmsle(inv_boxcox(y_test, lmbda_price), inv_boxcox(y_pred, lmbda_price))))
Регрессия тоже растет
*Данные были подготовлены отдельно через обработку сырых таблиц и отчетности.
Уровень преступности - агрегированан из отчета с сайта шерифа округа, где фиксируются все преступления по времени, по типу и месту
Школы - агрегирован, если правильно помню, из отчета по вакцинации учебных заведений округа (нач. школы и дет. сады).
Клиники - таблица с дислокацией мед. учереждений
Магазины, кафе, рестораны - агрегирован из отчета испекторов санитарных служб, которые занимаюца лецензированием таких заведений.
Сами данные брались из публичных источников и обрабатывались для выявления местонахождения объектов. Идея состоит в том, что близость к определенным объектам, наличие объектов для базовых потребностей(школа, больница, магазин), скопление объектов или событий(криминал) - все это влияет на конечную стоимость жилья. Проверим...
# Уровень преступности
crime = pd.read_csv('crimeKC.csv', index_col=0)
# Местоположения по индексу образовательных учереждений (нач. школа и дет. сады)
school = pd.read_csv('schoolKC.csv', index_col=0)
# Местоположения по индексу медицинских учереждений
clinic = pd.read_csv('clinicsKC.csv', index_col=0)
# Местоположения по индексу продуктовых магазинов, кафе, ресторанов
shop_horeca = pd.read_csv('shop_cafeKC.csv', index_col=0)
crime.head()
school.head()
clinic.head()
shop_horeca.head()
Новые данные представляют из себя 4 разреженных матрицы, которые я сгенерировал из различных публичных государственных источников округа Кинг
В данной таблице собраны суммарные происшествия за несколько лет по типам и по zipcode. Несмотря на то что данные агрегированны за период больший периоду продаж недвижимости, они скорее отражают тренд, выделяя более "преступные" локации.
work_data = work_data.merge(crime, how='left', on='zipcode')
#Некоторые данные не будут иметь совпадений из таблицы crime
work_data = work_data.fillna(0)
work_data.head()
%%time
# Проверим регрессию
df_train = work_data.drop(['id', 'price', 'date'], axis=1)
X_train, X_test, y_train, y_test = train_test_split(df_train, target, test_size=0.25, random_state=17)
scalar = StandardScaler()
clf = Lasso(random_state=17)
lasso_base = Pipeline([('scalar', scalar), ('Lasso', clf)])
lasso_base.set_params(Lasso__alpha=0.1).fit(X_train, y_train)
y_pred = lasso_base.predict(X_test)
print("R2_score: {}, RMSLE: {}".format(lasso_base.score(X_test, y_test),
rmsle(inv_boxcox(y_test, lmbda_price), inv_boxcox(y_pred, lmbda_price))))
Напомню до этого результат был: R2_score: 0.859201120304148, RMSLE: 0.1985696581535856 - т.е. новые данные улучшили модель.
%%time
# Сравним с RandomForest
clf = RandomForestRegressor(random_state=17)
rnd_base = Pipeline([('scalar', scalar), ('rnd', clf)])
rnd_base.fit(X_train, y_train)
y_pred = rnd_base.predict(X_test)
print("R2_score: {}, RMSLE: {}".format(rnd_base.score(X_test, y_test),
rmsle(inv_boxcox(y_test, lmbda_price), inv_boxcox(y_pred, lmbda_price))))
RandomForestRegressor немного ухудшился от прошлого значения: R2_score: 0.9169665913882967, RMSLE: 0.15419440365697273
work_data = work_data.merge(clinic, how='left', on='zipcode')
work_data = work_data.fillna(0)
%%time
df_train = work_data.drop(['id', 'price', 'date'], axis=1)
X_train, X_test, y_train, y_test = train_test_split(df_train, target, test_size=0.25, random_state=17)
scalar = StandardScaler()
clf = Lasso(random_state=17)
lasso_base = Pipeline([('scalar', scalar), ('Lasso', clf)])
lasso_base.set_params(Lasso__alpha=0.1).fit(X_train, y_train)
y_pred = lasso_base.predict(X_test)
print("R2_score: {}, RMSLE: {}".format(lasso_base.score(X_test, y_test),
rmsle(inv_boxcox(y_test, lmbda_price), inv_boxcox(y_pred, lmbda_price))))
%%time
df_train = work_data.drop(['id','price', 'date'], axis=1)
X_train, X_test, y_train, y_test = train_test_split(df_train, target, test_size=0.25, random_state=17)
clf = RandomForestRegressor(random_state=17)
rnd_base = Pipeline([('scalar', scalar), ('rnd', clf)])
rnd_base.fit(X_train, y_train)
y_pred = rnd_base.predict(X_test)
print("R2_score: {}, RMSLE: {}".format(rnd_base.score(X_test, y_test),
rmsle(inv_boxcox(y_test, lmbda_price), inv_boxcox(y_pred, lmbda_price))))
Снова небольшой рост на обоих моделях.
work_data = work_data.merge(school, how='left', on ='zipcode')
work_data = work_data.fillna(0)
%%time
df_train = work_data.drop(['id', 'price', 'date'], axis=1)
X_train, X_test, y_train, y_test = train_test_split(df_train, target, test_size=0.25, random_state=17)
scalar = StandardScaler()
clf = Lasso(random_state=17)
lasso_base = Pipeline([('scalar', scalar), ('Lasso', clf)])
lasso_base.set_params(Lasso__alpha=0.1).fit(X_train, y_train)
y_pred = lasso_base.predict(X_test)
print("R2_score: {}, RMSLE: {}".format(lasso_base.score(X_test, y_test),
rmsle(inv_boxcox(y_test, lmbda_price), inv_boxcox(y_pred, lmbda_price))))
%%time
clf = RandomForestRegressor(random_state=17)
rnd_base = Pipeline([('scalar', scalar), ('rnd', clf)])
rnd_base.fit(X_train, y_train)
y_pred = rnd_base.predict(X_test)
print("R2_score: {}, RMSLE: {}".format(rnd_base.score(X_test, y_test),
rmsle(inv_boxcox(y_test, lmbda_price), inv_boxcox(y_pred, lmbda_price))))
work_data = work_data.merge(shop_horeca, how='left', on ='zipcode')
work_data = work_data.fillna(0)
%%time
df_train = work_data.drop(['id', 'price', 'date'], axis=1)
X_train, X_test, y_train, y_test = train_test_split(df_train, target, test_size=0.25, random_state=17)
scalar = StandardScaler()
clf = Lasso(random_state=17)
lasso_base = Pipeline([('scalar', scalar), ('Lasso', clf)])
lasso_base.set_params(Lasso__alpha=0.1).fit(X_train, y_train)
y_pred = lasso_base.predict(X_test)
print("R2_score: {}, RMSLE: {}".format(lasso_base.score(X_test, y_test),
rmsle(inv_boxcox(y_test, lmbda_price), inv_boxcox(y_pred, lmbda_price))))
%%time
clf = RandomForestRegressor(random_state=17)
rnd_base = Pipeline([('scalar', scalar), ('rnd', clf)])
rnd_base.fit(X_train, y_train)
y_pred = rnd_base.predict(X_test)
print("R2_score: {}, RMSLE: {}".format(rnd_base.score(X_test, y_test),
rmsle(inv_boxcox(y_test, lmbda_price), inv_boxcox(y_pred, lmbda_price))))
# Воспользуемся функцией одной из лекции курса
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)
%%time
steps = [('scalar', StandardScaler()),
('Lasso', Lasso(random_state=17))]
pipeline = Pipeline(steps)
parameters = {'Lasso__alpha': [0.001, 0.01, 0.1, 1]}
cvLasso = GridSearchCV(pipeline,param_grid=parameters,cv=3)
cvLasso.fit(X_train.values,y_train)
print("Accuracy: {}".format(cvLasso.score(X_test, y_test)))
print("Tuned Model Parameters: {}".format(cvLasso.best_params_))
cvLasso.cv_results_
#Посмотрим, что у нас получилось на кросс-валидации для Lasso.
plt.figure(figsize=(8, 8))
test_results = np.array([cvLasso.cv_results_['split0_test_score'],cvLasso.cv_results_['split1_test_score'],
cvLasso.cv_results_['split2_test_score']]).T;
train_results = np.array([cvLasso.cv_results_['split0_train_score'],cvLasso.cv_results_['split1_train_score'],
cvLasso.cv_results_['split2_train_score']]).T;
plot_with_err([0.001, 0.01, 0.1, 1], test_results, label='validation scores')
plot_with_err([0.001, 0.01, 0.1, 1], train_results, label='train scores')
plt.xlabel('alphas'); plt.ylabel('R2')
plt.legend()
Качество падает линейно с ростом коэффициента регуляризации. Регрессия очень не стабильна.
steps = [('scalar', StandardScaler()),
('RandomForestRegressor', RandomForestRegressor(random_state=17))]
pipeline = Pipeline(steps)
parameters = {'RandomForestRegressor__n_estimators': [50, 100, 150],
'RandomForestRegressor__max_depth' : [15, 20, 25]}
cvRf = GridSearchCV(pipeline,param_grid=parameters,cv=3)
cvRf.fit(X_train,y_train)
print("Accuracy: {}".format(cvRf.score(X_test, y_test)))
print("Tuned Model Parameters: {}".format(cvRf.best_params_))
Внимание: при n_estimators более 250 (однажды выставил такое значение) - Считалось около 3.5 часов. Случайный лес с повышением при увеличение числа деревьев и максимальной глубины, наращивает качество, но становится очень медленным и начинает проигрывать по скорости регрессии очень сильно. Вообще я запускал перебор только 2 раза, а признаки добавлял новые гораздо чаще, так что при запуске score будет отличным от текущего.
cvRf.cv_results_
#Посмотрим, что у нас получилось на кросс-валидации для RandomForest.
plt.figure(figsize=(8, 8))
test_results = np.array([cvRf.cv_results_['split0_test_score'],cvRf.cv_results_['split1_test_score'],
cvRf.cv_results_['split2_test_score']]).T;
train_results = np.array([cvRf.cv_results_['split0_train_score'],cvRf.cv_results_['split1_train_score'],
cvRf.cv_results_['split2_train_score']]).T;
plot_with_err(['1-15_50', '2-15_100', '3-15_150', '4-20_50', '5-20_100', '6-20_150',
'7-25_50', '8-25_100', '9-25_150'], test_results, label='validation scores')
plot_with_err(['1-15_50', '2-15_100', '3-15_150', '4-20_50', '5-20_100', '6-20_150',
'7-25_50', '8-25_100', '9-25_150'], train_results, label='train scores')
plt.xlabel('depth-n_estiators'); plt.ylabel('R2')
plt.legend()
RandomForest намного стабильнее Lasso, причем уже при n_estimators = 100 и глубине = 15 уже дает отличный результат
steps = [('scalar', StandardScaler()),
('lgb.LGBMRegressor', lgb.LGBMRegressor(random_state=17))]
pipeline = Pipeline(steps)
parameters = {'lgb.LGBMRegressor__n_estimators': [700, 800,900],
'lgb.LGBMRegressor__max_depth' : [3, 4, 5],
'lgb.LGBMRegressor__learning_rate' : [0.03, 0.06, 0.09]}
cvLgb = GridSearchCV(pipeline,param_grid=parameters,cv=3)
cvLgb.fit(X_train,y_train)
print("Accuracy: {}".format(cvLgb.score(X_test, y_test)))
print("Tuned Model Parameters: {}".format(cvLgb.best_params_))
cvLgb.cv_results_
#Посмотрим, что у нас получилось на кросс-валидации для lgb.
plt.figure(figsize=(15, 15))
test_results = np.array([cvLgb.cv_results_['split0_test_score'],cvLgb.cv_results_['split1_test_score'],
cvLgb.cv_results_['split2_test_score']]).T;
train_results = np.array([cvLgb.cv_results_['split0_train_score'],cvLgb.cv_results_['split1_train_score'],
cvLgb.cv_results_['split2_train_score']]).T;
plot_with_err(['11-0.03_3_700', '11-0.03_3_800', '11-0.03_3_900', '12-0.03_4_700', '12-0.03_4_800', '12-0.03_4_900',
'13-0.03_5_700', '13-0.03_5_800', '13-0.03_5_900', '21-0.06_3_700', '21-0.06_3_800', '21-0.06_3_900',
'22-0.06_4_700', '22-0.06_4_800', '22-0.06_4_900','23-0.06_5_700', '23-0.06_5_800', '23-0.06_5_900',
'31-0.09_3_700', '31-0.09_3_800', '31-0.09_3_900', '32-0.09_4_700', '32-0.09_4_800', '32-0.09_4_900',
'33-0.09_5_700', '33-0.09_5_800', '33-0.09_5_900'], test_results, label='validation scores')
plot_with_err(['11-0.03_3_700', '11-0.03_3_800', '11-0.03_3_900', '12-0.03_4_700', '12-0.03_4_800', '12-0.03_4_900',
'13-0.03_5_700', '13-0.03_5_800', '13-0.03_5_900', '21-0.06_3_700', '21-0.06_3_800', '21-0.06_3_900',
'22-0.06_4_700', '22-0.06_4_800', '22-0.06_4_900','23-0.06_5_700', '23-0.06_5_800', '23-0.06_5_900',
'31-0.09_3_700', '31-0.09_3_800', '31-0.09_3_900', '32-0.09_4_700', '32-0.09_4_800', '32-0.09_4_900',
'33-0.09_5_700', '33-0.09_5_800', '33-0.09_5_900'], train_results, label='train scores')
plt.xlabel('learning rate - depth - n_estiators'); plt.ylabel('R2')
plt.xticks(rotation= 90)
plt.legend()
%%time
clf = RandomForestRegressor(random_state=17)
rf_base = Pipeline([('scalar', scalar), ('RandomForestRegressor', clf)])
rf_base.set_params(RandomForestRegressor__n_estimators=150, RandomForestRegressor__max_depth=15).fit(X_train, y_train)
y_pred = lgb_base.predict(X_test)
print("R2_score: {}, RMSLE: {}".format(rf_base.score(X_test, y_test),
rmsle(inv_boxcox(y_test, lmbda_price), inv_boxcox(y_pred, lmbda_price))))
Неплохо, но по времени он протигрывает lgb
%%time
clf = lgb.LGBMRegressor(random_state=17)
lgb_base = Pipeline([('scalar', scalar), ('lgb', clf)])
lgb_base.set_params(lgb__learning_rate=0.06, lgb__n_estimators=900, lgb__max_depth=4).fit(X_train, y_train)
y_pred = lgb_base.predict(X_test)
print("R2_score: {}, RMSLE: {}".format(lgb_base.score(X_test, y_test),
rmsle(inv_boxcox(y_test, lmbda_price), inv_boxcox(y_pred, lmbda_price))))
Лучший результат и отличное время
%%time
clf = Lasso(random_state=17)
lasso_base = Pipeline([('scalar', scalar), ('Lasso', clf)])
lasso_base.set_params(Lasso__alpha=0.001).fit(X_train, y_train)
y_pred = lasso_base.predict(X_test)
print("R2_score: {}, RMSLE: {}".format(lasso_base.score(X_test, y_test),
rmsle(inv_boxcox(y_test, lmbda_price), inv_boxcox(y_pred, lmbda_price))))
Правда, я не очень верю в стабильность регрессии с alpha = 0.001 - скорее всего на отложенной выборке такая регуляризация будет пагубно влиять на результат, если отложенные данные имеют хоть немного отличное статистическое распределение.
Воспользуемся авторской реализацие стеккинга и блендинга от Александра Дьяконова отсюда:
https://github.com/Dyakonov/ml_hacks/blob/master/dj_stacking.ipynb
Отличная статья на эту же тему от Александра Дьяконова:
# Класс для стеккинга и блендинга
class DjStacking(BaseEstimator, ClassifierMixin):
"""Стэкинг моделей scikit-learn"""
def __init__(self, models, ens_model):
"""
Инициализация
models - базовые модели для стекинга
ens_model - мета-модель
"""
self.models = models
self.ens_model = ens_model
self.n = len(models)
self.valid = None
def fit(self, X, y=None, p=0.25, cv=3, err=0.001, random_state=17):
"""
Обучение стекинга
p - в каком отношении делить на обучение / тест
если p = 0 - используем всё обучение!
cv (при p=0) - сколько фолдов использовать
err (при p=0) - величина случайной добавки к метапризнакам
random_state - инициализация генератора
"""
if (p > 0): # делим на обучение и тест
# разбиение на обучение моделей и метамодели
train, valid, y_train, y_valid = train_test_split(X, y, test_size=p, random_state=random_state)
# заполнение матрицы для обучения метамодели
self.valid = np.zeros((valid.shape[0], self.n))
for t, clf in enumerate(self.models):
clf.fit(train, y_train)
self.valid[:, t] = clf.predict(valid)
# обучение метамодели
self.ens_model.fit(self.valid, y_valid)
else: # используем всё обучение
# для регуляризации - берём случайные добавки
self.valid = err*np.random.randn(X.shape[0], self.n)
for t, clf in enumerate(self.models):
# это oob-ответы алгоритмов
self.valid[:, t] += cross_val_predict(clf, X, y, cv=cv, n_jobs=-1, method='predict')
# но сам алгоритм надо настроить
clf.fit(X, y)
# обучение метамодели
self.ens_model.fit(self.valid, y)
return self
def predict(self, X, y=None):
"""
Работа стэкинга
"""
# заполение матрицы для мета-классификатора
X_meta = np.zeros((X.shape[0], self.n))
for t, clf in enumerate(self.models):
X_meta[:, t] = clf.predict(X)
a = self.ens_model.predict(X_meta)
return (a)
%%time
# Обучу различные модели с разными гиперпараметрами и посмотрю на их отдельные результаты.
df_train = work_data.drop(['id','price', 'date'], axis=1)
X_train, X_test, y_train, y_test = train_test_split(df_train, target, test_size=0.25, random_state=17)
las0 = Lasso(alpha=0.001)
las0.fit(X_train, y_train)
y_pred = las0.predict(X_test)
print("las0, R2_score: {}, RMSLE: {}".format(las0.score(X_test, y_test),
rmsle(inv_boxcox(y_test, lmbda_price), inv_boxcox(y_pred, lmbda_price))))
las1 = Lasso(alpha=0.1)
las1.fit(X_train, y_train)
y_pred = las1.predict(X_test)
print("las1, R2_score: {}, RMSLE: {}".format(las1.score(X_test, y_test),
rmsle(inv_boxcox(y_test, lmbda_price), inv_boxcox(y_pred, lmbda_price))))
las2 = Lasso(alpha=0.01)
las2.fit(X_train, y_train)
y_pred = las2.predict(X_test)
print("las2, R2_score: {}, RMSLE: {}".format(las2.score(X_test, y_test),
rmsle(inv_boxcox(y_test, lmbda_price), inv_boxcox(y_pred, lmbda_price))))
rf1 = RandomForestRegressor(n_estimators=50)
rf1.fit(X_train, y_train)
y_pred = rf1.predict(X_test)
print("rf1, R2_score: {}, RMSLE: {}".format(rf1.score(X_test, y_test),
rmsle(inv_boxcox(y_test, lmbda_price), inv_boxcox(y_pred, lmbda_price))))
rf2 = RandomForestRegressor(n_estimators=250, max_depth=20)
rf2.fit(X_train, y_train)
y_pred = rf2.predict(X_test)
print("rf2, R2_score: {}, RMSLE: {}".format(rf2.score(X_test, y_test),
rmsle(inv_boxcox(y_test, lmbda_price), inv_boxcox(y_pred, lmbda_price))))
gbm1 = lgb.LGBMRegressor(boosting_type='gbdt', learning_rate=0.06, max_depth=4, n_estimators=900, nthread=-1, objective='regression')
gbm1.fit(X_train, y_train)
y_pred = gbm1.predict(X_test)
print("gbm1, R2_score: {}, RMSLE: {}".format(gbm1.score(X_test, y_test),
rmsle(inv_boxcox(y_test, lmbda_price), inv_boxcox(y_pred, lmbda_price))))
gbm2 = lgb.LGBMRegressor(boosting_type='gbdt', learning_rate=0.05, max_depth=3, n_estimators=800, nthread=-1, objective='regression')
gbm2.fit(X_train, y_train)
y_pred = gbm2.predict(X_test)
print("gbm2, R2_score: {}, RMSLE: {}".format(gbm2.score(X_test, y_test),
rmsle(inv_boxcox(y_test, lmbda_price), inv_boxcox(y_pred, lmbda_price))))
gbm3 = lgb.LGBMRegressor(boosting_type='gbdt', learning_rate=0.05, max_depth=5, n_estimators=1000, nthread=-1, objective='regression')
gbm3.fit(X_train, y_train)
y_pred = gbm3.predict(X_test)
print("gbm3, R2_score: {}, RMSLE: {}".format(gbm3.score(X_test, y_test),
rmsle(inv_boxcox(y_test, lmbda_price), inv_boxcox(y_pred, lmbda_price))))
# Набор моделей, в качестве мета-алгоритма возьму Lasso регрессию.
models = [las2,las0, rf1, rf2, gbm1, gbm2, gbm3]
ens_model = Lasso(alpha = 0.1, random_state = 17)
# Запустим блендинг и стеккинг
s1 = DjStacking(models, ens_model)
s1.fit(X_train, y_train)
y_pred1=s1.predict(X_test)
print("R2_score: {}, RMSLE: {}".format(r2_score(y_test, y_pred1),
rmsle(inv_boxcox(y_test, lmbda_price), inv_boxcox(y_pred1, lmbda_price))))
s2 = DjStacking(models, ens_model)
s2.fit(X_train, y_train, p=-1)
y_pred2 = s2.predict(X_test)
print("R2_score: {}, RMSLE: {}".format(r2_score(y_test, y_pred2),
rmsle(inv_boxcox(y_test, lmbda_price), inv_boxcox(y_pred2, lmbda_price))))
Стеккинг и блендинг не позволили улучшить качество модели. И в этом нет ничего не обычного, посокльку у нас среди моделей есть очень хорошая модель, которая дает очень высокий результат. Скорее всего улучшений можно было бы добиться больше "лучших" разннобразных моделей, которые в итоге бы дали небольшой прирост качства. Но, в целом, практическая польза от такого метода возможна, разве что, в соревнованиях, поскольку с ростом числа объектов в выборке и роста числа базовых тяжелых алгоритмов скорость работы будет еще больше увеличиваться относительно 1 алгоритма, например того же lgb. Еще, как вариант, можно исключать из стеккинга "тяжелые алгоритмы" и пробовать добиться высокого качества на простых, но польза от такого опять же только в соревнованиях.
# Мержим дополнительные таблицы
validation_data = validation_data.merge(clinic, how='left', on='zipcode')
validation_data = validation_data.merge(shop_horeca, how='left', on ='zipcode')
validation_data = validation_data.merge(school, how='left', on ='zipcode')
validation_data = validation_data.merge(crime, how='left', on='zipcode')
validation_data = validation_data.fillna(0)
df_train = work_data.drop(['id','price', 'date'], axis=1)
df_test = validation_data.drop(['id','price', 'date'], axis=1)
X_train, X_test, y_train, y_test = train_test_split(df_train, target, test_size=0.25, random_state=17)
"""
Проверяем лучшую модель. Тренируем на всей work_data и предсказываем для validation_data
дополнительно немного переберем параметры, дабы понять насколько мы переобучаемся на тестовых данных,
когда настраивали гиперпараметры
"""
test_R2, holdout_R2, test_rmsle, holdout_rmsle = [], [],[], []
n_estimator = [100,200,300,500,700,900]
for i in n_estimator:
lgb_final = lgb.LGBMRegressor(max_depth=4, lgb__learning_rate=0.06, n_estimators=i)
lgb_final.fit(X_train, y_train)
test_R2.append(r2_score(y_test, lgb_final.predict(X_test)))
test_rmsle.append(rmsle(inv_boxcox(y_test, lmbda_price), inv_boxcox(lgb_final.predict(X_test), lmbda_price)))
holdout_R2.append(r2_score(target_val, lgb_final.predict(df_test)))
holdout_rmsle.append(rmsle(inv_boxcox(target_val, lmbda_price), inv_boxcox(lgb_final.predict(df_test), lmbda_price)))
plt.figure(figsize=(15, 15))
plt.plot(n_estimator, test_R2, label='test_R2')
plt.plot(n_estimator, holdout_R2, label='holdout_R2')
plt.title('Ошибка lgb на отложенной выборке R2')
plt.legend();
plt.figure(figsize=(15, 15))
plt.plot(n_estimator, test_rmsle, label='test_rmsle')
plt.plot(n_estimator, holdout_rmsle, label='holdout_rmsle')
plt.title('Ошибка lgb на отложенной выборке RSMLE')
plt.legend();
# Теперь посмотрим на качество при изменении глубины
test_R2, holdout_R2, test_rmsle, holdout_rmsle = [], [],[], []
n_estimator = [1,2,3,4,5,6]
for i in n_estimator:
lgb_final = lgb.LGBMRegressor(max_depth=i, lgb__learning_rate=0.06, n_estimators=500)
lgb_final.fit(X_train, y_train)
test_R2.append(r2_score(y_test, lgb_final.predict(X_test)))
test_rmsle.append(rmsle(inv_boxcox(y_test, lmbda_price), inv_boxcox(lgb_final.predict(X_test), lmbda_price)))
holdout_R2.append(r2_score(target_val, lgb_final.predict(df_test)))
holdout_rmsle.append(rmsle(inv_boxcox(target_val, lmbda_price), inv_boxcox(lgb_final.predict(df_test), lmbda_price)))
plt.figure(figsize=(15, 15))
plt.plot(n_estimator, test_R2, label='test_R2')
plt.plot(n_estimator, holdout_R2, label='holdout_R2')
plt.title('Ошибка lgb на отложенной выборке R2')
plt.legend();
plt.figure(figsize=(15, 15))
plt.plot(n_estimator, test_rmsle, label='test_rmsle')
plt.plot(n_estimator, holdout_rmsle, label='holdout_rmsle')
plt.title('Ошибка lgb на отложенной выборке RSMLE')
plt.legend();
gbm_win = lgb.LGBMRegressor(learning_rate=0.06, max_depth=7, n_estimators=500, random_state=17)
gbm_win.fit(X_train, y_train)
y_pred = gbm_win.predict(df_test)
print("gbm_win, R2_score: {}, RMSLE: {}".format(gbm_win.score(df_test, target_val),
rmsle(inv_boxcox(target_val, lmbda_price), inv_boxcox(y_pred, lmbda_price))))