Описание проекта¶
Допустим, вы работаете в добывающей компании «ГлавРосГосНефть». Нужно решить, где бурить новую скважину.
Вам предоставлены пробы нефти в трёх регионах: в каждом 10 000 месторождений, где измерили качество нефти и объём её запасов. Постройте модель машинного обучения, которая поможет определить регион, где добыча принесёт наибольшую прибыль. Проанализируйте возможную прибыль и риски техникой Bootstrap.
Шаги для выбора локации:
- В избранном регионе ищут месторождения, для каждого определяют значения признаков;
- Строят модель и оценивают объём запасов;
- Выбирают месторождения с самым высокими оценками значений. Количество месторождений зависит от бюджета компании и стоимости разработки одной скважины;
- Прибыль равна суммарной прибыли отобранных месторождений.
Условия задачи:¶
- Для обучения модели подходит только линейная регрессия (остальные — недостаточно предсказуемые).
- При разведке региона проводится исследование 500 точек.
- Бюджет на разработку месторождений — 10 млрд рублей, стоимость бурения одной скважины — 50 млн рублей.
- Один баррель сырья приносит 4500 рублей прибыли.
- Не рассматривать регионы, в которых риск убытков выше 2.5%. Из оставшихся выбирается регион с наибольшей средней прибылью.
Оглавление¶
- I. Загрузка и подготовка данных
- II. Обучение и проверка модели
- III. Подготовка к расчёту прибыли
- IV. Расчёт прибыли и рисков
1) Объявляем класс Explorer и наполняем его методами¶
In [1]:
class Explorer:
def histogram(self, data, n_bins, range_start, range_end, grid, cumulative=False, x_label = "", y_label = "", title = ""):
"""
data - датасет
n_bins - количество корзин
range_start - минимальный икс для корзины
range_end - максимальный икс для корзины
grid - рисовать сетку или нет (False / True)
histogram(data, n_bins, range_start, range_end, grid, x_label = "", y_label = "", title = "")
Пример:
histogram(df, 100, 0, 150, True, 'Количество иксов', 'Количество игриков', 'Заголовок')
"""
# Создаем объект - график
_, ax = plt.subplots()
# Задаем параметры
ax.hist(data, bins = n_bins, range = (range_start, range_end), cumulative = cumulative, color = '#4169E1')
# Добавляем сетку
if grid == True:
ax.grid(color='grey', linestyle='-', linewidth=0.5)
else:
pass
# Добавляем медиану, среднее и квартили
ax.axvline(data.median(),linestyle = '--', color = '#FF1493', label = 'median')
ax.axvline(data.mean(),linestyle = '--', color = 'orange', label = 'mean')
ax.axvline(data.quantile(0.1),linestyle = '--', color = 'yellow', label = '1%')
ax.axvline(data.quantile(0.99),linestyle = '--', color = 'yellow', label = '99%')
ax.legend()
ax.set_ylabel(y_label)
ax.set_xlabel(x_label)
ax.set_title(title)
def scatterplot(self, x_data, y_data, x_label="", y_label="", title="", color = "r", yscale_log=False, figsize = (8, 6)):
# Создаем объект - график
_, ax = plt.subplots(figsize = (8, 6))
# Задаем параметры для графика, определяем размер (s), цвет и прозрачность точек на графике
ax.scatter(x_data, y_data, s = 10, color = color, alpha = 0.75)
if yscale_log == True:
ax.set_yscale('log')
# Создаем описание осей и заголовок для графика
ax.set_title(title)
ax.set_xlabel(x_label)
ax.set_ylabel(y_label)
def overlaid_histogram(self, data1, data2, n_bins = 0, data1_name="", data1_color="#539caf", data2_name="", data2_color="#7663b0", x_label="", y_label="", title=""):
# Устанавливаем границы для корзин так чтобы оба распределения на графике были соотносимы
max_nbins = 10
data_range = [min(min(data1), min(data2)), max(max(data1), max(data2))]
binwidth = (data_range[1] - data_range[0]) / max_nbins
if n_bins == 0:
bins = np.arange(data_range[0], data_range[1] + binwidth, binwidth)
else:
bins = n_bins
# Рисуем график
_, ax = plt.subplots(figsize=(10,8))
ax.hist(data1, bins = bins, color = data1_color, alpha = 0.65, label = data1_name)
ax.hist(data2, bins = bins, color = data2_color, alpha = 0.65, label = data2_name)
ax.axvline(data1.mean(),linestyle = '--', color = 'lime', label = 'mean for data 1')
ax.axvline(data2.mean(),linestyle = '--', color = 'coral', label = 'mean for data 2')
ax.set_ylabel(y_label)
ax.set_xlabel(x_label)
ax.set_title(title)
ax.legend(loc = 'best')
def corr_diagram(self, x):
plt.figure(figsize=(12,10), dpi= 80)
sns.heatmap(x.corr(), xticklabels=x.corr().columns, yticklabels=x.corr().columns, cmap='RdYlGn', center=0, annot=True)
plt.title('Диаграмма корреляции', fontsize=22)
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
plt.show()
def highlight_max(self, data, color='#00FF00'):
'''
highlight the maximum in a Series or DataFrame
'''
attr = 'background-color: {}'.format(color)
#remove % and cast to float
data = data.replace('%','', regex=True).astype(float)
data[data == 1] = None
if data.ndim == 1: # Series from .apply(axis=0) or axis=1
is_max = (data == data.abs().max()) & (data !=1)
return [attr if v else '' for v in is_max]
else: # from .apply(axis=None)
is_max = (data == data.abs().max()) & (data !=1)
return pd.DataFrame(np.where(is_max, attr, ''),
index=data.index, columns=data.columns)
def highlight_sorted_corr(self, data, color='#00FF00'):
'''
highlight the maximum in a Series or DataFrame
'''
attr = 'background-color: {}'.format(color)
#remove % and cast to float
data = data.replace('%','', regex=True).astype(float)
data[data == 1] = None
if data.ndim == 1: # Series from .apply(axis=0) or axis=1
is_max = (data > 0.1) & (data !=1)
return [attr if v else '' for v in is_max]
else: # from .apply(axis=None)
is_max = (data == data.abs().max()) & (data !=1)
return pd.DataFrame(np.where(is_max, attr, ''),
index=data.index, columns=data.columns)
def lineplot(self, x_data, y_data, x_label="", y_label="", title=""):
# Создаем объект - график
_, ax = plt.subplots(figsize=(8, 6))
# Задаем параметры для линии: ширину (lw), цвет и прозрачность (alpha)
ax.plot(x_data, y_data, lw = 2, color = '#539caf', alpha = 1)
# Даем имена осям и заголовок для графика
ax.set_title(title)
ax.set_xlabel(x_label)
ax.set_ylabel(y_label)
def double_lineplot(self, x_data_1, y_data_1, x_data_2, y_data_2, x_label="", y_label="", title="", label_one="", label_two=""):
# Создаем объект - график
_, ax = plt.subplots(figsize=(8, 6))
# Задаем параметры для линии: ширину (lw), цвет и прозрачность (alpha)
ax.plot(x_data_1, y_data_1, lw = 2, color = '#6400e4', alpha = 1, label = label_one)
ax.plot(x_data_2, y_data_2, lw = 2, color = '#ffc740', alpha = 1, label = label_two)
# Даем имена осям и заголовок для графика
ax.set_title(title)
ax.set_xlabel(x_label)
ax.set_ylabel(y_label)
ax.legend(loc = 'best')
def hexbin(self, data, x, y):
data.plot(x = x, y = y, kind='hexbin', gridsize=20, figsize=(15, 10), sharex=False, grid=True)
def bar_plotter(self, data):
data.plot.bar(rot=0, figsize = (16, 5))
def categorical_counter_plot(self, data, column, x = '', y = ''):
if x == '' or y == '':
plt.rcParams["figure.figsize"] = (15, 10)
else:
plt.rcParams["figure.figsize"] = (x, y)
order = data[column].value_counts().index
ax = sns.countplot(data[column], order = order)
ax.set_xticklabels(ax.get_xticklabels(), fontsize=11)
plt.xticks(rotation=90)
def sns_scatterplot(self, data, x="", y="", hue="", size="", palette=""):
sns.set(style="whitegrid")
f, ax = plt.subplots(figsize=(15, 10))
if palette == True:
sns.scatterplot(ax = ax, x=x, y=y, palette="ch:r=-.2,d=.3_r",
hue=hue, size=size, sizes=(1, 200), linewidth=0, data=data)
else:
sns.scatterplot(ax = ax, x=x, y=y,
hue=hue, size=size,
sizes=(1, 200), linewidth=0, data=data)
def sns_catplot(self, data, x="", y="", hue=""):
sns.set(style='whitegrid')
sns.catplot(x=x, y=y, hue=hue, kind='bar', errwidth=0,
data=data, height=5, aspect=3)
def sns_factorplot(self, data, x='', hue=''):
sns.axes_style('white')
g = sns.factorplot("exited", data=df, aspect=1, kind='count',
hue='hascrcard')
def squared_ratio(self, df, grouper, title=''):
df = df.groupby(grouper).size().reset_index(name='counts')
labels = df.apply(lambda x: str(x[0]) + "\n (" + str(x[1]) + ")", axis=1)
sizes = df['counts'].values.tolist()
colors = [plt.cm.Spectral(i/float(len(labels))) for i in range(len(labels))]
plt.figure(figsize=(10,6), dpi= 80)
squarify.plot(sizes=sizes, label=labels, color=colors, alpha=.8)
plt.title(title)
plt.axis('off')
plt.show()
def sorted_corr(self, data, attr):
correlated = pd.DataFrame(data.corr()[attr].sort_values(ascending = False))
return correlated
def transformer(self, data, name, grouper, func):
"""
data - датасет
name - столбец в котором меняем значения
grouper - столбец по которому группируем
func - пременяемая функция mean, median и т.д.
"""
name = name
data.loc[data[name].isnull(), name] = data.groupby(grouper)[name].transform(func)
def pr_curve(self, model, features_valid, target_valid):
probabilities_valid = model.predict_proba(features_valid)
precision, recall, thresholds = precision_recall_curve(target_valid, probabilities_valid[:, 1])
plt.figure(figsize=(6, 6))
plt.step(recall, precision, where='post')
plt.xlabel('Recall')
plt.ylabel('Precision')
plt.ylim([0.0, 1.05])
plt.xlim([0.0, 1.0])
plt.title('Кривая Precision-Recall')
plt.show()
def roc_curve(self, model, features_valid, target_valid):
probabilities_valid = model.predict_proba(features_valid)
probabilities_one_valid = probabilities_valid[:, 1]
fpr, tpr, thresholds = roc_curve(target_valid, probabilities_one_valid)
plt.figure()
plt.plot(fpr, tpr)
plt.plot([0, 1], [0, 1], linestyle='--')
plt.xlim(0,1)
plt.ylim(0,1)
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('ROC-кривая')
plt.show()
def metrics_plot(self, model, features_valid, target_valid):
probabilities_valid = model.predict_proba(features_valid)
precision, recall, thresholds = precision_recall_curve(target_valid, probabilities_valid[:, 1])
fpr, tpr, thresholds = roc_curve(target_valid, probabilities_valid[:, 1])
fig, ax = plt.subplots(1, 2, figsize=(15, 6))
fig.subplots_adjust(hspace=0.4, wspace=0.4)
#fig, ax = plt.subplots(ncols=3)
#fig.subplots_adjust(hspace=0.4, wspace=0.4)
sns.lineplot(recall, precision, drawstyle='steps-post', ax=ax[0])
ax[0].set_xlabel('Recall')
ax[0].set_ylabel('Precision')
ax[0].set_ylim([0.0, 1.05])
ax[0].set_xlim([0.0, 1.0])
ax[0].set_title('Кривая Precision-Recall')
sns.lineplot(fpr, tpr, ax=ax[1])
ax[1].plot([0, 1], [0, 1], linestyle='--')
ax[1].set_xlim(0,1)
ax[1].set_ylim(0,1)
ax[1].set_xlabel('False Positive Rate')
ax[1].set_ylabel('True Positive Rate')
ax[1].set_title('ROC-кривая')
def auc_roc(self, model, features_valid, target_valid):
probabilities_valid = model.predict_proba(features_valid)
probabilities_one_valid = probabilities_valid[:, 1]
auc_roc = roc_auc_score(target_valid, probabilities_one_valid)
return auc_roc
def upsample(self, features, target, repeat):
features_zeros = features[target == 0]
features_ones = features[target == 1]
target_zeros = target[target == 0]
target_ones = target[target == 1]
features_upsampled = pd.concat([features_zeros] + [features_ones] * repeat)
target_upsampled = pd.concat([target_zeros] + [target_ones] * repeat)
features_upsampled, target_upsampled = shuffle(features_upsampled, target_upsampled, random_state=42)
return features_upsampled, target_upsampled
def downsample(self, features, target, fraction):
features_zeros = features[target == 0]
features_ones = features[target == 1]
target_zeros = target[target == 0]
target_ones = target[target == 1]
features_downsampled = pd.concat([features_zeros.sample(frac=fraction, random_state=42)] + [features_ones])
target_downsampled = pd.concat([target_zeros.sample(frac=fraction, random_state=42)] + [target_ones])
features_downsampled, target_downsampled = shuffle(features_downsampled, target_downsampled, random_state=42)
return features_downsampled, target_downsampled
def firstsight(self, data):
head = data.head(5)
tail = data.tail(5)
sample = data.tail(5)
return head, tail, sample
def smape(self, y_test, y_predict):
y_test, y_predict = np.array(y_test), np.array(y_predict)
return np.mean((np.abs((y_predict - y_test)) / np.mean(np.abs(y_test) + np.abs(y_predict + 0.1**99))) * 100)
def mape(self, y_test, y_predict):
y_test, y_predict = np.array(y_test), np.array(y_predict)
return np.median((np.abs((y_test - y_predict)) / (y_test + 0.1**99)) * 100)
def df_split(self, data, features_drop, target, test_size, random_state):
feature = data.drop(features_drop, axis=1)
target = data[target]
x_train, x_test, y_train, y_test = train_test_split(feature, target, test_size=test_size, random_state = random_state)
return x_train, x_test, y_train, y_test
def grid_search(self, model, param_grid, cv, x, y):
grid_model = GridSearchCV(model, param_grid=param_grid, cv=cv, verbose=1, n_jobs=-1)
grid_model.fit(x, y)
best_estimator = grid_model.best_estimator_
return best_estimator
class Display(object):
"""Выводит HTML представление нескольких объектов"""
template = """<div style="float: left; padding: 10px;">
<p style='font-family:"Courier New", Courier, monospace'>{0}</p>{1}
</div>"""
def __init__(self, *args):
self.args = args
def _repr_html_(self):
return '\n'.join(self.template.format(a, eval(a)._repr_html_())
for a in self.args)
def __repr__(self):
return '\n\n'.join(a + '\n' + repr(eval(a))
for a in self.args)
In [2]:
explorer = Explorer()
(... загрузка библиотек и данных ...)
4) Выведем первые пять строк каждой из таблиц, общую информацию¶
- регион 1¶
In [5]:
head, tail, sample = explorer.firstsight(df_geo_1)
explorer.Display('head', 'tail', 'sample')
Out[5]:
In [6]:
df_geo_1.info()
⚡ Вывод: в датасете с информацией по первому региону 100 000 строк, без пропусков, без проблем с типами данных.
- регион 2¶
In [7]:
head, tail, sample = explorer.firstsight(df_geo_2)
explorer.Display('head', 'tail', 'sample')
Out[7]:
In [8]:
df_geo_2.info()
⚡ Вывод: в датасете с информацией по второму региону 100 000 строк, без пропусков, без проблем с типами данных.
- регион 3¶
In [9]:
head, tail, sample = explorer.firstsight(df_geo_3)
explorer.Display('head', 'tail', 'sample')
Out[9]:
In [10]:
df_geo_3.info()
⚡ Вывод: в датасете с информацией по третьему региону 100 000 строк, без пропусков, без проблем с типами данных.
5) Выведем общую статистку по датасетам, проверим наличие дубликатов и уникальные значения¶
- регион 1¶
In [11]:
df_geo_1.describe().T
Out[11]:
- регион 2¶
In [12]:
df_geo_2.describe().T
Out[12]:
- регион 3¶
In [13]:
df_geo_3.describe().T
Out[13]:
С учетом характера данных - общая статистика малоинформативна без точного понимания, какие именно атрибуты наблюдений содержаться в столбцах f0, f1, f2.
- проверим наличие полных дубликатов¶
In [14]:
(df_geo_1.duplicated() == True).sum()
Out[14]:
In [15]:
(df_geo_2.duplicated() == True).sum()
Out[15]:
In [16]:
(df_geo_3.duplicated() == True).sum()
Out[16]:
Полных дубликатов нет
- проверим уникальные значения¶
In [17]:
df_geo_1.nunique()
Out[17]:
In [18]:
df_geo_2.nunique()
Out[18]:
In [19]:
df_geo_3.nunique()
Out[19]:
⚡ Вывод:
- Описательная статистика неинформативна без точного описания данных без точного понимания, какие именно атрибуты наблюдений содержаться в столбцах f0, f1, f2.
- В датасетах нет полных дубликатов.
- В каждом из датасетов в столбце id есть повторяющиеся значения: df_geo_1 - 10 из 100 000 (0,0001%), df_geo_2 - 4 из 100 000 (0,00004%), df_geo_3 - 4 из 100 000 (0,00004%). Вероятно, из некторых скважин брали несколько замеров, скорее всего на разных глубинах.
- Исходя из описания данных, столбец id кандидат на удаление, поскольку уникальный идентификатор месторождения - не универсальный признак.
Вывод¶
Исходя из представленных данных и описаний данных, можно сделать следующие выводы:
- В датасетах представлены данные георазведки из трех регионов по 100 000 уникальных наблюдений в каждом.
- Во всех датасетах нет проблем с типами данных, полные дубликаты отсутствуют, нет пропусков.
- Есть повторяющиеся значения в столбцах с id: df_geo_1 - 10 из 100 000 (0,0001%), df_geo_2 - 4 из 100 000 (0,00004%), df_geo_3 - 4 из 100 000 (0,00004%).
- Поскольку нет возможности уточнить у источника данных - исходим из буквального описания датасета: есть наблюдения из одних тех же скважин.
- Столбец id не будет использован для обучения модели, поскольку является неуниверсальным атрибутом.
- создадим список в котом будем фиксировать обозначение модели, результаты метрик, средний запас сырья по предсказанию модели¶
In [20]:
results = []
In [21]:
lr = LinearRegression(n_jobs=-1)
- параметры для модели будут подбираться с помощью GridSearchCV, поэтому определим словарь с параметрами¶
In [22]:
param_grid={'fit_intercept':[True, False], 'normalize':[True, False], 'copy_X':[True, False]}
1) Разделим данные на обучающую и валидационную выборки в соотношении 75:25¶
In [23]:
x_train, x_test, y_train, y_test = explorer.df_split(df_geo_1, ['id', 'product'], 'product', 0.25, 42)
- проверим форму полученных выборок¶
In [24]:
x_train.shape, y_train.shape
Out[24]:
In [25]:
x_test.shape, y_test.shape
Out[25]:
2) Обучим модель и сделаем предсказания на валидационной выборке¶
In [26]:
lr_geo_one = explorer.grid_search(lr, param_grid, 5, x_train, y_train)
- сохраним предсказания в переменной¶
In [27]:
predicted_product = list(lr_geo_one.predict(x_test))
3) Cохраним предсказания модели и правильные ответы на валидационной выборке¶
In [28]:
region_1_report = pd.DataFrame(y_test)
- сохраним предсказания модели по каждому месторождению и пометим регион, для которого сделали предсказание¶
In [29]:
region_1_report['predicted'] = predicted_product
region_1_report['region'] = 1
- выведем для проверки первые пять строк¶
In [30]:
region_1_report.head(5)
Out[30]:
- выведем график с сотами, чтобы посмотреть пересечения между предсказаниями и правильными ответами¶
In [31]:
explorer.hexbin(region_1_report, 'predicted', 'product')
4) Выведем на экране средний запас сырья, RMSE модели и другие метрики¶
- средний запас сырья в регионе, предсказанный моделью¶
In [32]:
predicted_mean_product = region_1_report.predicted.mean()
predicted_mean_product
Out[32]:
- mean squared error (MSE)¶
In [33]:
mse = mean_squared_error(y_test, predicted_product)
mse
Out[33]:
- root-mean-square error (RMSE)¶
In [34]:
rmse = np.sqrt(mse)
rmse
Out[34]:
- mean absolute percentage error (MAPE)¶
In [35]:
mape = explorer.mape(y_test, predicted_product)
mape
Out[35]:
5) Дополним отчет¶
In [36]:
results.append(('lr_geo_one', lr_geo_one.score(x_test, y_test), rmse, mape, predicted_mean_product, lr_geo_one))
⚡ Вывод:¶
- Cредний запас сырья в регионе 1, предсказанный моделью - 92,4 тыс. баррелей.
- Среднеквадратичная ошибка (MSE) - 1426.
- Квадратный корень из средней квадратичной ошибки (RMSE) - 37.76.
- Средний процент отклонения (MAPE) - 28 %. Модель может ошибаться в среднем на 28%, что довольно много.
1) Разделим данные на обучающую и валидационную выборки в соотношении 75:25¶
In [37]:
x_train, x_test, y_train, y_test = explorer.df_split(df_geo_2, ['id', 'product'], 'product', 0.25, 42)
- проверим форму полученных выборок¶
In [38]:
x_train.shape, y_train.shape
Out[38]:
In [39]:
x_test.shape, y_test.shape
Out[39]:
2) Обучим модель и сделаем предсказания на валидационной выборке¶
In [40]:
lr_geo_two = explorer.grid_search(lr, param_grid, 5, x_train, y_train)
- сохраним предсказания в переменной¶
In [41]:
predicted_product = list(lr_geo_two.predict(x_test))
3) Cохраним предсказания модели и правильные ответы на валидационной выборке¶
In [42]:
region_2_report = pd.DataFrame(y_test)
- сохраним предсказания модели по каждому месторождению и пометим регион, для которого сделали предсказание¶
In [43]:
region_2_report['predicted'] = predicted_product
region_2_report['region'] = 2
In [44]:
region_2_report.head(5)
Out[44]:
- выведем график с сотами, чтобы посмотреть пересечения между предсказаниями и правильными ответами¶
In [45]:
explorer.hexbin(region_2_report, 'predicted', 'product')
4) Выведем на экране средний запас сырья и RMSE модели¶
In [46]:
predicted_mean_product = region_2_report.predicted.mean()
predicted_mean_product
Out[46]:
In [47]:
mse = mean_squared_error(y_test, predicted_product)
mse
Out[47]:
In [48]:
rmse = np.sqrt(mse)
rmse
Out[48]:
In [49]:
mape = explorer.mape(y_test, predicted_product)
mape
Out[49]:
5) Дополним отчет¶
In [50]:
results.append(('lr_geo_two', lr_geo_two.score(x_test, y_test), rmse, mape, predicted_mean_product, lr_geo_two))
⚡ Вывод:¶
- Cредний запас сырья в регионе 2, предсказанный моделью - 68.71 тыс. баррелей.
- Среднеквадратичная ошибка (MSE) - 0.79.
- Квадратный корень из средней квадратичной ошибки (RMSE) - 0.89.
- Средний процент отклонения (MAPE) - 1 %.
- Высокий результат модели связан со странностями распределния выборки по второму региону 2. Требуется уточнение у поставщика данных, почему целевые показатели сгруппированы относительно друг друга практически в равных интервалах.
- С учетом исследования, результаты модели следует рассматривать с крайней осторожностью. Может быть, даже вообще не рассматривать регион 2.
1) Разделим данные на обучающую и валидационную выборки в соотношении 75:25¶
In [51]:
x_train, x_test, y_train, y_test = explorer.df_split(df_geo_3, ['id', 'product'], 'product', 0.25, 42)
- проверим форму полученных выборок¶
In [52]:
x_train.shape, y_train.shape
Out[52]:
In [53]:
x_test.shape, y_test.shape
Out[53]:
2) Обучим модель и сделаем предсказания на валидационной выборке¶
In [54]:
lr_geo_three = explorer.grid_search(lr, param_grid, 5, x_train, y_train)
- сохраним предсказания в переменной¶
In [55]:
predicted_product = list(lr_geo_three.predict(x_test))
3) Cохраним предсказания модели и правильные ответы на валидационной выборке¶
In [56]:
region_3_report = pd.DataFrame(y_test)
- сохраним предсказания модели по каждому месторождению и пометим регион, для которого сделали предсказание¶
In [57]:
region_3_report['predicted'] = predicted_product
region_3_report['region'] = 3
In [58]:
region_3_report.head(5)
Out[58]:
- выведем график с сотами, чтобы посмотреть пересечения между предсказаниями и правильными ответами¶
In [59]:
explorer.hexbin(region_3_report, 'predicted', 'product')
4) Выведем на экране средний запас сырья и RMSE модели¶
In [60]:
predicted_mean_product = region_3_report.predicted.mean()
predicted_mean_product
Out[60]:
In [61]:
mse = mean_squared_error(y_test, predicted_product)
mse
Out[61]:
In [62]:
rmse = np.sqrt(mse)
rmse
Out[62]:
In [63]:
mape = explorer.mape(y_test, predicted_product)
mape
Out[63]:
5) Дополним отчет¶
In [64]:
results.append(('lr_geo_three', lr_geo_three.score(x_test, y_test), rmse, mape, predicted_mean_product, lr_geo_three))
⚡ Вывод:¶
- Cредний запас сырья в регионе 1, предсказанный моделью - 94.77 тыс. баррелей.
- Среднеквадратичная ошибка (MSE) - 1611.
- Квадратный корень из средней квадратичной ошибки (RMSE) - 40.15.
- Средний процент отклонения (MAPE) - 28.9 %.
- Результаты модели невысокие, чуть хуже модели, обученной на данных для первого региона.
5. Результаты моделей для каждого из регионов¶
In [65]:
results_df = pd.DataFrame(results, columns=['model', 'r2', 'rmse', 'mape', 'predicted_mean_product', 'params'])
In [66]:
results_df[['model', 'r2', 'rmse', 'mape', 'predicted_mean_product']].style.highlight_max(color='lightgreen', axis = 0)
Out[66]:
⚡ Вывод:
6. Сравнение правильного среднего и среднего, предсказанного моделью¶
In [67]:
region_report = pd.concat([region_1_report, region_2_report, region_3_report])
region_report.pivot_table(index='region')
Out[67]:
- построим гистограммы с одинаковыми границами бинов для каждого из регионов¶
In [68]:
explorer.overlaid_histogram(region_1_report['product'], region_1_report['predicted'], 150,
data1_name='product-data1', data1_color='red',
data2_name='predicted-data2', data2_color='royalblue',
title='Регион 1')
In [69]:
explorer.overlaid_histogram(region_2_report['product'], region_2_report['predicted'], 150,
data1_name='product-data1', data1_color='red',
data2_name='predicted-data2', data2_color='royalblue',
title='Регион 2')
In [70]:
explorer.overlaid_histogram(region_3_report['product'], region_3_report['predicted'], 150,
data1_name='product-data1', data1_color='red',
data2_name='predicted-data2', data2_color='royalblue',
title='Регион 3')
⚡ Вывод: числовые средние двух выборок близки друг к другу до степени совпадения. Равенство средних еще не говорит о принадлежности выборок одной и той же генеральной совокупности. С учетом значений rmse и r2 моделей, принимая во внимание близость средних предсказанных и правильных ответов, опираясь на визуализированные данные, полагаю необходимым проверить однородность распределений для каждой из выборок.
1) Проверим гипотезу об однородности распределений по критерию Смирнова в отношении региона 1¶
- сформулируем нулевую гипотезу¶
1) Нулевую гипотезу, которую проверяет тест, сформулируем так: "Выборка предсказанных моделью объемов сырья и выборка правильных объемов сырья принадлежат одному закону распределения".
2) Альтернативную гипотезу, соответственно, сформулируем так: "Выборка предсказанных моделью объемов сырья и выборка правильных объемов сырья принадлежат разным законам распределения".
- примем уровень значимости в 5%.¶
In [71]:
alpha = 0.05
- проведем тест и получим p-значение¶
In [72]:
ks_stat, pval_ks = st.ks_2samp(region_1_report['product'], region_1_report['predicted'])
pval_ks
Out[72]:
- cравним p-значение с уровнем значимости¶
In [73]:
pval_ks < alpha
Out[73]:
⚡ Вывод: исходя из представленных данных, на уровне значимости 5% - есть основания отвергнуть нулевую гипотезу в пользу альтернативы. Выборка предсказанных моделью объемов сырья и выборка правильных объемов сырья в отношении региона 1 принадлежат разным законам распределения.
2) Проверим гипотезу об однородности распределений по критерию Смирнова в отношении региона 2¶
- сформулируем нулевую гипотезу¶
1) Нулевую гипотезу, которую проверяет тест, сформулируем так: "Выборка предсказанных моделью объемов сырья и выборка правильных объемов сырья принадлежат одному закону распределения".
2) Альтернативную гипотезу, соответственно, сформулируем так: "Выборка предсказанных моделью объемов сырья и выборка правильных объемов сырья принадлежат разным законам распределения".
- примем уровень значимости в 5%.¶
In [74]:
alpha = 0.05
- проведем тест и получим p-значение¶
In [75]:
ks_stat, pval_ks = st.ks_2samp(region_2_report['product'], region_2_report['predicted'])
pval_ks
Out[75]:
- cравним p-значение с уровнем значимости¶
In [76]:
pval_ks < alpha
Out[76]:
⚡ Вывод: исходя из представленных данных, на уровне значимости 5% - есть основания отвергнуть нулевую гипотезу в пользу альтернативы. Выборка предсказанных моделью объемов сырья и выборка правильных объемов сырья в отношении региона 2 принадлежат разным законам распределения.
3) Проверим гипотезу об однородности распределений по критерию Смирнова в отношении региона 3¶
- сформулируем нулевую гипотезу¶
1) Нулевую гипотезу, которую проверяет тест, сформулируем так: "Выборка предсказанных моделью объемов сырья и выборка правильных объемов сырья принадлежат одному закону распределения".
2) Альтернативную гипотезу, соответственно, сформулируем так: "Выборка предсказанных моделью объемов сырья и выборка правильных объемов сырья принадлежат разным законам распределения".
- примем уровень значимости в 5%.¶
In [77]:
alpha = 0.05
- проведем тест и получим p-значение¶
In [78]:
ks_stat, pval_ks = st.ks_2samp(region_3_report['product'], region_3_report['predicted'])
pval_ks
Out[78]:
- cравним p-значение с уровнем значимости¶
In [79]:
pval_ks < alpha
Out[79]:
⚡ Вывод: исходя из представленных данных, на уровне значимости 5% - есть основания отвергнуть нулевую гипотезу в пользу альтернативы. Выборка предсказанных моделью объемов сырья и выборка правильных объемов сырья в отношени региона 3 принадлежат разным законам распределения.
- В отношении каждого из регионов построены модели линейной регрессии со следующими результатами RMSE:
- регион 1 - 37.76,
- регион 2 - 0.89,
- регион 3 - 40.14.
- Результаты работы моделей в отношении региона 1 и региона 2 - неудовлетворительные. Результаты работы модели в отношении региона 2 связан с особенностями распределения данных: через 6 областей через равные интервалы.
- Поэтому к предсказаниям всех трех моделей следует относится с крайней осторожностью и, возможно, собрать дополнительные признаки для перерасчета.
- Средний запас сырья, основываясь на предсказаниях моделей:
- регион 1 - 92.4 тыс. баррелей,
- регион 2 - 68.71 тыс. баррелей,
- регион 3 - 94.77 тыс. баррелей.
- Средние значения запасов сырья, основываясь на предсказаниях моделей, близки к средним значениям, посчитанным на основе настоящих данных:
- регион 1 - 92.33 тыс. баррелей,
- регион 2 - 68.73 тыс. баррелей,
- регион 3 - 95.15 тыс. баррелей.
- В отношении кадой из выборок были проведены статистические тесты об однородности распределений по критерию Смирнова в результате которого, исходя из имеющихся данных, при уровне значимости 0.05, были отвегнуты нулевые гипотезы о принадлежности выборок предсказанных моделью объемов сырья и выборок правильных объемов сырья одному закону распределения. Принята альтерантивная гипотеза - выборка предсказанных моделью объемов сырья и выборка правильных объемов сырья в отношени регионов 1,2 и 3 принадлежат разным законам распределения.
- Исходя из полученных результатов, с осторожностью можно утверждать, что модели обученные на представленных данных - малопригодны для прогонозирования объемов запасов сырья с разумной степенью достоверности. Использование результатов работы модели для принятия решений является высокорискованным.
1) Выделим ключевые значения для расчетов¶
- общий бюджет на разработку месторождений - 10 млрд. рублей¶
In [80]:
budget = 10**10
budget
Out[80]:
- стоимость бурения одной скважины — 50 млн рублей¶
In [81]:
well_price = 50 * 10**6
well_price
Out[81]:
- прибыль от продажи барреля сырья - 4500 рублей¶
In [82]:
price_per_barrel = 4500
price_per_barrel
Out[82]:
2) Подсчитаем минимальный средний объём сырья в месторождениях региона, достаточный для его разработки¶
- подсчитаем предельное число скважин, исходя из общего бюджета¶
In [83]:
num_of_wells = int(budget / well_price)
num_of_wells
Out[83]:
- минимальный средний объём сырья на скважину, для этого разделим стоимость разоаботки на стоимость одного барреля¶
In [84]:
mininmum_mean_product = well_price/price_per_barrel
mininmum_mean_product
Out[84]:
⚡ Вывод:
- Предельное число скважин исходя из бюджета - 200.
- Минимальный средний объём сырья на скважину - 11 тыс. баррелей.
- Во всех представленных регионах подавляющее большинство меторождений превышают минимальный требующийся объем.
3) Подготовим функцию для расчёта прибыли по набору отобранных месторождений и предсказаний модели¶
- при расчете прибыли необходимо учитывать только те скважины, стоимость бурения которых окупиться, то есть с запасом свыше 11 тыс. баррелей.
- в представленных данных запасы месторождений таковы, что проверка на минимальный запас является излишней, но принимая во внимание принцип безконстантности: функция для расчета выручки будет включать такую проверку.
- Выделены ключевые значения для расчетов:
- бюджет: 10 млрд. рублей,
- стоимость бурения скважины - 50 млн. рублей,
- прибыль от продажи барреля нефти - 4500 рублей.
- Предельное число скважин, исходя из бюджета - 200 штук.
- Минимальный средний объём сырья на скважину 11111 баррелей.
- Подготовлена функция для расчета выручки, с учетом проверки достаточности запасов для его разработки.
- создадим список, в который будем вносить результаты¶
In [87]:
revenue_report = []
- регион 1¶
In [88]:
bootstrap_region_one = bootstrap(region_1_report['predicted'], 500)
- подсчитаем среднюю прибыль, 2.5%-квантиль, 95% доверительный интервал¶
In [89]:
mean_rev = bootstrap_region_one.mean()/10**9
lower = bootstrap_region_one.quantile(0.025)/10**9
confidence_interval = st.t.interval(0.95, len(bootstrap_region_one)-1,
mean_rev, bootstrap_region_one.sem()/10**9)
- визуализируем распредление, 2.5%-квантиль¶
In [90]:
ax = sns.distplot(bootstrap_region_one)
ax.axvline(np.quantile(bootstrap_region_one, 0.025),linestyle = '--',color='r')
Out[90]:
- дополним отчет¶
In [91]:
revenue_report.append(('region_1', mean_rev, lower, confidence_interval))
- регион 2¶
In [92]:
bootstrap_region_two = bootstrap(region_2_report['predicted'], 500)
- подсчитаем среднюю прибыль и 2.5%-квантиль, 95% доверительный интервал¶
In [93]:
mean_rev = bootstrap_region_two.mean()/10**9
lower = bootstrap_region_two.quantile(0.025)/10**9
confidence_interval = st.t.interval(0.95, len(bootstrap_region_two)-1,
mean_rev, bootstrap_region_two.sem()/10**9)
- визуализируем распредление, 2.5%-квантиль¶
In [94]:
ax = sns.distplot(bootstrap_region_one)
ax.axvline(np.quantile(bootstrap_region_two, 0.025),linestyle = '--',color='r')
Out[94]:
- дополним отчет¶
In [95]:
revenue_report.append(('region_2', mean_rev, lower, confidence_interval))
- регион 3¶
In [96]:
bootstrap_region_three = bootstrap(region_3_report['predicted'], 500)
- подсчитаем среднюю прибыль, 2.5%-квантиль, 95% доверительный интервал¶
In [97]:
mean_rev = bootstrap_region_three.mean()/10**9
lower = bootstrap_region_three.quantile(0.025)/10**9
confidence_interval = st.t.interval(0.95, len(bootstrap_region_three)-1,
mean_rev, bootstrap_region_three.sem()/10**9)
- визуализируем распредление, 2.5%-квантиль¶
In [98]:
ax = sns.distplot(bootstrap_region_three)
ax.axvline(np.quantile(bootstrap_region_three, 0.025),linestyle = '--',color='r')
Out[98]:
- дополним отчет¶
In [99]:
revenue_report.append(('region_3', mean_rev, lower, confidence_interval))
- выведем revenue_report, содержащий данные о средней прибыли, 95-% доверительный интервал и 2.5%-квантили для прибыли¶
In [100]:
revenue_report_df = pd.DataFrame(revenue_report, columns=['region', 'mean_revenue', '2.5%-quantile', '95%_conf_int'])
revenue_report_df.style.highlight_max(color='lightgreen', axis = 0)
Out[100]:
- Исходя из представленных данных, принимая во внимание результат исследования - наилучший регион для разработки: регион № 1, несмотря на то что лучшая средняя выручка у региона № 2.
- Hезультаты в отношении региона № 2 детерменированы харакетром представленных данных, поэтому не могут приниматься в расчет до валидации у поставщика данных.
- Принимая во внимание близость качества моделей, предсказавших запасы в регионе № 1 и 3: средняя выручка в регионе № 1 (103.3 млрд. рублей) больше выручки в регионе № 3 (102.69), равно как и 2.5%-квантиль (101.1 млрд. рублей в первом реионе против 100.77 млрд. рублей в регионе № 3).
Среднняя выручка по регионам:
- регион № 1: 103.3 млрд. рублей,
- регион № 2: 104.5 млрд. рублей,
- регион № 3: 102.69 млрд. рублей.
Исследование показало, что исходя из представленных данных, риск убытков отсутвует, 2.5%-квантиль для регионов:
- регион № 1: 101.1 млрд. рублей,
- регион № 2: 100.6 млрд. рублей,
- регион № 3: 100.8 млрд. рублей.
- Исходя из полученных результатов, с осторожностью можно утверждать, что модели обученные на представленных данных - малопригодны для прогонозирования объемов запасов сырья с разумной степенью достоверности. Использование результатов работы моделей для принятия решений является высокорискованным.