Описание проекта¶
Подготовьте прототип модели машинного обучения для «Цифры». Компания разрабатывает решения для эффективной работы промышленных предприятий.
Модель должна предсказать коэффициент восстановления золота из золотосодержащей руды. В вашем распоряжении данные с параметрами добычи и очистки.
Модель поможет оптимизировать производство, чтобы не запускать предприятие с убыточными характеристиками.
Вам нужно:
- Подготовить данные;
- Провести исследовательский анализ данных;
- Построить и обучить модель.
Чтобы выполнить проект, обращайтесь к библиотекам pandas, matplotlib и sklearn. Вам поможет их документация.
(... загрузка данных и библиотек ...)
1. Подготовка данных¶
Посмотрим какие признаки есть в трейне, но нет в тесте
set(train.columns) - set(test.columns)
В тесте у нас нет признаков output, так как это целевые признаки.
Так же в тесте нет признаков calculation для этапа rougher, видимо эти данные рассчитываются позднее и не доступны во время процесса.
a = [1,2,3,4,5]
b = [2,3,4,6,8]
set(b) - set(a)
Проверим расчет recovery в train:
- Вычислим recovery из сырых значений, проверим сколько из них совпало с предоставленными
- Вычислим MAE для рассчитанных значений и исходых
def recovery(C, F, T):
numerator = (C*(F-T))
denominator =(F*(C-T))
rec = numerator / denominator * 100
#nonzero_denominator = denominator != 0
#rec[~nonzero_denominator] = 0
# так как мы не застрахованы от очень больших и очень маленьких значений, то заполним их
rec[rec<0] = np.nan
rec[rec>100] = np.nan
return rec
t,f,c = train['rougher.output.tail_au'], train['rougher.input.feed_au'], train['rougher.output.concentrate_au']
rec = recovery(c, f, t)
right_recovery_sum = np.isclose(train['rougher.output.recovery'], rec).sum()
right_recovery_sum, train.shape[0]-right_recovery_sum
rec.describe()
Совпало 14287 рассчетных с иходными, не совпало 2573, возможно там пропуски.
Для расчета MAE, заполним пропуски нулями =)
mean_absolute_error(train['rougher.output.recovery'].interpolate(method='time'),
rec.interpolate(method='time'))
rec.isna().sum()
Как мы видим средняя абсолютная ошибка достаточно мала, считаем считаем что recovery посчитан верно
Посмотрим на пропуски в данных в трейне и тесте:
def show_na(df):
na_info = (df.isna() | df.isnull()).sum()
res = (pd.concat([na_info / df.shape[0], na_info], axis=1, keys=['percent', 'abs'])
.sort_values('percent', ascending=False))
return res
show_na(train).head(10)
show_na(test).head(10)
Пропуски в тесте не значительны
Пропуски в трейне более значитиельны, большая часть из них в целевых признаках .output., так же в rougher.output.recovery пропусков 2573 что совпадает с количеством несовпавших расчетных значений
Пока не будем обрабатывть пропуски, при машинном обучении на кросс-валидации проверим какой способ заполнения или удаления пропусков покажет себя лучше
Подготовка данных¶
Для подготовки данных к машинному обучению
- Удалим из train признаки которых нет в тесте
- Разделим данные не целевой признак и все остальные
- Так же разделим датасет на 2 части rougher и cleaner, так как процессы флотации и очистки не зависимы, то и предсказывать их тоже логично независимо.
Возможно предсказзывать recovery вообще не верный путь, так как он зависит от одного из обучающих параметров, может лучше для каждого процесса предсказыать output.concentrate и output.tail и затем на их основе уже считать recovery.
Проверим это на этапе машинного обучения, пока что просто запишем эти данные в отдельные переменные
X_train = train[test.columns]
rougher_cols = X_train.columns.str.contains('rougher')
X_train_rougher = X_train.loc[:, rougher_cols]
X_test_rougher = test.loc[:, rougher_cols]
y_train_rougher = train[['rougher.output.tail_au', 'rougher.output.concentrate_au']]
final_cols = X_train.columns.str.contains('cleaner')
X_train_cleaner = X_train.loc[:, final_cols]
X_test_cleaner = test.loc[:, final_cols]
y_train_cleaner = train[['final.output.tail_au', 'final.output.concentrate_au']]
2. Анализ данных¶
Построим распределения концентраций металов на выходе каждого процесса и в хвостах
process = ['rougher.input.feed',
'rougher.output.concentrate',
'primary_cleaner.output.concentrate',
'final.output.concentrate']
process_tail = ['rougher.input.feed',
'rougher.output.tail',
'primary_cleaner.output.tail',
'final.output.tail']
metals = ['au', 'ag', 'pb']
fig, axs = plt.subplots(1, len(process), figsize=(20, 6), constrained_layout=True)
fig.suptitle('Metal concentrations in output by stage', fontsize=24)
for stage, ax in zip(process, axs):
ax.set_title(stage)
for metal in metals:
cols = train.columns.str.contains(stage+'_'+metal)
sns_ax = sns.distplot(train.loc[:, cols].dropna(), label=metal, ax=ax)
plt.legend()
fig, axs = plt.subplots(1, 4, figsize=(20, 6), constrained_layout=True)
fig.suptitle('Metal concentrations in tail by stage', fontsize=24)
for stage, ax in zip(process_tail, axs):
ax.set_title(stage)
for metal in metals:
cols = train.columns.str.contains(stage+'_'+metal)
sns_ax = sns.distplot(train.loc[:, cols].dropna(), label=metal, ax=ax)
plt.legend()
plt.show()
Из графиков видно что в процессе очистики концентрация золота возрастает, отметим достаточно большое количетсво нулевых значений
Так же отмети различную кноцентрацию металлов в хвостовых отвалах в после различных этапов
Доролнительно посторим графики для концентрации кадого металла на разны процессах на выходе и в хвосте
fig, axs = plt.subplots(1, len(metals), figsize=(20, 6), constrained_layout=True)
fig.suptitle('Metal concentrations in output by metal', fontsize=24)
for metal, ax in zip(metals, axs):
ax.set_title(metal)
for stage in process:
cols = train.columns.str.contains(stage+'_'+metal)
sns_ax = sns.distplot(train.loc[:, cols].dropna(), label=stage, ax=ax)
plt.legend()
fig, axs = plt.subplots(1, len(metals), figsize=(20, 6), constrained_layout=True)
fig.suptitle('Metal concentrations in tail by metal', fontsize=24)
for metal, ax in zip(metals, axs):
ax.set_title(metal)
for stage in process_tail:
cols = train.columns.str.contains(stage+'_'+metal)
sns_ax = sns.distplot(train.loc[:, cols].dropna(), label=stage, ax=ax)
plt.legend()
plt.show()
Отметим, что концентрация золота по мере очистки значительно вырастает, так же вырастают концентрации серебра и свинца но не так значительно.
Самая большая концентрация холота в хвостах после флотации, при очиске концентрация золота св хвоствх чуть меньше. Тоже самое характерно для других металлов.
Проверим размер гранул сырья в трейне и в тесте
sns.distplot(train['rougher.input.feed_size'].dropna(), label='train')
sns.distplot(test['rougher.input.feed_size'].dropna(), label='test')
plt.legend()
Как мы видим распределения очень похожи
Исследуем суммарню концерацию в разных процессах в трейне и в тесте
PS не очень понимаю зачем это нужно, видимо что бы найти большое количество нулевых значений и что-то с ними сделать
fig, axs = plt.subplots(1, len(process), figsize=(20, 6), constrained_layout=True)
fig.suptitle('Summary metal concentrations in train and test by process', fontsize=24)
for stage, ax in zip(process, axs):
ax.set_title(stage)
sum_train = train[stage+ '_ag'] + train[stage+ '_au'] + train[stage+ '_pb']
try:
sum_test = test[stage+ '_ag'] + test[stage+ '_au'] + test[stage+ '_pb']
sns.distplot(sum_test.dropna(), label='test', ax=ax)
except KeyError:
pass
sns.distplot(sum_train.dropna(), label='train', ax=ax)
plt.legend()
plt.show()
3. Модель¶
Перейдем к машинному обучению, стратегия следующая:
- Готовим пайплайн: Заполнение пропусков -> Скалер(для регрессии) -> Модель
- Готовим параметры для поиска по сетке
- Заполним пропуски в целевой переменной и удалим нули из нее. только для train
- Делаем поиск по сетке, выбираем лучшую модель для флотации по метрике sMAPE
- То же самое для очистки
- Прогоняем пайплайн с лушими параметрами по всем данным, получаем предсказания для теста выходной и хвостовой концентраций
- Расчитываем recovery, для каждой стадидии.
- Подглядываем в full и смотрим конечное значение метрики
pipe = Pipeline([
('imp', SimpleImputer(missing_values=np.nan)),
('scaler', StandardScaler()),
('model', RandomForestRegressor(n_estimators=100, random_state=SEED))
])
params = [
{
'imp__strategy': ['mean', 'median', 'most_frequent'],
'model': [RandomForestRegressor(n_estimators=10, random_state=SEED)],
'model__max_features': np.linspace(0.1, 1, 10)
}, {
'imp__strategy': ['mean', 'median', 'most_frequent'],
'model': [LinearRegression()]
}, {
'imp__strategy': ['mean', 'median', 'most_frequent'],
'model': [Ridge(random_state=SEED)],
'model__alpha': np.logspace(-3, 1, 10)
}, {
'imp__strategy': ['mean', 'median', 'most_frequent'],
'model': [Lasso(random_state=SEED)],
'model__alpha': np.logspace(-3, 1, 10)
}
]
#переделал на заполнение интерполяцией
def fill_target_nan(y):
y = y.interpolate(method='time')
return y
def drop_target_zeros(X, y):
y = y[(y != 0).all(1)]
X = X.loc[y.index, :]
return X, y
print('Shapes before:')
print(X_train_rougher.shape, y_train_rougher.shape,
X_train_cleaner.shape, y_train_cleaner.shape)
y_train_rougher = fill_target_nan(y_train_rougher)
y_train_cleaner = fill_target_nan(y_train_cleaner)
X_train_rougher, y_train_rougher = drop_target_zeros(X_train_rougher, y_train_rougher)
X_train_cleaner, y_train_cleaner = drop_target_zeros(X_train_cleaner, y_train_cleaner)
print('Shapes after:')
print(X_train_rougher.shape, y_train_rougher.shape,
X_train_cleaner.shape, y_train_cleaner.shape)
#cv = KFold(n_splits=5, shuffle=False, random_state=SEED)
cv =TimeSeriesSplit(n_splits=3)
# возможно тут надо применить TimeSeriesSplit
# но из услови задачи не понятно, насколько процессы зависимы от времени и их длительность
grid_rougher = GridSearchCV(pipe, param_grid=params, cv=cv, scoring=neg_smape, n_jobs=-1)
temp = np.ones(train.shape[0])
plt.plot(pd.Series(data=temp, index=train.index), 'bo', c='blue', label='train index')
temp = np.zeros(test.shape[0])
plt.plot(pd.Series(data=temp, index=test.index), 'bo', c='red', label='test index')
plt.title('Train and test index')
plt.xticks(rotation=90)
plt.legend()
plt.show()
%%time
grid_rougher.fit(X_train_rougher, y_train_rougher)
grid_rougher.best_params_, grid_rougher.best_score_
grid_cleaner = GridSearchCV(pipe, param_grid=params, cv=cv, scoring=neg_smape, n_jobs=-1)
%%time
grid_cleaner.fit(X_train_cleaner, y_train_cleaner)
grid_cleaner.best_params_, grid_cleaner.best_score_
input_au = test['rougher.input.feed_au'].fillna(0)
pipe_rougher = grid_rougher.best_estimator_
if type(pipe_rougher.steps[2][1]) is type(RandomForestRegressor()):
pipe_rougher.steps[2][1].n_estimators = 100 #с увеличением числа деревьев качество должно улучшаться
pipe_rougher.fit(X_train_rougher, y_train_rougher)
y_pred_rougher_tail, y_pred_rougher_conc = pipe_rougher.predict(X_test_rougher).T
rougher_recovery = recovery(y_pred_rougher_conc, input_au, y_pred_rougher_tail)
smape_rougher = sMAPE(full.loc[X_test_rougher.index, 'rougher.output.recovery'].interpolate(method='time'),
rougher_recovery.interpolate(method='time'))
smape_rougher
pipe_cleaner = grid_cleaner.best_estimator_
if type(pipe_cleaner.steps[2][1]) is type(RandomForestRegressor()):
pipe_cleaner.steps[2][1].n_estimators = 100 #с увеличением числа деревьев качество должно улучшаться
pipe_cleaner.fit(X_train_cleaner, y_train_cleaner)
y_pred_cleaner_tail, y_pred_cleaner_conc = pipe_cleaner.predict(X_test_cleaner).T
cleaner_recovery = recovery(y_pred_cleaner_conc, input_au, y_pred_cleaner_tail)
smape_cleaner = sMAPE(full.loc[X_test_rougher.index, 'final.output.recovery'].fillna(0),
cleaner_recovery.fillna(0))
smape_cleaner
final_smape = 0.25*smape_rougher + 0.75*smape_cleaner
final_smape
PS Пришло еще несколько идей:
- Предсказывать для rougher все остальные output и подавать их как фичи в пайплайн для cleaner
- Возможно сделать отдельные модели для primari cleaner и second cleaner