Объем предварительной работы, которую требовал этот проект, составил, пожалуй, один из самых сложных, но полезных опытов, которые я получил этим летом. Благодаря бесчисленным итерациям подхода, извлечения данных и окончательного выбора модели для построения красивых геокарт, предварительной обработки и, наконец, настройки параметров, каждый шаг создавал препятствие, похожее на POC, для преодоления которого требовалось время, но в конечном итоге это привело к плодотворному профессиональному развитию. Удивительно, хотя и полезно, но трехмесячный курс по науке о данных, который я прошел ранее в этом году, меркнет по сравнению с ростом, который я испытал, проводя пятничные вечера, работая над этим проектом в счастливом одиночестве.
Решающим фактором с точки зрения темы является опыт, который я получил в прошлом году в качестве аналитика кредитных рисков в банке. Как и ожидалось, эта роль требовала обработки многих TFSA, RRSP, инвестиционных кредитов, а также заявок на жилье. Это пробудило во мне любопытство к рынку жилья Онтарио, особенно к прогнозированию цен без участия человека. Многие вопросы, например, как разные области соотносятся друг с другом с точки зрения стоимости, цены, демографического распределения и т. д., остались без ответа даже после моего ухода из финансов в прошлом году. Поскольку банковским учреждениям трудно полностью отказаться от своего аналогового/ручного прошлого, работа с новыми эффективными технологиями становится необыкновенной роскошью. С моим нынешним набором навыков я решил, что у меня есть инструменты, чтобы попытаться найти решение.
Особое примечание: проект был основан на реальных данных из базы данных банка. Чтобы избежать каких-либо проблем с конфиденциальностью, я решил воссоздать проект и использовать те же методы, которые я использовал там, но с общедоступными данными. Я воспользовался тем, что правительство Канады ввело принудительные меры в отношении компаний, занимающихся недвижимостью, чтобы раскрыть информацию о продажах жилья для общественности.
Связанные вопросы:
- Источники и средства для получения набора данных о жилье
- Разработка новых предикторов с использованием разных источников
- Основные выводы из текущего рынка жилья
- Выбор модели МЛ
- Лучший подход к оптимизации гиперпараметров
Ключевые выводы проекта:
- Эффективный просмотр веб-страниц
- Самостоятельно построенный набор данных и разработка функций
- Визуализация географических данных
- Стандартизация и мощное преобразование данных
- Обработка недостающих данных
- Создание пользовательских целей и оценок для модели ML
- Автоматическая настройка гиперпараметров
- Усреднение моделей и суммирование моделей
Наборы данных:
- Набор данных о домах ручной работы (временной горизонт — 1 год)
- Географические данные районов GTA
- Профили соседей (2 — года)
- Географические данные линий / станций метро Торонто
Ограничение: несмотря на то, что последний уровень дохода в Торонто был недоступен, реляционной информации из набора данных за 2017 год было достаточно.
# Load dependencies
import warnings
warnings.filterwarnings('ignore')
import pandas as pd
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm
import urllib.request
import time
import bs4 as bs
from geopy.geocoders import Nominatim
import geopandas as gpd
import seaborn as sns
from scipy import stats
# To make pyproj conversions work
import os
os.environ['PROJ_LIB']=r"C:\Users\spiris\AppData\Local\Continuum\anaconda3\Library\share"
Пожалуйста, обратитесь к проекту на GitHub для получения подробной информации и создания набора данных. Я пропущу несколько шагов в начале.
# https://open.toronto.ca/dataset/neighbourhoods/
# https://open.toronto.ca/dataset/neighbourhood-profiles/
t_geodata = gpd.read_file('datasets/opendata_neighbourhoods/Neighbourhoods.shp')
# rename columns
t_geodata_columns = pd.read_csv('datasets/opendata_neighbourhoods/Neighbourhoods_fields.csv')
t_geodata.columns = t_geodata_columns.name.values
# Convert data to lat, long type
t_geodata['geometry'] = t_geodata['geometry'].to_crs({'init': 'EPSG:4269'})
houses = pd.read_csv('houses.csv', index_col='index')
houses = houses.dropna(subset=['city_district'])
houses.drop(index=[14, 14388], inplace=True)
print(houses.shape)
Ушел: (15234, 17)
Библиотека Bokeh предоставляет интерактивные геокарты высокого разрешения, что делает ее идеальным инструментом для проверки и визуализации наших геоданных.
Bokeh – это интерактивная библиотека визуализации, предназначенная для презентаций в современных веб-браузерах. Его цель — обеспечить элегантную и лаконичную конструкцию универсальной графики и расширить эту возможность за счет высокопроизводительной интерактивности для очень больших или потоковых наборов данных. Боке может помочь всем, кто хотел бы быстро и легко создавать интерактивные графики, информационные панели и приложения для работы с данными.
import json
from bokeh.io import output_notebook, show, output_file
from bokeh.plotting import figure
from bokeh.models import GeoJSONDataSource, LinearColorMapper, ColorBar
#Read data to json.
t_geodata_json = json.loads(t_geodata.to_json())
#Convert to String like object.
json_data = json.dumps(t_geodata_json)
#Input GeoJSON source that contains features for plotting.
geosource = GeoJSONDataSource(geojson = json_data)
p = figure(title = 'GTA neighbourhoods', plot_height = 600 , plot_width = 950)
p.xgrid.grid_line_color = None
p.ygrid.grid_line_color = None
#Add patch renderer to figure.
p.patches('xs','ys', source = geosource,
line_color = 'white', line_width = 0.25, fill_alpha = 1)
#Display figure in Jupyter Notebook.
output_notebook()
#Display figure.
show(p)
Третий набор данных, представляющий данные о районе, оказался более подробным, чем ожидалось. Образование, этническая принадлежность, иммиграция/гражданство, доход, население и т. д. Потребуются только данные о доходах и коды городов/имена.
#get Toronto neighbourhood data and extract district codes with their average incomes
t_data = pd.read_csv('datasets/opendata_neighbourhoods/neighbourhood-profiles-2016-csv.csv', )
# Potential Income prospects, ids: 1030, 2273, 2274
t_income = t_data[(t_data.Characteristic == 'Neighbourhood Number') | (t_data.Characteristic == 'Total income: Average amount ($)')]
t_income = t_income.T.reset_index()
t_income.columns = ['area', 'area_code', 'average_income']
t_income = t_income.iloc[6:,:]
print(t_income.head())
print('-' * 60)
#Looks like we need to convert income and area code columns into integers
print(t_income.info())
print('-' * 60)
t_income['average_income'] = t_income['average_income'].apply(lambda a: a.replace(',','')).astype('int64')
t_income['area_code'] = t_income['area_code'].astype('int64')
print(t_income.info())
print('-' * 60)
# Calculate income mean
mean_t_income = t_income['average_income'].mean()
print('GTA mean income: ',mean_t_income)
# Merge districts geodata and Toronto income
all_t_data = t_geodata.merge(t_income, left_on='AREA_SHORT_CODE', right_on='area_code')
useless_columns_all_t_data = ['_id', 'AREA_ID', 'AREA_ATTR_ID', 'PARENT_AREA_ID', 'AREA_SHORT_CODE', 'AREA_LONG_CODE',
'AREA_NAME', 'AREA_DESC', 'X', 'Y']
all_t_data.drop(useless_columns_all_t_data, axis=1, inplace=True)
all_t_data.head()
Для каждого дома в наборе данных требуется присвоенный ему средний доход и имя района. Для этого каждый дом перепроверяется путем проецирования на карту.
Некоторые области, которые никуда не принадлежат, также лежат за пределами нашей целевой карты. После тщательных размышлений я решил, что эти области выходят за рамки проекта и могут быть безопасно удалены.
from shapely.geometry import Point, Polygon
def find_district(point, df):
'''
Returns District name, Area code, and District average income
'''
in_shape = []
for sh in df.geometry:
within = point.within(sh)
in_shape.append(within)
return df.loc[in_shape.index(True), 'area'], df.loc[in_shape.index(True), 'area_code'], df.loc[in_shape.index(True), 'average_income']
houses['mean_district_income'] = 0
houses['district_code'] = 0
not_assigned = []
for row, value in houses.iterrows():
try:
result = find_district(Point(houses.loc[row, 'long'], houses.loc[row, 'lat']), all_t_data)
houses.loc[row, 'city_district'] = result[0]
houses.loc[row, 'mean_district_income'] = result[2]
houses.loc[row, 'district_code'] = result[1]
time.sleep(0.1)
except:
not_assigned.append(row)
Исследование данных
Наконец, давайте изучим данные.
В приведенном ниже примере мы ясно видим, что наши данные о ценах и доходах сильно искажены, что, вероятно, влияет на результаты на карте. (Данные содержат много выбросов => многие цвета даже не будут отображаться на карте). Из доступных вариантов Bokeh или Seaborne я решил изучить пакет scipy.stats, чтобы найти подходящий силовой трансформатор.
Логарифмическое преобразование — отличный метод преобразования данных с положительной асимметрией в нормальные, но здесь я в основном сосредоточился на автоматическом преобразовании мощности по имени Boxcox.
Преобразование Бокса-Кокса — это способ преобразования ненормальных зависимых переменных в нормальную форму. Нормальность является важным допущением для многих статистических методов; если ваши данные не являются нормальными, применение Box-Cox означает, что вы можете запустить большее количество тестов. Преобразование Бокса-Кокса названо в честь статистиков Джорджа Бокса и сэра Дэвида Роксби Кокса, которые совместно работали над статьей 1964 года и разработали этот метод.
Однопараметрические преобразования Бокса – Кокса определяются как:
plt.figure(1, figsize=[16,3])
plt.subplots_adjust(hspace=0.4)
sns.distplot(all_t_data['average_income'], fit=stats.norm)
plt.title('Distribution of Income in GTA')
(mu, sigma) = stats.norm.fit(all_t_data['average_income'])
plt.legend(['Normal dist. ($\mu=$ {:.2f} and $\sigma=$ {:.2f} )'.format(mu, sigma)],
loc='best')
plt.figure(2, figsize=[16,3])
ax1=plt.subplot(1,2,1)
x = all_t_data['average_income']
stats.probplot(x, dist=stats.norm, plot=ax1)
ax1.set_xlabel('')
ax1.set_title('Probplot against normal distribution')
ax2=plt.subplot(1,2,2)
boxcox_lambda = stats.boxcox_normmax(x, brack=(-1.9, 2.0), method='mle')
transformed_array = stats.boxcox(x, boxcox_lambda)
# To transform back
# np.exp(np.log(boxcox_lambda*transformed_array+1)/boxcox_lambda)
stats.probplot(transformed_array, dist=stats.norm , plot=ax2)
ax2.set_title('Probplot after Box-Cox transformation')
all_t_data['average_income_transformed'] = transformed_array
plt.figure(3, figsize=[16,3])
sns.distplot(all_t_data['average_income_transformed'], fit=stats.norm)
plt.title('Distribution of Transformed Income in GTA')
(mu, sigma) = stats.norm.fit(all_t_data['average_income_transformed'])
plt.legend(['Normal dist. ($\mu=$ {:.2f} and $\sigma=$ {:.2f} )'.format(mu, sigma)],
loc='best')
plt.figure(4, figsize=[16,3])
sns.distplot(houses['final_price'], fit=stats.norm)
plt.title('Distribution of Housing Prices in GTA')
(mu, sigma) = stats.norm.fit(houses['final_price'])
plt.legend(['Normal dist. ($\mu=$ {:.2f} and $\sigma=$ {:.2f} )'.format(mu, sigma)],
loc='best')
plt.figure(5, figsize=[16,3])
ax3=plt.subplot(1,2,1)
x = houses['final_price']
stats.probplot(x, dist=stats.norm, plot=ax3)
ax3.set_xlabel('')
ax3.set_title('Probplot against normal distribution')
ax4=plt.subplot(1,2,2)
boxcox_lambda = stats.boxcox_normmax(x, brack=(-1.9, 2.0), method='mle')
transformed_array = stats.boxcox(x, boxcox_lambda)
stats.probplot(transformed_array, dist=stats.norm , plot=ax4)
ax4.set_title('Probplot after Box-Cox transformation')
houses['final_price_transformed'] = transformed_array
plt.figure(6, figsize=[16,3])
sns.distplot(houses['final_price_transformed'], fit=stats.norm)
plt.title('Distribution of Transformed Housing Prices in GTA')
(mu, sigma) = stats.norm.fit(houses['final_price_transformed'])
plt.legend(['Normal dist. ($\mu=$ {:.2f} and $\sigma=$ {:.2f} )'.format(mu, sigma)],
loc='best')
all_t_data['average_income_log']= np.log(all_t_data['average_income'])
houses['final_price_log']= np.log(houses['final_price'])
plt.show()
Имея нормализованные данные, мы готовы нанести их на график и проанализировать, как они соотносятся друг с другом.
from bokeh.palettes import brewer
from bokeh.models import Label, HoverTool, AjaxDataSource
#Read data to json.
t_geodata_json = json.loads(all_t_data.to_json())
#Convert to String like object.
json_data = json.dumps(t_geodata_json)
#Input GeoJSON source that contains features for plotting.
geosource = GeoJSONDataSource(geojson = json_data)
#Add hover tool
hover = HoverTool(tooltips = [ ('District','@area'),('Average Income', '@average_income')])
p = figure(title = 'GTA Neighbourhood Income', plot_height = 600 , plot_width = 950, tools = [hover])
p.xgrid.grid_line_color = None
p.ygrid.grid_line_color = None
# Add patch renderer to figure.
palette = brewer['Spectral'][11] #[::-1]
color_mapper = LinearColorMapper(palette = palette,
low = all_t_data['average_income_transformed'].min(),
high = all_t_data['average_income_transformed'].max())
p.patches('xs','ys', source = geosource,
line_color = 'white', line_width = 0.25, fill_alpha = 1,
fill_color = {'field' :'average_income_transformed', 'transform' : color_mapper})
#Display figure in Jupyter Notebook.
output_notebook()
#Display figure.
show(p)
Нанесение данных о жилье на карту.
Вы можете увидеть два уровня:
- Нижний уровень представляет районы и среднегодовую заработную плату в районах.
- Средний уровень представляет линии и станции метро Торонто.
- Вверху вверху представлены цены на жилье
# upload TTC stations and lines
slines_geodata = gpd.read_file('datasets/transit/TTC_SubwayLines_20120119/subwayLines.shp')
sstations_geodata = gpd.read_file('datasets/transit/TTC_SubwayStns_20120119/subwayStations.shp')
slines_geodata['colors'] = ['limegreen', 'plum', 'dodgerblue', 'gold']
for dataset in [slines_geodata, sstations_geodata]:
dataset['geometry'] = dataset['geometry'].to_crs({'init': 'EPSG:4269'})
slines_geodata
from bokeh.models.glyphs import MultiLine
from bokeh.models import Cross,NumeralTickFormatter
# Input GeoJSON source that contains features for plotting.
geosource = GeoJSONDataSource(geojson = all_t_data.to_json())
slines_source = GeoJSONDataSource(geojson = slines_geodata.to_json())
sstation_source = GeoJSONDataSource(geojson = sstations_geodata.to_json())
p = figure(title = 'GTA housing prices projected on average district wages', plot_height = 680 , plot_width = 950,
active_scroll='wheel_zoom')
p.xgrid.grid_line_color = None
p.ygrid.grid_line_color = None
# Add patch renderer to figure - Areas and Average Income.
palette_1 = brewer['RdBu'][11] #[::-1]
color_mapper_1 = LinearColorMapper(palette = palette_1,
low = all_t_data['average_income_log'].min(),
high = all_t_data['average_income_log'].max())
p.patches('xs','ys', source = geosource,
line_color = 'white', line_width = 0.25, fill_alpha = 0.7,
fill_color = {'field' :'average_income_log', 'transform' : color_mapper_1})
#Add subway lines
p.multi_line('xs', 'ys', source=slines_source, color='colors', line_width=2.5)
#Define custom tick labels for color bar.
tick_labels = {'10': '22K','10.5':'35K','11':'60K','11.5':'100k', '12':'162K', '12.5':'268k'}
#Create color bar 1.
color_bar_1 = ColorBar(color_mapper=color_mapper_1, label_standoff=8,width = 500, height = 15,
border_line_color=None,location = (0,0), orientation = 'horizontal', major_label_overrides = tick_labels)
#Specify figure layout.
p.add_layout(color_bar_1, 'below')
# Add houses to the map
palette = brewer['Spectral'][11] #[::-1]
color_mapper = LinearColorMapper(palette = palette,
low = houses['final_price_log'].min(),
high = houses['final_price_log'].max())
geosource_2 = AjaxDataSource(houses)
p1 = p.circle('long', 'lat', alpha=0.3, radius=0.0006, color={'field': 'final_price_log', 'transform': color_mapper}, source=geosource_2)
#Define custom tick labels for color bar.
tick_labels = {'12': '160K','13':'440K','14':'1,200K','15':'3,270k', '16':'8,890K'}
#Create color bar 2.
color_bar = ColorBar(color_mapper=color_mapper, label_standoff=8,width = 500, height = 15,
border_line_color=None,location = (0,0), orientation = 'horizontal', major_label_overrides = tick_labels)
#Add hover tool
hover = HoverTool(renderers=[p1],
tooltips =[('District','@city_district'),
('District Income', '$ @mean_district_income{0,0}'),
('Price','$ @final_price{0,0}'),
('Type', '@type'),
('MLS','@mls')])
p.add_tools(hover)
#Specify figure layout.
p.add_layout(color_bar, 'below')
# Add subway stations
p.circle('x', 'y', source=sstation_source, radius=0.001, color="white", line_color='black', line_width=0.3)
#Display figure in Jupyter Notebook.
output_notebook()
#Display figure.
show(p)
Основываясь на сюжете, необходимо сделать несколько важных замечаний:
- Средний доход и конечная цена коррелируют друг с другом
- Большинство квартир в кондо (в основном зеленые) сосредоточены в центре города и за линиями метро в направлении Северного Йорка. (от 500к до 700к)
- Большинство областей, отмеченных желтым/красным цветом, ограничены отдельными домами. Для которых цена варьируется от 800k до 1500k
- Логично предположить, что в центре города было продано непропорционально больше квартир, чем в любом другом районе.
- В большинстве районов с низким доходом, таких как Скарборо, есть кооперативные квартиры и кондоминиумы в диапазоне от 200 до 400 тысяч. Три самых доходных / самых дорогих объекта: York Mills, Rosedale-Moore Park, Forest Hill Sout.
Дальнейший анализ и преобразования
Цели:
- Функция разделения спальни на уровни выше (bedrooms_ag) и ниже уровня (bedrooms_bg). Спальни выше класса могут быть более ценными, чем ниже уровня
- Установить максимальное ограничение для парковки
- Преобразовать квадратные метры в целые числа
- Объединить типы второстепенных групп
# Before we make any transformations to our data let's make a copy
houses_edited = houses.copy()
def convert_bedrooms(n_plus_n, above = True):
removed_beds_str = n_plus_n.replace(' beds', '')
if above == True:
return int(removed_beds_str[0])
else:
if len(removed_beds_str) > 2:
return int(removed_beds_str[-1])
else: return 0
def convert_parking(parking):
parking = int(parking.replace(' parking', '').replace('no', '0'))
if parking > 10:
return 11
else: return parking
def convert_sqft(value):
value = value.replace(' sq. ft.', '').replace('0-499', '0–499')
value = value.split('–')
try: return (int(value[0]) + int(value[1]))/2
except: 0
houses_edited['bedrooms_ag'] = houses_edited.bedrooms.apply(convert_bedrooms)
houses_edited['bedrooms_bg'] = houses_edited.bedrooms.apply(convert_bedrooms, above = False)
houses_edited.bathrooms = houses_edited.bathrooms.apply(lambda a: int(a.replace(' baths', '')))
houses_edited.parking = houses_edited.parking.apply(convert_parking)
houses_edited.sqft = houses_edited.sqft.apply(convert_sqft).round(0)
houses_edited.type.replace({
'Duplex': 'Plex',
'Triplex': 'Plex',
'Fourplex': 'Plex',
'Multiplex': 'Plex',
'Det Condo': 'Condo Apt',
'Leasehold Condo': 'Condo Apt'
}, inplace=True)
houses_edited.to_csv('houses_edited.csv')
houses_edited.head(2)
Создайте фиктивные столбцы для категориальных переменных
houses_edited = pd.read_csv('houses_edited.csv', index_col='index')
needed_columns = ['sqft', 'parking' , 'mean_district_income', 'bedrooms_ag', 'bedrooms_bg', 'type', 'final_price']
# houses_edited.info()
houses_dummies = pd.get_dummies(houses_edited[needed_columns])
predictors = houses_dummies.drop('final_price', axis=1).columns
houses_dummies.head(5)
Недавно я наткнулся на еще один замечательный инструмент под названием Профилирование Panda.
Он генерирует отчеты о профилях из панд DataFrame. Функция pandas df.describe() великолепна, но немного проста для серьезного исследовательского анализа данных. pandas_profiling расширяет DataFrame pandas с помощью df.profile_report() для быстрого анализа данных.
Для каждого столбца в интерактивном HTML-отчете представлены следующие статистические данные, если они относятся к типу столбца:
- Основные сведения: тип, уникальные значения, отсутствующие значения.
- Квантильные статистические данные, такие как минимальное значение, Q1, медиана, Q3, максимум, диапазон, межквартильный диапазон
- Описательные статистические данные, такие как среднее значение, мода, стандартное отклонение, сумма, среднее абсолютное отклонение, коэффициент вариации, эксцесс, асимметрия.
- Наиболее часто встречающиеся значения
- Гистограмма
- Корреляции выделение переменных с высокой степенью корреляции, матриц Спирмена, Пирсона и Кендалла
- Матрица отсутствующих значений, количество, тепловая карта и дендрограмма отсутствующих значений
Чтобы определить, какие функции наиболее важны для нашего анализа, а какие значения требуют условного исчисления, нам нужно построить матрицу важности функций и построить ее с помощью Matplotlib.
Существует множество методов оценки важности функций, здесь мы воспользуемся простым: Tree Repressor.
from sklearn.ensemble import ExtraTreesRegressor
sns.set_style("whitegrid")
plt.figure(figsize=[13,5])
df = houses_dummies.dropna(subset=['sqft'])
model = ExtraTreesRegressor(n_estimators=10)
model.fit(df[predictors], df.final_price)
# display the relative importance of each attribute
sorted_feature_importance = sorted(zip(model.feature_importances_, df[predictors].columns))
x = [a[0] for a in sorted_feature_importance]
y = [a[1] for a in sorted_feature_importance]
plt.scatter(x,y)
plt.title('Feature Importances')
plt.show()
Похоже, что одной из самых ценных функций является sqft. Поскольку нам не хватает некоторых из этих данных, я создал подсчет значений для существующих и отсутствующих наблюдений объектов sqft.
house_type_notna, house_type_na = [], []
houses_types = pd.unique(houses_edited.type)
for house_type in houses_types:
house_type_notna.append(sum((houses_edited.type == house_type) & (pd.notna(houses_edited.sqft))))
house_type_na.append(sum((houses_edited.type == house_type) & (pd.isna(houses_edited.sqft))))
ind = np.arange(houses_types.shape[0])
plt.figure(figsize=[15,5])
width = 0.5
p1 = plt.barh(ind, house_type_notna, color='lightgreen')
p2 = plt.barh(ind, house_type_na, left=house_type_notna, color='lightcoral')
for no in ind:
percent_missing = np.round((house_type_na[no] / (house_type_notna[no] + house_type_na[no]))*100, 1)
plt.text(house_type_notna[no] + house_type_na[no], ind[no], str(percent_missing) + '%')
plt.xlim([0,8000])
plt.legend(['Present sqft', 'Missing sqft'])
plt.yticks(ind, houses_types)
plt.title('Sqft value counts in terms of missing and present data for all property types', fontsize=13)
plt.show()
К сожалению, большинство частных домов не имеют площади в квадратных футах, поэтому вменение этих данных имеет приоритет.
# Create features Heatmap to see how feature are correlated with final_price
def correlantion_heatmap(df, center = None):
fig, ax = plt.subplots(figsize = [14,12])
colormap = sns.diverging_palette(10, 220, as_cmap=True)
fig = sns.heatmap(df.corr(),
cmap = colormap,
center = center,
annot = True,
linewidths = 0.1,
annot_kws={'fontsize':9})
correlantion_heatmap(houses_dummies, center=0)
sns.set_style("whitegrid")
plt.subplots(figsize=[16,8])
sorted_medians = sorted([(houses_edited[houses_edited.type==st].final_price.median(), st) for st in houses_types], reverse=True)
sns.boxplot(y = 'type', x = 'final_price', data = houses_edited, orient='h', showfliers=False, meanline =True, showmeans=True,
order=[b for a,b in sorted_medians])
plt.title('Sorted distribution of Final Price among all property types', fontsize=16)
plt.show()
plt.subplots(figsize=[16,8])
sorted_medians = sorted([(houses_edited[houses_edited.type==st].sqft.median(), st) for st in houses_types], reverse=True)
sns.boxplot(y = 'type', x = 'sqft', data = houses_edited, orient='h', showfliers=True, meanline =True, showmeans=True,
order=[b for a,b in sorted_medians])
plt.title('Sorted distribution of Sqft among all property types', fontsize=16)
plt.show()
# Plot pairwise relationships
sns.set_style("white")
pp = sns.pairplot(houses_dummies[pd.notna(houses_dummies.sqft)][predictors[:5]], diag_kind = 'hist', plot_kws=dict(s=14, alpha=.2, linewidth=0))
pp.set(xticklabels=[])
plt.show()
# Check how final_price and sqft are related to each other among all property types
sns.set_style("whitegrid")
plt.figure(figsize=[15,8])
sns.scatterplot(x='sqft', y='final_price' ,data=houses_edited,
hue = houses_edited.type.replace({
'Plex': 'Other',
'Link': 'Other',
'Co-Ownership Apt': 'Other',
'Co-Op Apt': 'Other',
'Store W/Apt/Offc': 'Other'}),
linewidth=0, alpha=0.3, palette='bright', s=30)
plt.title('Final price and Sqft relationship among all property types', fontsize=16)
plt.show()
Чтобы заменить отсутствующие значения, мы создадим собственный импутер. Чтобы обеспечить его совместимость с конвейером, мы унаследуем такие базовые классы, как BaseEstimator и TransformerMixin.
Рассмотрены вменения:
- Случайный лес. Вмените пропущенные значения в данных предиктора, используя близость от randomForest
- Кнн. Используйте непараметрический алгоритм, называемый k-ближайшими соседями (KNN), для замены пропущенных значений.
- Мыши. Используйте многомерное вменение с помощью цепных уравнений
- Среднее. Замените отсутствующие значения на Среднее значение соответствующего свойства time
- Медиана. Замените отсутствующие значения медианой соответствующего времени свойства
- Режим. Замените отсутствующие значения режимом соответствующего времени свойства.
После бесчисленных сравнений перекрестной проверки различных алгоритмов я выбрал MICE из-за его стабильности.
from sklearn.base import BaseEstimator, TransformerMixin
class Impute_sqft(BaseEstimator, TransformerMixin):
def __init__(self, how='random forest'):
self.how = how
def fit(self, X, y=None):
return self
def transform(self, X):
import sys
sys.setrecursionlimit(100000) #Increase the recursion limit of the OS
from impyute.imputation.cs import fast_knn, mice
result = X.copy()
if self.how == 'random forest':
train_X = houses_dummies.dropna(subset=['sqft']).drop(columns=['sqft'])
train_Y = houses_dummies.dropna(subset=['sqft'])['sqft']
test_X = houses_dummies[pd.isna(houses_dummies.sqft)].drop(columns=['sqft'])
rf = RandomForestRegressor()
rf.fit(train_X, train_Y)
pred_Y = rf.predict(test_X)
result.loc[test_X.index,'sqft'] = pred_Y
if self.how == 'knn':
knn_n = 30
result = fast_knn(houses_dummies, k=knn_n)
result.columns = houses_dummies.columns
result.index = houses_dummies.index
result = result.loc[X.index,:]
if self.how == 'mice':
result = mice(houses_dummies)
result.columns = houses_dummies.columns
result.index = houses_dummies.index
result = result.loc[X.index,:]
if self.how == 'mean':
result['sqft'] = houses_edited.groupby('type')['sqft'].transform(lambda x: x.fillna(x.mean(skipna=True)))
if self.how == 'median':
result['sqft'] = houses_edited.groupby('type')['sqft'].transform(lambda x: x.fillna(x.median(skipna=True)))
if self.how == 'mode':
result['sqft'] = houses_edited.groupby('type')['sqft'].transform(lambda x: x.fillna(x.mode()[0]))
return result[predictors]
imputer = Impute_sqft(how='mice')
houses_nona = imputer.fit_transform(houses_dummies[predictors])
Ранее мы уже преобразовывали данные с помощью Boxcox, поэтому в этом следующем шаге нет необходимости. Одна проблема с Boxcox
Неожиданная проблема, с которой я столкнулся при работе с Boxcox, заключалась в удобном возврате преобразованных значений. Он не может преобразовывать отрицательные значения, что явно конфликтует при попытке обработать широкий спектр моделей, которые дают как положительные, так и отрицательные прогнозы. По этой причине вместо этого было рассмотрено преобразование Йео-Джонсона.
Преобразование Йео-Джонсона можно рассматривать как расширение преобразования Бокса-Кокса. Оно обрабатывает как положительные, так и отрицательные значения, тогда как преобразование Бокса-Кокса обрабатывает только положительные значения. И то, и другое можно использовать для преобразования данных, чтобы улучшить нормальность.
Опять же, давайте проверим асимметрию признаков и исправим самые верхние из них.
numeric_feats = ['sqft', 'mean_district_income', 'final_price']
skewness = pd.DataFrame({
'Skew': houses_dummies[numeric_feats].apply(lambda x: stats.skew(x.dropna())).sort_values(ascending=False)})
skewness = skewness[abs(skewness.Skew) > 0.75]
skewness
from sklearn.preprocessing import PowerTransformer
from sklearn.metrics import mean_squared_error, mean_absolute_error
houses_nona_copy = houses_nona.copy()
houses_nona_copy['mean_district_income'] = houses_nona.mean_district_income/10
houses_nona_copy['final_price'] = houses_dummies.final_price/100
# np.exp(np.log(boxcox_lambda*transformed_array+1)/boxcox_lambda)
yeo_johnson_lambdas = []
# for feature in skewness.index:
# boxcox_lambda = stats.boxcox_normmax(houses_nona_copy[feature], brack=(-1.9, 2.0), method='mle')
# houses_nona_copy[feature] = boxcox1p(houses_nona_copy[feature], boxcox_lambda)
# boxcox_lambdas.append(boxcox_lambda)
# for feature in skewness.index:
# yeo_lambda = stats.yeojohnson_normmax(houses_nona_copy[feature].astype('float64'), brack=[-2.9,20])
# houses_nona_copy[feature] = stats.yeojohnson(houses_nona_copy[feature].astype('float64'), lmbda=yeo_lambda)
# boxcox_lambdas.append(yeo_lambda)
pt = PowerTransformer(method='yeo-johnson')
for numeric_feat in numeric_feats:
houses_nona_copy[numeric_feat] = pt.fit_transform(houses_nona_copy[[numeric_feat]])
yeo_johnson_lambdas.append(pt.lambdas_[0])
# rf = RandomForestRegressor()
# rf.fit(houses_nona_copy[predictors], houses_nona_copy['final_price'])
# prediction = rf.predict(houses_nona_copy[predictors])
# y = pt.inverse_transform(houses_nona_copy[['final_price']]) * 100
# y_pred = pt.inverse_transform(pd.DataFrame(prediction)) * 100
# mean_absolute_error(y, y_pred)
Разделите массивы или матрицы на случайные поезда и тестовые подмножества
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(houses_nona_copy[predictors], houses_nona_copy.final_price, train_size=0.7,
stratify=houses_nona_copy['type_Condo Apt'].values)
Изучите модели и отметьте лучшие из них
Создавайте собственные функции оценки. Во-первых, преобразуйте как прогнозируемые, так и фактические значения в доллары и оцените одну из этих четырех метрик:
def custom_mae_(y, y_pred, **kwargs):
y_pred[y_pred > 5] = 5
y = pt.inverse_transform(pd.DataFrame(y)) * 100
y_pred = pt.inverse_transform(pd.DataFrame(y_pred)) * 100
return mean_absolute_error(y, y_pred)
custom_mae = make_scorer(custom_mae_, greater_is_better=False)
def custom_mape_(y, y_pred, **kwargs):
y_pred[y_pred > 5] = 5
y = pt.inverse_transform(pd.DataFrame(y)) * 100
y_pred = pt.inverse_transform(pd.DataFrame(y_pred)) * 100
return np.mean(np.abs((y - y_pred) / y)) * 100
custom_mape = make_scorer(custom_mape_, greater_is_better=False)
def custom_rmse_(y, y_pred, **kwargs):
y_pred[y_pred > 5] = 5
y = pt.inverse_transform(pd.DataFrame(y)) * 100
y_pred = pt.inverse_transform(pd.DataFrame(y_pred)) * 100
return np.sqrt(mean_squared_error(y, y_pred))
custom_rmse = make_scorer(custom_rmse_, greater_is_better=False)
Теперь мы можем создать пространство для оценки и предварительно выбрать 20 моделей, которые могут помочь в нашем решении. Затем мы пропустим каждый алгоритм через наш набор поездов, проверим результаты и сохраним все в один фрейм данных.
from sklearn.ensemble import AdaBoostRegressor, BaggingRegressor, ExtraTreesRegressor, GradientBoostingRegressor, RandomForestRegressor
from sklearn.kernel_ridge import KernelRidge
from sklearn.linear_model import LinearRegression, LassoCV, RidgeCV, ElasticNetCV, ARDRegression, BayesianRidge, HuberRegressor, LarsCV, OrthogonalMatchingPursuitCV
from sklearn.neighbors import KNeighborsRegressor
from sklearn.svm import NuSVR, SVR
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.tree import DecisionTreeRegressor
from sklearn.model_selection import ShuffleSplit, cross_validate, GridSearchCV, KFold
from sklearn.isotonic import IsotonicRegression
from sklearn.cross_decomposition import PLSRegression
from sklearn.preprocessing import RobustScaler
from sklearn.pipeline import make_pipeline
from sklearn.metrics import make_scorer
import lightgbm as lgb
import xgboost as xgb
#Machine Learning Algorithm (MLA) Selection and Initialization
MLA = [
# #Ensemble Methods
AdaBoostRegressor(),
BaggingRegressor(),
ExtraTreesRegressor(),
GradientBoostingRegressor(),
RandomForestRegressor(),
lgb.LGBMRegressor(),
xgb.XGBRegressor(objective='reg:squarederror'),
# #Kernel ridge regression
# KernelRidge(),
# #GLM
LinearRegression(),
make_pipeline(RobustScaler(), LassoCV()),
make_pipeline(RobustScaler(), ElasticNetCV()),
RidgeCV(),
LarsCV(),
BayesianRidge(),
HuberRegressor(),
OrthogonalMatchingPursuitCV(),
# #Decomposition
PLSRegression(),
# #KNN
KNeighborsRegressor(),
#SVM
NuSVR(),
SVR(),
# #Tree Models
DecisionTreeRegressor()
]
#split dataset in cross-validation with this splitter class: http://scikit-learn.org/stable/modules/generated/sklearn.model_selection.ShuffleSplit.html#sklearn.model_selection.ShuffleSplit
#note: this is an alternative to train_test_split
cv_split = KFold(n_splits=4, shuffle=True, random_state=10) # run model 10x with 60/30 split intentionally leaving out 10%
#create table to compare MLA metrics
MLA_columns = ['MLA Name', 'MLA Parameters','MLA Train Accuracy Mean', 'MLA Test Accuracy Mean', 'MLA Test Accuracy 3*STD' ,'MLA Time']
MLA_compare = pd.DataFrame(columns = MLA_columns)
#index through MLA and save performance to table
row_index = 0
for alg in MLA:
#set name and parameters
MLA_name = alg.__class__.__name__
print(MLA_name + '_'*10, end="\r")
MLA_compare.loc[row_index, 'MLA Name'] = MLA_name
MLA_compare.loc[row_index, 'MLA Parameters'] = str(alg.get_params())
#score model with cross validation: http://scikit-learn.org/stable/modules/generated/sklearn.model_selection.cross_validate.html#sklearn.model_selection.cross_validate
cv_results = cross_validate(alg, X_train, y_train, cv = cv_split,
return_train_score=True, scoring=custom_mae)
MLA_compare.loc[row_index, 'MLA Time'] = cv_results['fit_time'].mean()
MLA_compare.loc[row_index, 'MLA Train Accuracy Mean'] = -cv_results['train_score'].mean()
MLA_compare.loc[row_index, 'MLA Test Accuracy Mean'] = -cv_results['test_score'].mean()
#if this is a non-bias random sample, then +/-3 standard deviations (std) from the mean, should statistically capture 99.7% of the subsets
MLA_compare.loc[row_index, 'MLA Test Accuracy 3*STD'] = cv_results['test_score'].std()*3 #let's know the worst that can happen!
row_index+=1
MLA_compare.sort_values('MLA Test Accuracy Mean', inplace = True, ascending = True)
MLA_compare
fig, ax = plt.subplots(figsize = [16,8])
sns.barplot(x = 'MLA Test Accuracy Mean', y = 'MLA Name', data = MLA_compare,)
plt.title('ML Algorithm Test Accuracy Score', fontsize=13)
plt.xlabel('Mean Error')
plt.ylabel('Algorithm')
plt.show()
Краткий список лучших оценщиков (также моих любимых):
- СветGBM
- XGBoost
- RandomForest
Настройка гиперпараметров
Структура этого абзаца:
- LightGBM с Hyperopt
- LightGBM с RandomSearchCV
- Hyperopt против RandomSearchCV
- XGboost с Hyperopt
- RandomForest с Hyperopt
LightGBM с Hyperopt
Существует множество подходов к настройке гиперпараметров. Некоторые примеры включают ручные пробы и ошибки, выполнение поиска по сетке/случайного поиска по сетке и более сложные автоматизированные подходы, такие как байесовская оптимизация (которые предоставляются Hyperopt). сильный> пакет).
Hyperopt – это пакет глобальной оптимизации с ограничениями, основанный на байесовском выводе и гауссовском процессе, который пытается найти максимальное значение неизвестной функции за минимально возможное количество итераций. Этот метод особенно подходит для оптимизации высокозатратных функций, когда важен баланс между разведкой и эксплуатацией.
В задаче байесовской оптимизации есть четыре части:
- Целевая функция: то, что мы хотим свести к минимуму, в данном случае ошибка проверки модели машинного обучения по отношению к гиперпараметрам.
- Доменное пространство: значения гиперпараметров для поиска.
- Алгоритм оптимизации: метод построения суррогатной модели и выбора следующих значений гиперпараметров для оценки.
- История результатов: сохраненные результаты оценок целевой функции, состоящей из гиперпараметров и потерь при проверке.
Начнем с определения показателей оценки и создания цели.
import csv
from hyperopt import STATUS_OK
from timeit import default_timer as timer
MAX_EVALS = 200
train_set = lgb.Dataset(X_train, label = y_train)
# f(preds: array, train_data: Dataset) -> name: string, eval_result: float, is_higher_better: bool
# How: https://github.com/Microsoft/LightGBM/blob/master/examples/python-guide/advanced_example.py
def custom_mae_lgb(y_pred, y):
# Replace too big values
y_pred[y_pred > 5] = 5
y = pt.inverse_transform(pd.DataFrame(y.get_label())) * 100
y_pred = pt.inverse_transform(pd.DataFrame(y_pred)) * 100
return 'custom_mae_lgb', mean_absolute_error(y, y_pred), False
def objective(params):
"""Objective function for Gradient Boosting Machine Hyperparameter Optimization"""
# Keep track of evals
global ITERATION
ITERATION += 1
# Retrieve the subsample if present otherwise set to 1.0
subsample = params['boosting_type'].get('subsample', 1.0)
# Extract the boosting type
params['boosting_type'] = params['boosting_type']['boosting_type']
params['subsample'] = subsample
# Make sure parameters that need to be integers are integers
for parameter_name in ['num_leaves', 'subsample_for_bin', 'min_child_samples']:
params[parameter_name] = int(params[parameter_name])
start = timer()
# Perform n_folds cross validation
cv_results = lgb.cv(params, train_set, num_boost_round = 1000, folds=cv_split,
early_stopping_rounds = 50, seed = 50, stratified=False, feval=custom_mae_lgb)
run_time = timer() - start
# Extract the best score, Loss must be minimized
best_score = np.min(cv_results['custom_mae_lgb-mean'])
loss = best_score
# Boosting rounds that returned the highest cv score
n_estimators = int(np.argmin(cv_results['custom_mae_lgb-mean']) + 1)
# Dictionary with information for evaluation
return {'loss': best_score, 'params': params, 'iteration': ITERATION,
'estimators': n_estimators,
'train_time': run_time, 'status': STATUS_OK}
Определить пространство поиска
from hyperopt import hp, tpe, Trials, fmin
from hyperopt.pyll.stochastic import sample
space = {
'class_weight': hp.choice('class_weight', [None, 'balanced']),
'boosting_type': hp.choice('boosting_type', [{'boosting_type': 'gbdt', 'subsample': hp.uniform('gdbt_subsample', 0.5, 1)},
{'boosting_type': 'dart', 'subsample': hp.uniform('dart_subsample', 0.5, 1)},
{'boosting_type': 'goss', 'subsample': 1.0}]),
'num_leaves': hp.quniform('num_leaves', 30, 150, 1),
'learning_rate': hp.loguniform('learning_rate', np.log(0.01), np.log(0.2)),
'subsample_for_bin': hp.quniform('subsample_for_bin', 20000, 300000, 20000),
'min_child_samples': hp.quniform('min_child_samples', 20, 500, 5),
'reg_alpha': hp.uniform('reg_alpha', 0.0, 1.0),
'reg_lambda': hp.uniform('reg_lambda', 0.0, 1.0),
'colsample_bytree': hp.uniform('colsample_by_tree', 0.6, 1.0)
}
# optimization algorithm
tpe_algorithm = tpe.suggest
# Keep track of results
bayes_trials = Trials()
# Global variable
global ITERATION
ITERATION = 0
%sc
# Run optimization
lgb_best = fmin(fn = objective, space = space, algo = tpe.suggest,
max_evals = MAX_EVALS, trials = bayes_trials, rstate = np.random.RandomState(10))
# Sort the trials with lowest loss (highest AUC) first
lgb_trials_results = sorted(bayes_trials.results, key = lambda x: x['loss'])
Преобразовать в фрейм данных
lgb_results = pd.DataFrame(lgb_trials_results) lgb_results.head()
# Save best parameters best_lgb_estimators= lgb_results.loc[0,'estimators'] best_lgb_params= lgb_results.loc[0,'params'] #Fit the best parameters to the model best_lgb_model = lgb.LGBMRegressor(n_estimators=best_lgb_estimators, n_jobs = -1, random_state = 50, **best_lgb_params) #Save parameters. *Just in case* # best_lgb_model = lgb.LGBMRegressor(boosting_type='gbdt', class_weight='balanced', # colsample_bytree=0.9414862936943954, importance_type='split', # learning_rate=0.07074010161873992, max_depth=-1, # min_child_samples=20, min_child_weight=0.001, min_split_gain=0.0, # n_estimators=976, n_jobs=-1, num_leaves=38, objective=None, # random_state=50, reg_alpha=0.6547434726382414, # reg_lambda=0.7404502093528875, silent=True, # subsample=0.5213499245813108, subsample_for_bin=180000, # subsample_freq=0)
LightGBM со случайным поиском
# Hyperparameter grid
param_grid = {
'class_weight': [None, 'balanced'],
'boosting_type': ['gbdt', 'goss', 'dart'],
'num_leaves': list(range(30, 150)),
'learning_rate': list(np.logspace(np.log(0.005), np.log(0.2), base = np.exp(1), num = 1000)),
'subsample_for_bin': list(range(20000, 300000, 20000)),
'min_child_samples': list(range(20, 500, 5)),
'reg_alpha': list(np.linspace(0, 1)),
'reg_lambda': list(np.linspace(0, 1)),
'colsample_bytree': list(np.linspace(0.6, 1, 10))
}
# Subsampling (only applicable with 'goss')
subsample_dist = list(np.linspace(0.5, 1, 100))
params = {key: np.random.choice(value) for key, value in param_grid.items()}
params
def random_objective(params, iteration, n_folds = N_FOLDS):
"""Random search objective function. Takes in hyperparameters
and returns a list of results to be saved."""
start = timer()
# Perform n_folds cross validation
cv_results = lgb.cv(params, train_set, num_boost_round=1000, folds=cv_split,
early_stopping_rounds=50, seed=50, stratified=False, feval=custom_mae_lgb)
end = timer()
best_score = np.min(cv_results['custom_mae_lgb-mean'])
# Loss must be minimized
loss = best_score
# Boosting rounds that returned the highest cv score
n_estimators = int(np.argmin(cv_results['custom_mae_lgb-mean']) + 1)
# Return list of results
return [loss, params, iteration, n_estimators, end - start]
np.random.seed(50)
# Dataframe to hold cv results
random_results = pd.DataFrame(columns = ['loss', 'params', 'iteration', 'estimators', 'time'],
index = list(range(MAX_EVALS)))
# Iterate through the specified number of evaluations
for i in range(MAX_EVALS):
print(i, end='\r')
# Randomly sample parameters for gbm
params = {key: np.random.choice(value) for key, value in param_grid.items()}
if params['boosting_type'] == 'goss':
# Cannot subsample with goss
params['subsample'] = 1.0
else:
# Subsample supported for gdbt and dart
params['subsample'] = np.random.choice(subsample_dist)
results_list = random_objective(params, i)
# Add results to next row in dataframe
random_results.loc[i, :] = results_list
# Sort results by best validation score
random_results.sort_values('loss', ascending = True, inplace = True)
random_results.reset_index(inplace = True, drop = True)
random_results.head()
Hyperopt против RandomSearchCV
Как видите, байесовская оптимизация занимает меньше времени, более эффективно ищет параметры, тем самым достигая меньших потерь после каждой итерации.
best_random_params = random_results.loc[0,'params']
best_random_estimators = random_results.loc[0,'estimators']
print('The best model through Bayesian search scores: ', int(lgb_results.loc[0,'loss']), 'found in ',
int(lgb_results.loc[0,'iteration']), ' Total time: ', int(lgb_results.train_time.sum()))
print('The best model through Random search scores: ', int(random_results.loc[0,'loss']), 'found in ',
int(random_results.loc[0,'iteration']), ' Total time: ', int(random_results.time.sum()))
Лучшая модель по результатам байесовского поиска: 95339 найдено за 164 Всего за время: 4816
Лучшая модель по результатам случайного поиска: 97548 найдено за 22 Всего за время: 6445
Сравнивая наши результаты с результатами случайного поиска, как в числах, так и в цифрах, мы можем начать понимать, как работает байесовская оптимизация. Во-первых, мы можем посмотреть на лучшие гиперпараметры (определенные по ошибке проверки) из обоих поисков.
# lgb_results = lgb_results.rename({'train_time':'time'}, axis=1).drop(columns=['status'], axis=1)
lgb_results['method'] = 'Bayesian search'
random_results['method'] = 'Random Search'
all_lgb_results = pd.concat([lgb_results, random_results])
all_lgb_results = pd.concat([all_lgb_results,
pd.DataFrame(columns=list(all_lgb_results.loc[0,'params'].to_dict()[0].keys()),index=all_lgb_results.index)],
axis=1)
all_lgb_results.reset_index(inplace=True)
# Create a new colmn for each parameter and copy-paste it from params column
for i, params in enumerate(all_lgb_results['params']):
for param in params.keys():
all_lgb_results.loc[i,param] = params[param]
all_lgb_results.head()
Похоже, что потери при проверке для байесовской оптимизации ниже, чем для случайного поиска. Однако это не обязательно приводит к лучшему результату тестирования!
Байесовская оптимизация должна улучшаться со временем. Давайте нанесем баллы на итерацию, чтобы определить, было ли какое-либо улучшение.
for param in ['iteration', 'estimators', 'num_leaves', 'subsample_for_bin', 'min_child_samples', 'loss']:
all_lgb_results[param] = all_lgb_results[param].astype('float')
lmplot_1 = sns.lmplot('iteration', 'loss', hue='method', data=all_lgb_results, aspect=2, height=7, order=3)
plt.title('Loss vs Iteration')
plt.show()
Случайный поиск использует типы бустинга с одинаковой частотой. Однако байесовская оптимизация могла бы определить, что один тип бустинга лучше другого для конкретной проблемы.
plt.figure(figsize=[15,5])
sns.countplot('boosting_type', hue='method', data=all_lgb_results)
plt.title('Boosting type count', fontsize=13)
plt.show()
Разумеется, байесовская оптимизация выбирала тип повышения дерева решений с градиентным усилением чаще, чем любой другой тип. Мы могли бы использовать эту информацию для информирования последующих поисков оптимальных гиперпараметров, сосредоточив внимание на меньших областях.
num_params = ['num_leaves', 'learning_rate', 'subsample_for_bin', 'min_child_samples',
'reg_alpha', 'reg_lambda', 'colsample_bytree']
plt.figure(figsize=[15,13])
plt.subplots_adjust(hspace=1)
for i, param in enumerate(num_params):
plt.subplot(4,2,i+1)
for method in ['Bayesian search', 'Random Search']:
sns.kdeplot(all_lgb_results[all_lgb_results.method==method][param], label=method)
plt.title(param + ' Distribution')
plt.xlabel(param), plt.ylabel('Density')
XGBoost с Hyperopt
from hyperopt.pyll.stochastic import sample
MAX_EVALS = 102
train_set = xgb.DMatrix(X_train.values, label = y_train.values)
def custom_mae_xgb(y_pred, y):
# Replace too big values
y_pred[y_pred > 5] = 5
y = pt.inverse_transform(pd.DataFrame(y.get_label())) * 100
y_pred = pt.inverse_transform(pd.DataFrame(y_pred)) * 100
return 'custom_mae_xgb', mean_absolute_error(y, y_pred)
def objective(params):
"""Objective function for Gradient Boosting Machine Hyperparameter Optimization"""
# Keep track of evals
global ITERATION
ITERATION += 1
# # Retrieve the subsample if present otherwise set to 1.0
# subsample = params['boosting_type'].get('subsample', 1.0)
# # Extract the boosting type
# params['boosting_type'] = params['boosting_type']['boosting_type']
# params['subsample'] = subsample
# Make sure parameters that need to be integers are integers
for parameter_name in ['max_depth', 'min_child_weight']:
params[parameter_name] = int(params[parameter_name])
start = timer()
# Perform n_folds cross validation
cv_results = xgb.cv(params, train_set, num_boost_round = 500, folds = cv_split,
early_stopping_rounds = 50, seed = 50, feval=custom_mae_xgb)
run_time = timer() - start
# Extract the best score, Loss must be minimized
best_score = np.min(cv_results['test-custom_mae_xgb-mean'])
loss = best_score
# Boosting rounds that returned the highest cv score
n_estimators = int(np.argmin(cv_results['test-custom_mae_xgb-mean']) + 1)
# Dictionary with information for evaluation
return {'loss': best_score, 'params': params, 'iteration': ITERATION,
'estimators': n_estimators,
'train_time': run_time, 'status': STATUS_OK}
# Define the search space
space = {
'booster': hp.choice('booster', ['gbtree', 'gblinear', 'dart']),
'max_depth': hp.quniform('max_depth', 1, 10, 1),
'learning_rate': hp.loguniform('learning_rate', np.log(0.01), np.log(0.2)),
'gamma': hp.quniform('gamma', 0.5,5,0.5),
'min_child_weight': hp.quniform('min_child_weight', 1, 20, 1),
'subsample': hp.uniform('subsample', 0.5, 1),
'colsample_bytree': hp.uniform('colsample_by_tree', 0.6, 1.0),
'reg_alpha': hp.uniform('reg_alpha', 0.0, 1.0),
'reg_lambda': hp.uniform('reg_lambda', 0.0, 1.0)
}
# optimization algorithm
tpe_algorithm = tpe.suggest
# Keep track of results
bayes_trials = Trials()
# Global variable
global ITERATION
ITERATION = 0
from hyperopt import fmin
lgb_best = fmin(fn = objective, space = space, algo = tpe.suggest,
max_evals = MAX_EVALS, trials = bayes_trials, rstate = np.random.RandomState(10))
xgb_trials_results = sorted(bayes_trials.results, key = lambda x: x['loss'])
xgb_results = pd.DataFrame(xgb_trials_results)
best_xgb_estimators= xgb_results.loc[0,'estimators']
best_xgb_params= xgb_results.loc[0,'params']
best_xgb_model = xgb.XGBRegressor(n_estimators=best_xgb_estimators, n_jobs = -1,
random_state = 50, **best_xgb_params)
RandomForest с Hyperopt
MAX_EVALS = 200
def objective(params):
"""Objective function for Gradient Boosting Machine Hyperparameter Optimization"""
# Keep track of evals
global ITERATION
ITERATION += 1
params['oob_score'] = params['bootstrap']['oob_score']
params['bootstrap'] = params['bootstrap']['bootstrap']
for parameter_name in ['min_samples_leaf', 'min_samples_split', 'n_estimators']:
params[parameter_name] = int(params[parameter_name])
start = timer()
rf = RandomForestRegressor(**params)
cv_results = cross_validate(rf, X_train, y_train, cv = cv_split, scoring=custom_mae)
run_time = timer() - start
best_score = -cv_results['test_score'].mean()
loss = best_score
return {'loss': best_score, 'params': params, 'iteration': ITERATION,
'train_time': run_time, 'status': STATUS_OK}
space = {
'bootstrap': hp.choice('bootstrap', [{'bootstrap': True, 'oob_score': True},
{'bootstrap': True, 'oob_score': False}]),
'max_features': hp.choice('max_features', ['auto', 'sqrt']),
'max_depth': hp.choice('max_depth', [None,1,2,3]),
'min_samples_split': hp.quniform('min_samples_split', 2,20,1),
'min_samples_leaf': hp.quniform('min_samples_leaf ', 1,10,1),
'n_estimators': hp.quniform('n_estimators ',20,1000,10)
}
tpe_algorithm = tpe.suggest
bayes_trials = Trials()
ITERATION = 0
rf_best = fmin(fn = objective, space = space, algo = tpe.suggest,
max_evals = MAX_EVALS, trials = bayes_trials, rstate = np.random.RandomState(10))
rf_trials_results = sorted(bayes_trials.results, key = lambda x: x['loss'])
rf_results = pd.DataFrame(rf_trials_results)
best_rf_params= rf_results.loc[0,'params']
best_rf_model = RandomForestRegressor(random_state = 50, **best_rf_params)
Модели усреднения и суммирования
Усредненный базовый класс моделей
Начнем с этого простого подхода для усреднения базовых моделей. Мы создаем новый класс, который расширяет scikit-learn с помощью нашей модели, что также позволяет нам повторно использовать код и использовать инкапсуляцию.
from sklearn.base import RegressorMixin, TransformerMixin, BaseEstimator, clone
class AveragingModels(BaseEstimator, RegressorMixin, TransformerMixin):
def __init__(self, models):
self.models = models
# we define clones of the original models to fit the data in
def fit(self, X, y):
self.models_ = [clone(x) for x in self.models]
# Train cloned base models
for model in self.models_:
model.fit(X, y)
return self
#Now we do the predictions for cloned models and average them
def predict(self, X):
predictions = np.column_stack([
model.predict(X) for model in self.models_
])
return np.mean(predictions, axis=1)
averaged_models = AveragingModels(models=(best_lgb_model, best_xgb_model, best_rf_model))
averaged_results = cross_validate(averaged_models, X_train, y_train, cv = cv_split, scoring=custom_mae)
print('Cross-Validated score on train data: ', -averaged_results['test_score'].mean())
Перекрестная оценка данных о поездах: 97875,9151127886.
StackingCVRegressor
# Read for details: http://rasbt.github.io/mlxtend/user_guide/regressor/StackingCVRegressor/#targetText=Stacking%20is%20an%20ensemble%20learning,for%20the%20level%2D2%20regressor.
class StackingAveragedModels(BaseEstimator, RegressorMixin, TransformerMixin):
def __init__(self, base_models, meta_model, cv=cv_split):
self.base_models = base_models
self.meta_model = meta_model
self.cv = cv
# We again fit the data on clones of the original models
def fit(self, X, y):
self.base_models_ = [list() for x in self.base_models]
self.meta_model_ = clone(self.meta_model)
# Train cloned base models then create out-of-fold predictions
# that are needed to train the cloned meta-model
out_of_fold_predictions = np.zeros((X.shape[0], len(self.base_models)))
for i, model in enumerate(self.base_models):
for train_index, holdout_index in cv_split.split(X, y):
instance = clone(model)
self.base_models_[i].append(instance)
instance.fit(X.iloc[train_index], y.iloc[train_index])
y_pred = instance.predict(X.iloc[holdout_index])
out_of_fold_predictions[holdout_index, i] = y_pred
# Now train the cloned meta-model using the out-of-fold predictions as new feature
self.meta_model_.fit(out_of_fold_predictions, y)
return self
#Do the predictions of all base models on the test data and use the averaged predictions as
#meta-features for the final prediction which is done by the meta-model
def predict(self, X):
meta_features = np.column_stack([
np.column_stack([model.predict(X) for model in base_models]).mean(axis=1)
for base_models in self.base_models_ ])
return self.meta_model_.predict(meta_features)
stacked_models = StackingAveragedModels(base_models=(best_lgb_model, best_xgb_model, best_rf_model),
meta_model = LassoCV())
stacked_models.fit(X_train, y_train)
stacked_models_prediction = stacked_models.predict(X_train)
print('Cross-Validated score on train data: ', custom_mae_(y_train, stacked_models_prediction))
Перекрестная проверка данных поезда: 75813,51533323177
Окончательный прогноз
stacked_models_prediction = stacked_models.predict(X_test) custom_mae_(y_test, stacked_models_prediction) 97889.44897178133
В заключение я считаю, что проект удался! Результаты тестовых данных намного превзошли результаты обученных с параметрами по умолчанию. Объединение моделей с отдельной метамоделью помогло нам получить относительно высокий балл.
Не стесняйтесь связаться со мной через:
Электронная почта: [email protected]
Linkedin: https://www.linkedin.com/in/slava-spirin/
Гитхаб: https://github.com/slavaspirin
