В задаче необходимо предсказать, продлит пользователь сервиса подписку или нет. Задача бинарной классификации.
** Список переменных **
1 - id идентификатор.
2 - taxactionSystem система налогообложения (категориальный).
3 - regdt Дата регистрации (число).
4 - workerCount количество сотрудников (число).
5 - fssdccount количество отправленных отчетов в ФСС из этой организации из БК за все время существования (число).
6 - pfrdcCount количество отправленных отчетов в ПФР из этой организации из БК за все время существования (число).
7 - fnsdcCount количество отправленных отчетов в ФНС из этой организации из БК за все время существования (число).
8 - hasCloudCryptCertificate был ли когда-либо в этой организации выпущен облачный сертификат (бинарный).
9 - OrgCreationDate дата добавления организации в БК Это не дата регистрации организации в ФНС и т.д., это дата, когда организация была добавлена (создана) в БК (число).
10 - documentsCount количество документов. Считает количество документов в системе (которые показываются на вкладке ""Все"") в этом количестве учитываются не все документы ) (число).
11 - cnt_users количество пользователей (число).
Целевая переменная 12 - is_prolong - продлится пользователь или нет. (бинарная: 1, 0)
Формат файла ответов:
id, is_prolong
Импортируем бииблиотеки
# -*- coding: utf-8 -*-
import pandas as pd
import math
import re
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import classification_report, accuracy_score,roc_auc_score, roc_curve, confusion_matrix
from sklearn.preprocessing import LabelEncoder, OneHotEncoder,StandardScaler,LabelBinarizer
from sklearn.datasets import fetch_20newsgroups, load_files
from sklearn.model_selection import GridSearchCV,StratifiedKFold
from sklearn.ensemble import RandomForestClassifier,GradientBoostingClassifier
from sklearn.metrics import roc_auc_score,precision_recall_curve,roc_curve
from sklearn.pipeline import Pipeline
from sklearn.neighbors import KNeighborsClassifier
from sklearn.learning_curve import validation_curve,learning_curve
import xgboost as xgb
from evolutionary_search import EvolutionaryAlgorithmSearchCV
%matplotlib inline
#import sys
#reload(sys)
#sys.setdefaultencoding('utf-8')
# настройка внешнего вида графиков в seaborn
sns.set_style("dark")
sns.set_palette("RdBu")
sns.set_context("notebook", font_scale = 1.5,
rc = { "figure.figsize" : (15, 5), "axes.titlesize" : 18 })
Загрузим данные
train=pd.read_csv('../../data/prolongation_service_train.csv',sep='\t',encoding='cp1251',parse_dates=['regdt','OrgCreationDate'])
test=pd.read_csv('../../data/prolongation_service_test.csv',sep='\t',encoding='cp1251',parse_dates=['regdt','OrgCreationDate'])
train.head()
Выведем основные харакетристики переменных:
print(train.shape)
train.describe(include = "all").T
Визуализируем данные
sns.pairplot(train)
print 'Распределение целевой переменной: \t'
sns.countplot(train.is_prolong)
Имеется небольшой дисбаланс классов, избавимся от него дальше.
sns.heatmap(train.corr(method='pearson'),xticklabels=True,yticklabels=True);
sns.heatmap(train.corr(method='spearman'),xticklabels=True,yticklabels=True);
Корреляция признаков "count" обусловлена тем, что чем чаще клиент пользуется сервисом и отправляет отчетность, тем больше он захочет продлить сервис.
df = pd.concat([train,test],ignore_index=True)
df.info()
#Удалим столбец id, так как он не несет никакой информации
df.drop(columns=['id'],axis=1,inplace=True)
df.fillna(value=0,inplace=True)
#Заменим строчки с кривой датой
#Дата 0001-01-01 обусловлена default value в sql server, т.е при вставке данных в таблицу не была указана дата
df['regdt'] = df['regdt'].replace(['0001-01-01 00:00:00.0000000'],u'2013-07-16') # 2013-07-16 медиана
#приведем дату в нужный формат
df[['regdt','OrgCreationDate']]=df[['regdt','OrgCreationDate']].apply(pd.to_datetime)
df['regdt'] = df['regdt'].fillna(value=pd._libs.tslib.Timestamp('2013-07-16 00:00:00'))
#Приведем колонки к формату int
columns = df.select_dtypes(['floating']).columns
df[columns] = df[columns].astype('int64')
#Работаем с датой, достанем год, месяц, день
df['regdt_year']=pd.DatetimeIndex(df['regdt']).year
df['regdt_month']=pd.DatetimeIndex(df['regdt']).month
df['regdt_day']=pd.DatetimeIndex(df['regdt']).day
df['OrgCreationDate_year']=pd.DatetimeIndex(df['OrgCreationDate']).year.astype('int64')
df['OrgCreationDate_month']=pd.DatetimeIndex(df['OrgCreationDate']).month.astype('int64')
df['OrgCreationDate_day']=pd.DatetimeIndex(df['OrgCreationDate']).day.astype('int64')
df['delta_year']=df['regdt_year']-df['OrgCreationDate_year']
#Попробуем добавить кол-во лет с момента регистрации
df['today-regdt_year']=2018-df['regdt_year']
df['today-OrgCreationDate_year']=2018-df['OrgCreationDate_year']
df.drop(axis=1,columns=['regdt','OrgCreationDate'],inplace=True)
df.info()
Достанем из признака taxactionSystem систему налогооблажения и величину ставки
new_tax=pd.DataFrame(columns=['tax','stavka'])
for idx,i in enumerate(df.taxactionSystem):
temp=i.split(', ')
if len(temp)!=2:
new_tax.loc[idx]=[temp[0],0]
else:
new_tax.loc[idx]=[temp[0],int(str(re.search(r'\d+%', temp[1]).group(0))[:-1])]
df['tax']=new_tax['tax']
df['stavka']=new_tax['stavka']
df.drop(columns='taxactionSystem',axis=1,inplace=True)
#train.dropna(axis=0, how='any',inplace=True)
df['stavka'] = df['stavka'].astype('int64')
df.head()
label_encoder = LabelEncoder()
mapped_education = pd.Series(label_encoder.fit_transform(np.hstack((df['tax']))))
mapped_education.value_counts().plot.barh()
# integer encode
integer_encoded = label_encoder.fit_transform(df['tax']).astype('int')
onehot_encoder = OneHotEncoder(sparse=False)
integer_encoded = integer_encoded.reshape(len(integer_encoded), 1)
onehot_encoded = onehot_encoder.fit_transform(integer_encoded)
print(onehot_encoded[:5])
# invert first example
inverted = label_encoder.inverse_transform([np.argmax(onehot_encoded[0, :])])
df=df.join(pd.DataFrame(onehot_encoded,columns=["tax"+str(i) for i in range(onehot_encoded.shape[1])]),how='outer')
df.drop(columns='tax',axis=1,inplace=True)
#train.dropna(axis=0, how='any',inplace=True)
df[df.columns] = df[df.columns].astype('int64')
df.head()
print len(train), len(test)
train = df[:][:len(train)]
test = df[:][len(train):]
Разделим выборку на обучающую и тестовую
target=train.is_prolong
train.drop(axis=1,columns=['is_prolong'],inplace=True)
train_X,test_X,train_y,test_y=train_test_split(train,target,shuffle=True,random_state=17,test_size=0.2)
Выровняем классы с помощью SMOTE
from imblearn.over_sampling import SMOTE
smote = SMOTE(kind='borderline1', random_state=17)
train_X, train_y = smote.fit_sample(train_X, train_y)
train_X=pd.DataFrame(train_X,columns=train.columns)
sns.countplot(train_y)
В качестве метрики будем использовать ROC_AUC, так как имеем дело с несбалансированной выборкой, и порог, равный 0.5, может оказывается не оптимальным
Напишем функцию для получения метрики roc-auc
def getROC_AUC(method,test_X,test_y):
precision, recall, thresholds = precision_recall_curve(test_y, method.predict_proba(test_X)[:,1])
plt.title("precision-recall")
plt.plot(recall,precision)
fpr, tpr, thresholds = roc_curve(test_y, method.predict_proba(test_X)[:,1])
plt.figure()
plt.title("ROC-AUC")
plt.plot(fpr,tpr)
Попробуем применить LR с подбором параметров с помощью EvolutionaryAlgorithmSearchCV
lr=LogisticRegression(n_jobs=-1)
lr_pipe = Pipeline([('scaler', StandardScaler()), ('lr', LogisticRegression())])
paramgrid = {"lr__C": np.logspace(-8, 1, 10),
"lr__penalty" : ['l1','l2']}
np.random.seed(17)
cv = EvolutionaryAlgorithmSearchCV(estimator=lr_pipe,
params=paramgrid,
scoring="roc_auc",
cv=StratifiedKFold(n_splits=4),
verbose=1,
population_size=10,
gene_mutation_prob=0.10,
gene_crossover_prob=0.5,
tournament_size=3,
generations_number=10,
n_jobs=4)
cv.fit(train_X, train_y)
Построим модель с лучшими параметрами, и посмотрим качество на отложенной выборке
lr = LogisticRegression(n_jobs=-1,penalty='l1',C=1)
scaler_lr=StandardScaler()
train_X_scale=scaler_lr.fit_transform(train_X)
test_X_scale=scaler_lr.transform(test_X)
lr.fit(train_X_scale, train_y)
lr_pred=lr.predict(test_X_scale)
print accuracy_score(test_y, lr_pred)
print roc_auc_score(test_y, lr.predict_proba(test_X_scale)[:,-1])
Оценим важность признаков
Лассо обнуляет веса ненужных признаков
features = pd.DataFrame(lr.coef_.reshape((20,1)),
index=train_X.columns,
columns=['Importance']).sort_values(['Importance'], ascending=False)
features
Видим, что наибольший вклад внёс признак hasCloudCryptCertificate. Компании, у которых имеется электронаая подпись, чаще продливают подписку, это подтверждено практикой.
Сделаем тоже самое для Случайного леса
rnd_forest=RandomForestClassifier(n_jobs=-1,random_state=17,oob_score=False)
forest_params = {'max_depth': [5,6,7,8,9,10]
,'max_features': [1,2,3,4,5,6,7]
,'min_samples_leaf':[1,2,3,5,7,8]
,'n_estimators':[70,80,90,100,110,125,135,150]
}
np.random.seed(17)
cv = EvolutionaryAlgorithmSearchCV(estimator=rnd_forest,
params=forest_params,
scoring="roc_auc",
cv=5,
verbose=1,
population_size=100,
gene_mutation_prob=0.10,
gene_crossover_prob=0.5,
tournament_size=3,
generations_number=10,
n_jobs=8)
cv.fit(train_X, train_y)
Построим лучшую модель
rnd_forest=RandomForestClassifier(n_jobs=-1,n_estimators=150,max_depth=10,max_features=5,min_samples_leaf=1,random_state=17)
rnd_forest.fit(train_X,train_y)
pred_forest=rnd_forest.predict(test_X)
print accuracy_score(test_y,pred_forest)
print roc_auc_score(test_y, rnd_forest.predict_proba(test_X)[:,-1])
features = pd.DataFrame(rnd_forest.feature_importances_,
index=train_X.columns,
columns=['Importance']).sort_values(['Importance'], ascending=False)
features
Случайный лес, также в качестве наиболее важного признака выбрал hasCloudCryptCertificate
График "важности" признаков
plt.plot(range(len(features.Importance.tolist())),
features.Importance.tolist())
def plot_with_std(x, data, **kwargs):
mu, std = data.mean(1), data.std(1)
lines = plt.plot(x, mu, '-', **kwargs)
plt.fill_between(x, mu - std, mu + std, edgecolor='none',
facecolor=lines[0].get_color(), alpha=0.2)
def plot_learning_curve(clf, X, y, scoring, cv=5):
train_sizes = np.linspace(0.05, 1, 20)
n_train, val_train, val_test = learning_curve(clf,
X, y, train_sizes, cv=cv,
scoring=scoring)
plot_with_std(n_train, val_train, label='training scores', c='green')
plot_with_std(n_train, val_test, label='validation scores', c='red')
plt.xlabel('Training Set Size'); plt.ylabel(scoring)
plt.legend()
plot_learning_curve(RandomForestClassifier(n_jobs=-1,n_estimators=150,
max_depth=10,max_features=5,
min_samples_leaf=1,random_state=17),
train_X, train_y, scoring=None, cv=5)
Построим валидационную кривую для данных параметров леса. В качестве параметра сложности будем использовать max_depth:
def plot_validation_curve(clf, X, y, cv_param_name,
cv_param_values, scoring):
val_train, val_test = validation_curve(clf, X, y, cv_param_name,
cv_param_values, cv=5,
scoring=scoring)
plot_with_std(cv_param_values, val_train,
label='training scores', c='green')
plot_with_std(cv_param_values, val_test,
label='validation scores', c='red')
plt.xlabel(cv_param_name); plt.ylabel(scoring)
plt.legend()
max_depth = range(1,10)
plot_validation_curve(RandomForestClassifier(n_jobs=-1,n_estimators=150,
max_depth=10,max_features=5,
min_samples_leaf=1,random_state=17), train_X, train_y,
cv_param_name='max_depth',
cv_param_values=max_depth,
scoring=None)
Найдём лучшие параметры для ещё нескольких алгоритмов и попробуем EnsembleClassifier
Попробуем градиентный бустинг над деревьями
g_boost=GradientBoostingClassifier(random_state=17,learning_rate=0.6)
boost_params = {'max_depth': range(1,10)
,'max_features': range(1,13)
,'min_samples_leaf':range(1,5)
,'n_estimators':range(10,150,10)
,'learning_rate':np.arange(0.01,1,0.05)
}
np.random.seed(17)
cv = EvolutionaryAlgorithmSearchCV(estimator=g_boost,
params=boost_params,
scoring="roc_auc",
cv=5,
verbose=1,
population_size=50,
gene_mutation_prob=0.10,
gene_crossover_prob=0.5,
tournament_size=3,
generations_number=10,
n_jobs=8)
cv.fit(X_train, train_y)
Построим лучшую модель
g_boost=GradientBoostingClassifier(random_state=17,n_estimators=120,max_depth=4,learning_rate=0.16,max_features=5,min_samples_leaf=4)
g_boost.fit(train_X,train_y)
pred_boost=g_boost.predict(test_X)
print accuracy_score(test_y,pred_boost)
print roc_auc_score(test_y, g_boost.predict_proba(test_X)[:,-1])
# x_boost=xgb.XGBClassifier()
# X_boost_params = {'max_depth': range(1,10)
# ,'n_estimators':range(10,150,10)
# ,'learning_rate':np.arange(0.01,0.2,0.05)
# }
# np.random.seed(17)
# cv = EvolutionaryAlgorithmSearchCV(estimator=x_boost,
# params=X_boost_params,
# scoring="accuracy",
# cv=5,
# verbose=1,
# population_size=50,
# gene_mutation_prob=0.10,
# gene_crossover_prob=0.5,
# tournament_size=3,
# generations_number=10,
# n_jobs=8)
# cv.fit(X_train, train_y)
Построим лучшую модель
# x_boost=xgb.XGBClassifier(n_estimators=70,learning_rate=0.06,max_depth=4)
# x_boost.fit(X_train,train_y)
# pred_x_boost=x_boost.predict(X_test)
# print accuracy_score(test_y,pred_x_boost)
# print roc_auc_score(test_y, x_boost.predict_proba(X_test)[:,-1])
knn_pipe = Pipeline([('scaler', StandardScaler()), ('knn', KNeighborsClassifier(n_jobs=-1,weights='distance',n_neighbors=14))])
n_neighbors=range(2, 15)
knn_params = {'knn__n_neighbors': n_neighbors}
knn_grid = GridSearchCV(knn_pipe, knn_params,
cv=5, n_jobs=-1,
verbose=True)
knn_grid.fit(train_X, train_y)
print (knn_grid.best_params_, knn_grid.best_score_)
print accuracy_score(test_y, knn_grid.predict(test_X))
print roc_auc_score(test_y, knn_grid.predict_proba(test_X)[:,-1])
from brew.base import Ensemble, EnsembleClassifier
from brew.stacking.stacker import EnsembleStack, EnsembleStackClassifier
from brew.combination.combiner import Combiner
from sklearn import clone
# Initializing Classifiers
clf1 = Pipeline([('scaler', StandardScaler()), ('knn', KNeighborsClassifier(n_jobs=-1,weights='distance',n_neighbors=14))])
clf2 = RandomForestClassifier(n_jobs=-1,n_estimators=150,max_depth=10,max_features=5,min_samples_leaf=1,random_state=17)
clf3 = GradientBoostingClassifier(random_state=17,n_estimators=110,max_depth=7,learning_rate=0.11,max_features=3,min_samples_leaf=3)
clf4 = Pipeline([('scaler', StandardScaler()), ('lr', LogisticRegression(n_jobs=-1,penalty='l2',C=10,random_state=17))])
clf5 = RandomForestClassifier(n_jobs=-1,n_estimators=100,max_depth=10,random_state=17)
# Creating Ensemble
ensemble = Ensemble([clf1, clf2, clf3, clf4, clf5])
eclf = EnsembleClassifier(ensemble=ensemble, combiner=Combiner('mean'))
# Creating Stacking
layer_1 = Ensemble([clf1, clf2, clf3, clf4])
layer_2 = Ensemble([clone(clf5)])
stack = EnsembleStack(cv=3)
stack.add_layer(layer_1)
stack.add_layer(layer_2)
sclf = EnsembleStackClassifier(stack)
# Loading some example data
X = X_train
y = train_y
d = {yi : i for i, yi in enumerate(set(y))}
y = np.array([d[yi] for yi in y])
sclf.fit(X,y)
print 'sclf ' + str(roc_auc_score(test_y,sclf.predict_proba(X_test)[:,-1]))
eclf.fit(X,y)
print 'eclf ' + str(roc_auc_score(test_y,eclf.predict_proba(X_test)[:,-1]))
getROC_AUC(sclf,X_test,test_y)
Результаты ну прям не очень, попробуем что-нибудь по-серьезнее:)
train_X,val_X,train_y,val_y=train_test_split(train,target,shuffle=True,random_state=17,test_size=0.2)
train_X1,train_X2,train_y1,train_y2=train_test_split(train_X,train_y,random_state=17,test_size=0.7)
data = pd.DataFrame()
lr=LogisticRegression(random_state=17)
lr = Pipeline([('scaler', StandardScaler()), ('lr', LogisticRegression(n_jobs=-1,penalty='l2',C=10,random_state=17))])
lr.fit(train_X1,train_y1)
data['lr'] = lr.predict_proba(train_X2)[:,-1]
random_forest = RandomForestClassifier(n_jobs=-1,n_estimators=70,max_depth=10,max_features=4,min_samples_leaf=3,random_state=17)
random_forest.fit(train_X1,train_y1)
data['random_forest'] = random_forest.predict_proba(train_X2)[:,-1]
gradient_boosting = GradientBoostingClassifier(random_state=17,n_estimators=140,max_depth=5,learning_rate=0.06,max_features=6,min_samples_leaf=4)
gradient_boosting.fit(train_X1,train_y1)
data['gradient_boosting'] = gradient_boosting.predict_proba(train_X2)[:,-1]
x_boost=xgb.XGBClassifier(n_estimators=70,learning_rate=0.06,max_depth=4)
x_boost.fit(train_X1,train_y1)
data['x_boost']=x_boost.predict_proba(train_X2)[:,-1]
kneighbors = Pipeline([('scaler', StandardScaler()), ('knn', KNeighborsClassifier(n_jobs=-1,weights='distance',n_neighbors=14))])
kneighbors.fit(train_X1,train_y1)
data['kneighbors']=kneighbors.predict_proba(train_X2)[:,-1]
data.head()
Функция, возвращающая предсказанные данные:
def get_predictions(models,train_X2):
data = pd.DataFrame()
for idx, i in enumerate(models):
data[str(idx)] = i.predict_proba(train_X2)[:,-1]
return data
train_X2 = get_predictions([lr,random_forest,gradient_boosting,x_boost,kneighbors], train_X2)
val_X = get_predictions([lr,random_forest,gradient_boosting,x_boost,kneighbors], val_X)
Обучим простую логистическую регрессию:
lr2 = LogisticRegression(C=0.5)
lr2.fit(train_X2,train_y2)
print accuracy_score(val_y,lr2.predict(val_X))
print roc_auc_score(val_y,lr2.predict_proba(val_X)[:,-1])
Результат тоже не впечатляет.
Наибольшее значение ROC_AUC у случайного леса, поэтому применим его для формирования итогового прогноза.
rnd_forest=RandomForestClassifier(n_jobs=-1,n_estimators=150,max_depth=10,max_features=5,min_samples_leaf=1,random_state=17)
rnd_forest.fit(train,target)
test.drop(axis=1,columns='is_prolong', inplace=True)
pd.DataFrame(np.stack((range(len(test)),rnd_forest.predict(test)), axis=1),columns=['id','is_prolong']).to_csv('submission.csv', index=False)
pd.read_csv('submission.csv').head()
К сожалению, оценить наш прогноз не получится, так как мне ещё не предоставили ответы. Будем надеяться, что оценка на тестовой выборке, будет близка к оценке на кросс валидации.