PESC - Programa de Engenharia de Sistemas e Computação
COPPE / UFRJ
Autores: Bernardo Souza e Rafael Lopes Conde dos Reis
GitHub: https://github.com/condereis/credit-card-default/
In [9]:
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
from scipy.stats import randint, uniform
from sklearn.decomposition import PCA
from sklearn.ensemble import GradientBoostingClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import cross_val_score, train_test_split, GridSearchCV, RandomizedSearchCV
from sklearn.metrics import accuracy_score, auc, confusion_matrix, f1_score, roc_curve
from sklearn.naive_bayes import GaussianNB
from sklearn.neighbors import KNeighborsClassifier
from sklearn.neural_network import MLPClassifier
from sklearn.svm import SVC
import tensorflow as tf
import tflearn
from tflearn.data_utils import to_categorical
%matplotlib inline
Este projeto utiliza a base de dados "Default of Credit Card Clients Dataset" disponibilizada pela UCI Machine Learning no site do Kaggle. "Default of Credit Card Clients Dataset" apresenta as informações de clientes de cartão de crédito e de seus pagamentos de fatura de cartão. Nela os clientes são clssificados em clientes que pagam uma fatura e clientes que atrazam ou dão calote em uma fatura do cartão de crédito.A cada entrada da base apresenta 25 variáveis. Abaixo é transcrita a descrição dada no Kaggle.
ID: ID of each client
LIMIT_BAL: Amount of given credit in NT dollars (includes individual and family/supplementary credit
SEX: Gender (1=male, 2=female)
EDUCATION: (1=graduate school, 2=university, 3=high school, 4=others, 5=unknown, 6=unknown)
MARRIAGE: Marital status (1=married, 2=single, 3=others)
AGE: Age in years
PAY_0: Repayment status in September, 2005 (-1=pay duly, 1=payment delay for one month, 2=payment delay for two months, ... 8=payment delay for eight months, 9=payment delay for nine months and above)
PAY_2: Repayment status in August, 2005 (scale same as above)
PAY_3: Repayment status in July, 2005 (scale same as above)
PAY_4: Repayment status in June, 2005 (scale same as above)
PAY_5: Repayment status in May, 2005 (scale same as above)
PAY_6: Repayment status in April, 2005 (scale same as above)
BILL_AMT1: Amount of bill statement in September, 2005 (NT dollar)
BILL_AMT2: Amount of bill statement in August, 2005 (NT dollar)
BILL_AMT3: Amount of bill statement in July, 2005 (NT dollar)
BILL_AMT4: Amount of bill statement in June, 2005 (NT dollar)
BILL_AMT5: Amount of bill statement in May, 2005 (NT dollar)
BILL_AMT6: Amount of bill statement in April, 2005 (NT dollar)
PAY_AMT1: Amount of previous payment in September, 2005 (NT dollar)
PAY_AMT2: Amount of previous payment in August, 2005 (NT dollar)
PAY_AMT3: Amount of previous payment in July, 2005 (NT dollar)
PAY_AMT4: Amount of previous payment in June, 2005 (NT dollar)
PAY_AMT5: Amount of previous payment in May, 2005 (NT dollar)
PAY_AMT6: Amount of previous payment in April, 2005 (NT dollar)
default.payment.next.month: Default payment (1=yes, 0=no)
In [10]:
train = pd.read_csv('../data/interim/train.csv', index_col=0)
del train.index.name
train.describe()
Out[10]:
In [3]:
sns.pairplot(train.sample(frac=0.1), diag_kind="hist", hue='default.payment.next.month')
Out[3]:
O pair plot acima mostra o plot de cada um dos atributos em relação ao outro e o histograma de cada atributo (isso ocm as classes ressaltadas nas cores). Pode-se observar que as classes apresentam grande sobreposição de valores, não havendo um valor óbvio entre quaisquer dois atributos que separe bem (mesmo que um subconjunto) os dados. Observou-se entretanto que os valores de classe 0 (não atraza pagamento) são mais esparços (apresentam mais outliers e pontos no mais distantes), isso não ajuda o problema de classificar clientes que atrazarão, mas permitiria identificar um subgrupo de clientes facilmente marcados como que não atrazam (o q poderia ser útil para o banco).
In [28]:
corrmat = train.corr()
sns.heatmap(corrmat, square=True)
Out[28]:
O mapa de calor acima exibe a correlação entre as diversas variáveis da base. Podemos observar uma correlação forte, devido a dependencia tempora, entre o atraso no pagamento nos últimos 6 meses (PAY_0 a PAY_6), entre as variáveis das faturas do cartão nos últimos 6 meses (BILL_AMT1 a BILL_AMT6) e entre os pagamentos feitos nos últimos 6 meses (PAY_AMT1 a PAY_AMT6). Existe também forte correlação entre alguns desses grupos, como entre o atraso e a fatura e a fatura e o pagamento. Entre as variáveis casamento e idade também há forte correlação.
As variáveis mais correlacionadas com o pagamento ou não da próxima fatura são, respectivamente, o atraso no pagamento dos últimos meses e o limite do cartão.
É possível observar, a partir dos plots abaixo, que a distribuição dos casos onde foi efetuado o pagamento possui um valor de crédito mais alto do que naqueles onde ouve default. Pode-se cogitar que valores mais altos de limite estão associados com menor risco de default no pagamento, uma vez que usualmente o lumite do cartão aumenta uma vez que o cliente se mostre comprometido com os pagamentos.
O primeiro é um boxplot, no qual podemos ver que a linha da mediana no caso de não pagamento se encontra abaixo da mediana no caso de pagamento. Isso é confirmado, a seguir pelos histogramas, onde o histograma de default se encontra acima do de pagamento para valores mais baixos de emprestimo e o oposto para valores mais altos. O ultimo mostra a taxa de não pagamento em cada quartil (separado por idade) de empréstimo. Este último mostra mais explicitamente a dependencia do valor do empréstimo na taxa de não pagamento.
In [5]:
sns.boxplot(x='default.payment.next.month', y='LIMIT_BAL', data=train);
In [6]:
sns.kdeplot(train[train['default.payment.next.month']==0].LIMIT_BAL, label='Pagamento')
sns.kdeplot(train[train['default.payment.next.month']==1].LIMIT_BAL, label='Default')
Out[6]:
In [7]:
categorical_limit_bal = pd.cut(train.LIMIT_BAL, [0,50000,140000,240000,1000000], labels=["Q1","Q2","Q3","Q4"])
sns.barplot(x=categorical_limit_bal, y=train['default.payment.next.month']);
Por fim, converteu-se o valor de limite do cartão, de dólares de Taiwan, para Reais (cotação do dia 04/12/2016 às 13:56), a fim de obeter uma noção melhoror do poder aquisitivo da base de usuarios. A descrição estatística do limite do cartão em reais mostra que apenas em torno de 1/4 das pessoas tem um limite entre, R\$1000 e R\$5500, ou seja, 2/3 da base tem limites superiores a R\$5500. Isso pode parecer um limite muito alto, mas se for normalizado pela relação entre o salário médio no Brasil e em Taiwan, obtemos a sequinte distribuição estatística, muito mais próxima do que se esperaria de uma base que contenha todas as faixas de poder aquisitivo.
In [5]:
(train.LIMIT_BAL/(9.17696 * 2.9975)).describe()
Out[5]:
É possível observar algo estranho nos dados PAY_x. Segundo a descrição do problema ela deveria ter apenas valores maiores que 1 referentes ao atraso e caso não tivesse havido atraso o valor da variável seria -1. Podemos observar os valores -2 e 0 nos dados, que, a princípio não eram esperados.
Como foi observado no gráfico de correlação entre variáveis notamos que as variáveis PAY_0 - PAY_6 são fortemente correlacionadas entre sí, portando analisamos abaixo apenas PAY_0 que, dentre elas, é a mais correlacionada com a saída.
Observamos que os valores -2 a 0 tem uma taxa de default muito menor do que dos demais valores. Como era possível esperar pessoas que com um grande atraso no pagamento tem uma taxa muita maior de default no próximo mês do que aquelas que tem os pagamentos em dia. A distribuição para os valores de atraso mostra que 0 é o valor mais frequente.
In [8]:
sns.barplot(x='PAY_0', y='default.payment.next.month', data=train);
In [9]:
sns.countplot(x=train.PAY_0);
Observamos no plot abaixo que university, graduate e high school são as classes mais representativas, respectivamente. Já o plote sequinte mostra que o valor médio da ocorrência de pagamento para cada categoria de escolaridade. Podemos observar que a taxa de ocorrencia de default é menor quento maior o grau de escolaridade.
In [10]:
sns.countplot(train.EDUCATION)
print '(0=unknown, 1=graduate school, 2=university, 3=high school, 4=others, 5=unknown, 6=unknown)'
In [11]:
sns.barplot(x='EDUCATION', y='default.payment.next.month', data=train);
print '(0=unknown, 1=graduate school, 2=university, 3=high school, 4=others, 5=unknown, 6=unknown)'
In [14]:
# sns.lmplot(x="size", y="default.payment.next.month", data=train);
sns.kdeplot(train[train['default.payment.next.month']==0].AGE, label='Pagamento')
sns.kdeplot(train[train['default.payment.next.month']==1].AGE, label='Default')
Out[14]:
In [18]:
categorical_age = pd.cut(train.AGE, [21, 28, 34, 41, 79], labels=["Q1","Q2","Q3","Q4"])
sns.barplot(x=categorical_age, y=train['default.payment.next.month']);
In [19]:
sex = train.SEX.replace(1,'Masculino').replace(2,'Feminino')
sns.barplot(x=sex, y=train['default.payment.next.month'])
Out[19]:
Pode-se observar uma diferença entre a taxa de pagamento de não casados (solteiros e outros) e casados, especialemnte ao sequmentar os dados nos quartis de faixa etária. Jovens (primeiro quartil) solteiros tem uma taxa de default média aproximadamente 8 pontos percentuais mais alta do que a de casados entre 28 e 34 anos (segundo quartil).
In [13]:
marrige = train.MARRIAGE.replace(0,u'Não casado')
marrige = marrige.replace(1,u'Não casado')
marrige = marrige.replace(3,u'Não casado')
marrige = marrige.replace(2,'Casado')
sns.barplot(x=marrige, y=train['default.payment.next.month'])
Out[13]:
In [20]:
sns.barplot(x=marrige, y=train['default.payment.next.month'], hue=categorical_age)
Out[20]:
In [19]:
corrmat = train[['PAY_AMT1', 'PAY_AMT2', 'PAY_AMT3', 'PAY_AMT4', 'PAY_AMT5', 'PAY_AMT6', 'BILL_AMT1', 'BILL_AMT2', 'BILL_AMT3', 'BILL_AMT4', 'BILL_AMT5', 'BILL_AMT6']].corr()
sns.heatmap(corrmat, vmax=.8, square=True)
Out[19]:
Tentou-se subtrair o valor pago do valor da fatura, a fim de reduzir o número de variáveis. Essas passam a indicar o valor que ficou em dívida da fatura do mês. Após esta redução, pode-se observar, que as variáveis encontradas tem uma correlação temporal ainda mais forte.
In [20]:
train['DEBT1'] = train.PAY_AMT1 - train.BILL_AMT2
train['DEBT2'] = train.PAY_AMT2 - train.BILL_AMT3
train['DEBT3'] = train.PAY_AMT3 - train.BILL_AMT4
train['DEBT4'] = train.PAY_AMT4 - train.BILL_AMT5
train['DEBT5'] = train.PAY_AMT5 - train.BILL_AMT6
# train['PAY_AMT6'] = train.PAY_AMT6 - train.BILL_AMT7
corrmat = train[['DEBT1', 'DEBT2', 'DEBT3', 'DEBT4', 'DEBT5']].sample(frac=0.1).corr()
sns.heatmap(corrmat, square=True)
Out[20]:
Por fim o conjunto das 6 variáveis de fatura e 6 de pagamento foram reduzidas para apenas uma, o saldo médio dos últimos meses. Ao relaciona-la com o limite do cartão obtemos as seguintes distribuições. Observamos que aqueles que pagaram a fatura seguinte tem um saldo médio, mais concentrado em 0. A concentração próxima a -1, ou seja um saldo negativo igual ao limite do cartão, é maior para o caso de default.
In [21]:
a = (train.PAY_AMT1+train.PAY_AMT2+train.PAY_AMT3+train.PAY_AMT4+train.PAY_AMT5+train.PAY_AMT6)
b = (train.BILL_AMT2+train.BILL_AMT3+train.BILL_AMT4+train.BILL_AMT5+train.BILL_AMT6)
train['percent_debt'] = (a-b)/train.LIMIT_BAL
sns.violinplot(y=train['percent_debt'], x=train['default.payment.next.month'])
Out[21]:
Os dados foram pré-processados de forma a:
In [2]:
train = pd.read_csv('../data/processed/train.csv', index_col=0)
del train.index.name
X = train.drop('default.payment.next.month', axis=1)
y = train['default.payment.next.month']
train.head()
Out[2]:
A acurácia de escolher a classe mais frequente é de 77,92%. Isso mostra que esse é um problema de classes desbalanceadas, de forma que a acurácia não é a melhor métrica para ser utilizada. Escolheu-se usar, então, como métrica a AUC (area under the curve), que representa a área sobre a curva ROC. A curva ROC, nada mais é do que uma curva, onde cada ponto corresponde a probabilidade de detecção e a probabilidade de falso alarme para um determinado valor de limiar de classificação.
In [74]:
print 'Acurácia:', max(y.mean(), 1-y.mean())
sns.countplot(y)
Out[74]:
In [75]:
model = GaussianNB()
scores = cross_val_score(model, X, y, cv=5, n_jobs=-1, scoring='roc_auc')
print 'Naive Bayes (AUC): %f +- %f' % (np.mean(scores), np.std(scores))
scores = cross_val_score(model, X, y, cv=5, n_jobs=-1)
print 'Naive Bayes (acurácia): %f +- %f' % (np.mean(scores), np.std(scores))
In [5]:
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.4)
model = GaussianNB()
y_pred = model.fit(X_train,y_train).predict(X_test)
labels = ['Pago', 'Default']
mtx = confusion_matrix(y_test, y_pred)
mtx = [x/float(sum(x)) for x in mtx]
sns.heatmap(pd.DataFrame(mtx, columns=labels, index=labels), annot=True, fmt=".2f", linewidths=.5)
plt.xlabel('Valor Classificado')
plt.ylabel('Valor Real')
Out[5]:
In [6]:
y_pred = model.fit(X_train, y_train).predict_proba(X_test)
y_pred = np.array([x[1] for x in y_pred])
fpr, tpr, thresholds = roc_curve(y_test+1, y_pred, pos_label=2)
roc_auc = auc(fpr, tpr)
plt.figure()
lw = 2
plt.plot(fpr, tpr, color='darkorange',
lw=lw, label='ROC curve (area = %0.2f)' % roc_auc)
plt.plot([0, 1], [0, 1], color='navy', lw=lw, linestyle='--')
plt.xlim([0.0, 1.0])
plt.ylim([0.0, 1.05])
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('Receiver operating characteristic example')
plt.legend(loc="lower right")
plt.show()
In [6]:
score_list = []
error_list = []
n_dims = range(1,train.shape[1])
model = KNeighborsClassifier()
for dim in n_dims:
reduced_data = PCA(n_components=dim).fit_transform(X)
scores = cross_val_score(model, reduced_data, y,
cv=5, n_jobs=-1, scoring='roc_auc')
score_list.append(np.mean(scores))
error_list.append(np.std(scores))
plt.errorbar(n_dims, score_list, yerr=error_list)
Out[6]:
In [10]:
X_red = CA(n_components=10).fit_transform(X)
Optou-se por utilizar um numero de vizinhos em torno do joelho da curva, que ocorreu em torno de 50. Com uma variação mais fina, chegou-se ao valor de 62 vizinhos.
In [16]:
k_list = [1, 10, 10, 50, 100, 150, 200]
score_list = []
error_list = []
for k in k_list:
model = KNeighborsClassifier(n_neighbors=k)
scores = cross_val_score(model, X_red, y, cv=5, n_jobs=-1, scoring='roc_auc')
score_list.append(np.mean(scores))
error_list.append(np.std(scores))
plt.errorbar(k_list, score_list, yerr=error_list)
Out[16]:
In [20]:
k_list = range(45,70,2)
score_list = []
error_list = []
for k in k_list:
model = KNeighborsClassifier(n_neighbors=k)
scores = cross_val_score(model, X_red, y, cv=5, n_jobs=-1, scoring='roc_auc')
score_list.append(np.mean(scores))
error_list.append(np.std(scores))
plt.errorbar(k_list, score_list, yerr=error_list)
Out[20]:
In [21]:
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.4)
model = KNeighborsClassifier(n_neighbors=62)
y_pred = model.fit(X_train, y_train).predict_proba(X_test)
y_pred = np.array([x[1] for x in y_pred])
fpr, tpr, thresholds = roc_curve(y_test+1, y_pred, pos_label=2)
roc_auc = auc(fpr, tpr)
plt.figure()
lw = 2
plt.plot(fpr, tpr, color='darkorange',
lw=lw, label='ROC curve (area = %0.2f)' % roc_auc)
plt.plot([0, 1], [0, 1], color='navy', lw=lw, linestyle='--')
plt.xlim([0.0, 1.0])
plt.ylim([0.0, 1.05])
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('Receiver operating characteristic example')
plt.legend(loc="lower right")
plt.show()
In [27]:
score_list = []
error_list = []
n_dims = range(1,train.shape[1])
model = LogisticRegression()
for dim in n_dims:
reduced_data = PCA(n_components=dim).fit_transform(X)
scores = cross_val_score(model, reduced_data, y,
cv=5, n_jobs=-1, scoring='roc_auc')
score_list.append(np.mean(scores))
error_list.append(np.std(scores))
plt.errorbar(n_dims, score_list, yerr=error_list)
Out[27]:
O paramentros de regularização C que escolheu-se utilizar foi C = 1, por ser um dos valores que resultaram no maior AUC e por estar próximo do joelho da curva.
In [36]:
c_list = [0.0001, 0.001, 0.01, 0.1, 1, 10, 100, 1000, 10000]
score_list = []
error_list = []
for c in c_list:
model = LogisticRegression(C=c)
scores = cross_val_score(model, X, y, cv=5, n_jobs=-1, scoring='roc_auc')
score_list.append(np.mean(scores))
error_list.append(np.std(scores))
plt.errorbar(c_list, score_list, yerr=error_list)
plt.xscale('log')
In [40]:
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.4)
model = LogisticRegression(C=0.1)
y_pred = model.fit(X_train, y_train).predict_proba(X_test)
y_pred = np.array([x[1] for x in y_pred])
fpr, tpr, thresholds = roc_curve(y_test+1, y_pred, pos_label=2)
roc_auc = auc(fpr, tpr)
plt.figure()
lw = 2
plt.plot(fpr, tpr, color='darkorange',
lw=lw, label='ROC curve (area = %0.2f)' % roc_auc)
plt.plot([0, 1], [0, 1], color='navy', lw=lw, linestyle='--')
plt.xlim([0.0, 1.0])
plt.ylim([0.0, 1.05])
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('Receiver operating characteristic example')
plt.legend(loc="lower right")
plt.show()
Pode-se observar que a redução do número de dimensões com PCA não gerou uma melhora do AUC. Entretretanto o score de AUC para 15 dimensões é bem próximo do valor máximo encontrado e poderia ser usado para gerar um modelo mais leve computacionalmente.
In [15]:
score_list = []
error_list = []
n_dims = range(1,train.shape[1])
model = GradientBoostingClassifier()
for dim in n_dims:
reduced_data = PCA(n_components=dim).fit_transform(X)
scores = cross_val_score(model, reduced_data, y,
cv=5, n_jobs=-1, scoring='roc_auc')
score_list.append(np.mean(scores))
error_list.append(np.std(scores))
plt.errorbar(n_dims, score_list, yerr=error_list)
Out[15]:
In [52]:
score_list = []
error_list = []
learning_rates = np.logspace(-4,1,20)
for learning_rate in learning_rates:
model = GradientBoostingClassifier(learning_rate=learning_rate)
scores = cross_val_score(model, X, y, cv=5, n_jobs=-1, scoring='roc_auc')
score_list.append(np.mean(scores))
error_list.append(np.std(scores))
plt.errorbar(learning_rates, score_list, yerr=error_list)
plt.xscale('log')
Corresponde a porcentagem de amostras que são utilizadas para fitar os classificadores do ensemble. Valores menores que 1 tendem a reduzir a variância porém aumentar o bias do modelo. Não se observou melhora significativa ao varia a taxa de subamostragem, portanto optou-se por mante-la igual a 1.
In [62]:
score_list = []
error_list = []
subsamples = np.linspace(0.1,1,10)
for subsample in subsamples:
model = GradientBoostingClassifier(subsample=subsample)
scores = cross_val_score(model, X, y, cv=5, n_jobs=-1, scoring='roc_auc')
score_list.append(np.mean(scores))
error_list.append(np.std(scores))
plt.errorbar(subsamples, score_list, yerr=error_list)
Out[62]:
Utilizar um número de features menor que o total costuma reduzir a variância do modelo em troca de um aumento do bias. Ao variar essa taxa não observou-se um valor estatisticamente superior ao outro, devido as grandes barras de erro. Optou-se por utilizar a taxa que apresentou maior AUC médio, no caso, 0.2.
In [63]:
score_list = []
error_list = []
max_features_list = np.linspace(0.1,1,10)
for max_features in max_features_list:
model = GradientBoostingClassifier(max_features=max_features)
scores = cross_val_score(model, X, y, cv=5, n_jobs=-1, scoring='roc_auc')
score_list.append(np.mean(scores))
error_list.append(np.std(scores))
plt.errorbar(max_features_list, score_list, yerr=error_list)
Out[63]:
In [68]:
score_list = []
error_list = []
max_depths = [1,2,3,4,5,6]
for max_depth in max_depths:
model = GradientBoostingClassifier(max_features=0.2, max_depth=max_depth)
scores = cross_val_score(model, X, y, cv=5, n_jobs=-1, scoring='roc_auc')
score_list.append(np.mean(scores))
error_list.append(np.std(scores))
plt.errorbar(max_depths, score_list, yerr=error_list)
Out[68]:
In [71]:
score_list = []
error_list = []
min_samples_splits = np.logspace(-5,0,10)
for min_samples_split in min_samples_splits:
model = GradientBoostingClassifier(max_features=0.2, min_samples_split=min_samples_split)
scores = cross_val_score(model, X, y, cv=5, n_jobs=-1, scoring='roc_auc')
score_list.append(np.mean(scores))
error_list.append(np.std(scores))
plt.errorbar(min_samples_splits, score_list, yerr=error_list)
plt.xscale('log')
In [72]:
score_list = []
error_list = []
min_samples_splits = [2,3,4,5,6,7,8,9,10]
for min_samples_split in min_samples_splits:
model = GradientBoostingClassifier(max_features=0.2, min_samples_split=min_samples_split)
scores = cross_val_score(model, X, y, cv=5, n_jobs=-1, scoring='roc_auc')
score_list.append(np.mean(scores))
error_list.append(np.std(scores))
plt.errorbar(min_samples_splits, score_list, yerr=error_list)
Out[72]:
In [83]:
score_list = []
error_list = []
min_samples_leafs = [1,10,20,30,40,50,60,70]
for min_samples_leaf in min_samples_leafs:
model = GradientBoostingClassifier(max_features=0.2, min_samples_leaf=min_samples_leaf)
scores = cross_val_score(model, X, y, cv=5, n_jobs=-1, scoring='roc_auc')
score_list.append(np.mean(scores))
error_list.append(np.std(scores))
plt.errorbar(min_samples_leafs, score_list, yerr=error_list)
Out[83]:
A busca de parâmetros até agora considerou os parâmetros como sendo independentes. Buscar a melhor combinação de todos os parâmetros testados seria computacionalmente inviável. Optou-se por comparar o melhor modelo obtido com a busca de parametros supondo independencia com um modelo obtido através de uma busca por combinações aleatóreas de parâmentros.
In [108]:
param_dist = {"max_depth": randint(1, 11),
"max_features": randint(1, 11),
"min_samples_split": randint(2, 11),
"min_samples_leaf": randint(1, 11),
"subsample": uniform(0.1, 0.9),
"learning_rate": uniform(0.001, 0.9),
"loss" : ['deviance', 'exponential']}
# run randomized search
model = GradientBoostingClassifier()
random_search = RandomizedSearchCV(model, param_distributions=param_dist,n_iter=50)
random_search.fit(X, y)
Out[108]:
In [133]:
model = GradientBoostingClassifier(max_features=0.2)
scores = cross_val_score(model, X, y, cv=5, n_jobs=-1, scoring='roc_auc')
print 'Busca Independente: %f +- %f' % (np.mean(scores), np.std(scores))
model = random_search.best_estimator_
scores = cross_val_score(model, X, y, cv=5, n_jobs=-1, scoring='roc_auc')
print 'Busca Aleatórea: %f +- %f' % (np.mean(scores), np.std(scores))
In [101]:
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.4)
model = random_search.best_estimator_
y_pred = model.fit(X_train, y_train).predict_proba(X_test)
y_pred = np.array([x[1] for x in y_pred])
fpr, tpr, thresholds = roc_curve(y_test+1, y_pred, pos_label=2)
roc_auc = auc(fpr, tpr)
plt.figure()
lw = 2
plt.plot(fpr, tpr, color='darkorange',
lw=lw, label='ROC curve (area = %0.2f)' % roc_auc)
plt.plot([0, 1], [0, 1], color='navy', lw=lw, linestyle='--')
plt.xlim([0.0, 1.0])
plt.ylim([0.0, 1.05])
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('Receiver operating characteristic example')
plt.legend(loc="lower right")
plt.show()
In [ ]:
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.4)
y_pred = model.fit(X_train,y_train).predict(X_test)
labels = ['Pago', 'Default']
mtx = confusion_matrix(y_test, y_pred)
mtx = [x/float(sum(x)) for x in mtx]
sns.heatmap(pd.DataFrame(mtx, columns=labels, index=labels), annot=True, fmt=".2f", linewidths=.5)
plt.xlabel('Valor Classificado')
plt.ylabel('Valor Real')
In [62]:
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.4)
with tf.Graph().as_default():
net = tflearn.input_data(shape=[None, X.shape[1]])
net = tflearn.fully_connected(net, 100, activation='relu')
net = tflearn.dropout(net, 0.5)
net = tflearn.fully_connected(net, 100, activation='relu')
net = tflearn.dropout(net, 0.5)
net = tflearn.fully_connected(net, 2, activation='softmax')
net = tflearn.regression(net, optimizer='adam', loss='categorical_crossentropy')
Y = np.array([[float(a), float(1-a)] for a in y_train])
lm = tflearn.DNN(net)
lm.fit(X_train.as_matrix(), Y, validation_set=0.3, show_metric=True, batch_size=200, n_epoch=50, snapshot_epoch=False)
In [61]:
y_pred = lm.predict(X_test)
y_pred = [a[0] for a in y_pred]
fpr, tpr, thresholds = roc_curve(y_test+1, y_pred, pos_label=2)
roc_auc = auc(fpr, tpr)
plt.figure()
lw = 2
plt.plot(fpr, tpr, color='darkorange',
lw=lw, label='ROC curve (area = %0.2f)' % roc_auc)
plt.plot([0, 1], [0, 1], color='navy', lw=lw, linestyle='--')
plt.xlim([0.0, 1.0])
plt.ylim([0.0, 1.05])
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title ('Receiver operating characteristic example')
plt.legend(loc="lower right")
plt.show()
Observando as curvas ROC dos modelos, plotadas juntas, observamos que a curva do Gradient Boosting Trees, se manteve assima das demais o tempo todos, ou seja, para qualquer valor que escolhermos como limiar de operação ele irá ter um desempenho superior aos demais.
Os demais modelos apresentaram um desempenho relativamente semelhante, de forma que o segundo melhor modelo dependo do ponto de operação que se utilizar.
In [2]:
train = pd.read_csv('../data/processed/train.csv', index_col=0)
test = pd.read_csv('../data/processed/test.csv', index_col=0)
del train.index.name
del test.index.name
X_train = train.drop('default.payment.next.month', axis=1)
X_train_red = PCA(n_components=10).fit_transform(X_train)
y_train = train['default.payment.next.month']
X_test = test.drop('default.payment.next.month', axis=1)
X_test_red = PCA(n_components=10).fit_transform(X_test)
y_test = test['default.payment.next.month']
def get_crossval(model, scoring='accuracy', knn=False):
result_list = cross_val_score(model, X_train, y_train, cv=5, n_jobs=-1, scoring=scoring)
if knn:
result_list = cross_val_score(model, X_train_red, y_train, cv=5, n_jobs=-1, scoring=scoring)
return '%.3f +- %.3f' % (np.mean(result_list), np.std(result_list))
def get_auc(model, knn=False):
y_pred = [y[1] for y in model.fit(X_train, y_train).predict_proba(X_test)]
if knn:
y_pred = [y[1] for y in model.fit(X_train_red, y_train).predict_proba(X_test_red)]
fpr, tpr, _ = roc_curve(y_test+1, y_pred, pos_label=2)
return auc(fpr, tpr)
def get_accuracy(model, knn=False):
y_pred = model.fit(X_train, y_train).predict(X_test)
if knn:
y_pred = model.fit(X_train_red, y_train).predict(X_test_red)
return accuracy_score(y_test, y_pred)
def plot_roc(model, label, knn=False):
y_pred = [y[1] for y in model.fit(X_train, y_train).predict_proba(X_test)]
if knn:
y_pred = [y[1] for y in model.fit(X_train_red, y_train).predict_proba(X_test_red)]
fpr, tpr, _ = roc_curve(y_test+1, y_pred, pos_label=2)
plt.plot(fpr, tpr, label=label)
plt.plot([0, 1], [0, 1], color='navy',linestyle='--')
plt.xlim([0.0, 1.0])
plt.ylim([0.0, 1.05])
plt.xlabel('Taxa de Falsos Positivos')
plt.ylabel('Taxa de Verdadeiros Positivos')
plt.title(u'Comparação dos Modelos')
plt.legend(loc="lower right")
def plot_conf_mtx(model, title, pos, knn=False):
y_pred = model.fit(X_train,y_train).predict(X_test)
if knn:
y_pred = model.fit(X_train_red,y_train).predict(X_test_red)
columns = ['Pago', 'Default']
index = ['Pago', 'Default']
if pos[0]==0:
columns = ['', '']
if pos[1]==1:
index = ['', '']
mtx = confusion_matrix(y_test, y_pred)
mtx = [x/float(sum(x)) for x in mtx]
sns.heatmap(pd.DataFrame(mtx, columns=columns, index=index), annot=True, fmt=".2f", linewidths=.5)
if pos[0]==1:
plt.xlabel('Valor Classificado')
if pos[1]==0:
plt.ylabel('Valor Real')
plt.title(title)
In [3]:
nb = GaussianNB()
lr = LogisticRegression(C=10)
knn = KNeighborsClassifier(n_neighbors=62)
gbt = GradientBoostingClassifier(max_features=0.2)
plt.figure()
plot_roc(nb, 'Naive Bayes')
plot_roc(lr, 'Logistic Regression')
plot_roc(knn, 'KNN', True)
plot_roc(gbt, 'Gradient Boosting Trees')
plt.show()
In [4]:
roc = pd.DataFrame()
roc.set_value('Naive Bayes', '5-fold CV', get_crossval(nb, 'roc_auc'))
roc.set_value('KNN', '5-fold CV', get_crossval(knn, 'roc_auc', True))
roc.set_value('Logistic Regression', '5-fold CV', get_crossval(lr, 'roc_auc'))
roc.set_value('Gradient Boosting Trees', '5-fold CV', get_crossval(gbt, 'roc_auc'))
roc.set_value('Naive Bayes', 'Conjunto de Teste', get_auc(nb))
roc.set_value('KNN', 'Conjunto de Teste', get_auc(knn, True))
roc.set_value('Logistic Regression', 'Conjunto de Teste', get_auc(lr))
roc.set_value('Gradient Boosting Trees', 'Conjunto de Teste', get_auc(gbt))
roc
Out[4]:
Por último são comparadas as eficiências utilizando como ponto de operação 0.5, ou seja, caso o valor retornado polo modelo esteja acima de 0.5 a amostra é classificada como 1 (default), caso seja abaixo é classificada como 0 (pago).
Podemos observar que nesta configuração Naive Bayes tem uma queda drástica no seu resultado, enquanto KNN e Logistic Regression ficaram mais próximos do resultado do Gradient Boosting Trees, embora este ainda tenha sido superior.
In [5]:
accuracy = pd.DataFrame()
accuracy.set_value('Naive Bayes', '5-fold CV', get_crossval(nb))
accuracy.set_value('Logistic Regression', '5-fold CV', get_crossval(lr))
accuracy.set_value('KNN', '5-fold CV', get_crossval(knn, 'accuracy', True))
accuracy.set_value('Gradient Boosting Trees', '5-fold CV', get_crossval(gbt))
accuracy.set_value('Naive Bayes', 'Conjunto de Teste', get_accuracy(nb))
accuracy.set_value('Logistic Regression', 'Conjunto de Teste', get_accuracy(lr))
accuracy.set_value('KNN', 'Conjunto de Teste', get_accuracy(knn, True))
accuracy.set_value('Gradient Boosting Trees', 'Conjunto de Teste', get_accuracy(gbt))
accuracy
Out[5]:
Ainda tratando da operação em 0.5, podemos observar que os modelos, com excessão do Naive Bayes, possuem uma taxa de acerto muito alta de detecção de pagadores, porém apresentam um resultado ruim na detecção de defaults. O Naive Bayes, apresentou a taxa mais alta dentre eles na detecção de defaults, porém teve um resultado bastante abaixo na detecção de pagadores.
Como foi visto anteriormente, Gradient Boosting Trees apresenta uma curva ROC superior do que a dos demais modelos. Portanto bastaria variar o limiar de corte para um valor que permita maior equilíbrio entre detecção e falso alarme. Para isso utilizou-se F1-score.
In [6]:
plt.subplot(221)
plot_conf_mtx(nb, 'Naive Bayes', pos=(0,0))
plt.subplot(222)
plot_conf_mtx(knn, 'KNN', pos=(0,1), knn=True)
plt.subplot(223)
plot_conf_mtx(lr, 'Logistic Regression', pos=(1,0))
plt.subplot(224)
plot_conf_mtx(gbt, 'Gradient Boosting Trees', pos=(1,1))
F1-score é uma métrica para classificadores binários que considera tanto a taxa de detecção quanto a de falso alarme no seu cálculo, fazendo uma especie de média harmónica de ambas. Ao calcular o f1-score para diferentes limiares de corte pode-se observar que ele possui um ponto máximo. Ulilizando este ponto como limiar de corte o modelo apresentou a matiz de confusão descrita abaixo, muito mais equilibranda que a primeira.
In [7]:
def get_f1(fpr, tpr):
return 2 * tpr * (1-fpr) / (tpr + 1-fpr)
y_pred = [y[1] for y in gbt.fit(X_train, y_train).predict_proba(X_test)]
fpr_l, tpr_l, tsh_l = roc_curve(y_test+1, y_pred, pos_label=2)
f1_list = []
for fpr, tpr in zip(fpr_l, tpr_l):
f1_list.append(get_f1(fpr, tpr))
print 'Máximo F1-score:', max(f1_list)
print 'Limiar de corte:', (max(f1_list))
plt.plot(tsh_l, f1_list)
Out[7]:
In [8]:
thd = 0.2
def get_value(y):
if y > thd:
return 1
else:
return 0
y_pred = [get_value(y) for y in y_pred]
labels = ['Pago', 'Default']
mtx = confusion_matrix(y_test, y_pred)
mtx = [x/float(sum(x)) for x in mtx]
sns.heatmap(pd.DataFrame(mtx, columns=labels, index=labels), annot=True, fmt=".2f", linewidths=.5)
plt.xlabel('Valor Classificado')
plt.ylabel('Valor Real')
plt.title('Gradient Boosting Trees - Maior F1-score')
Out[8]:
Dos modelos testados o Gradient Boosting Trees apresentou o melhor resultado para qualquer ponto de operação. Os dados foram pré-processados de forma a:
A melhor configuração do modelo para estes dados, utilizou os seguintes parâmetros sem nenhuma redução de dimensionalidade: