Práctica de Aprendizaje Supervisado

Se ha abordado el ejemplo Give Me Some Credit. La descripción del problema, así como los datos empleados se encuentran disponibles aquí.

Lo primero ha sido el estudio de las features dsiponibles, en concreto nos preocupaba la cantidad de datos disponibles y de buena calidad (ausencia de NaN) así como las relaciones entre las variables, para descubrir estructuras iniciales.

Se ha probado la combinación y creación de diferentes features, pero no se adjuntan en esta memoria, pues se ha comprobado que no mejoraban los resultados respecto al uso del conjunto inicial.

Veamos las features de las que disponemos:

Análisis de features

Se adjunta un diagrama que muestra la correlación y distribución entre las distintas features:



In [77]:

    
corrmat = trainB[np.r_[INS,['MonthlyIncome']]].corr()
axs = scatter_matrix(trainB[np.r_[INS,['MonthlyIncome']]],diagonal='kde',figsize=(16,6))
for ax in axs[:,0]: ax.set_ylabel(ax.get_ylabel(),rotation=0)
for ax in axs[-1,:]: ax.set_xlabel(ax.get_xlabel(),rotation=10)
names = ['var {}'.format(i) for i in xrange(len(corrmat.columns))]
corrmat.index, corrmat.columns = names, names



In [79]:

    
corrmat

Se observa una alta correlación entre la estructura de Dues, por otro lado se desea conocer también que cantidad de datos ausentes existen en los conjuntos de test y de training:



In [70]:

    
f, axs = plt.subplots(1,2,figsize=(16,3))
plot_NaNs_info(1,axs[0])
plot_NaNs_info(2,axs[1])

Y si están distribuidos de la misma forma en función de las clases en la variable de salida:



In [83]:

    
f, axs = plt.subplots(1,1,figsize=(8,3))
plot_NaNs_info(3,axs)

Se observa que los conjuntos tienen la misma proporción de NaNs y en las mismas clases (MonthlyIncome y #Dependents), Sin embargo los conjuntos no son iguales en absoluto, para empezar porque son de diferente tamaño y luego presentan diferencias en los datos (veamos por ejemplo la variable MonthlyIncome):



In [66]:

    
a = test.MonthlyIncome.copy()
b = traininig.MonthlyIncome.copy()
a.sort(), b.sort()
(a-b).plot(title="Monthly Income test-train differences")









    Out[66]:





<matplotlib.axes._subplots.AxesSubplot at 0x7fc0746783d0>

Debido a estas características se ha decidido eliminar aquellos registros en los que el #Dependents sea NaN, Creando un conjunto de datos de entrenamiento que será el utilizado para entrenar a los modelos, Veamos cómo cargamos los datos:



In [84]:

    
import os,re
import urllib2
import numpy as np
import pandas as pd
from bs4 import  BeautifulSoup as bs4
from variables.variables import VarRepl
import matplotlib.pyplot as plt

ROOT = '/home/x/Documentos/'
if 'OS' in list(os.environ) and any(re.findall('^[wW]in.+',os.environ['OS'])):
#    ROOT = 'C:/Users/Xavi/Documents/u-tad/Modulo12/'
    ROOT = 'D:/Users/jserrano/Documents/u-tad/Modulo12/'
path = 'https://xavi783.github.io/data/GiveMeSomeCredit/'

# Cargamos los datos y la agrupación de datos por categoría.

categories = ['dues','loans','familiar','out']
memory = bs4(urllib2.urlopen('https://github.com/xavi783/data/blob/master/GiveMeSomeCredit/memoria.md'))
#Creamos el diccionario de variables:
dv = {c:np.r_[map(lambda y: y.text.strip(),x.select('li'))] for c,x in zip(categories,memory.select("article[itemprop] .task-list")[2:])}
dv['familiarB'] = np.r_[list(set(dv['familiar']) - {u'MonthlyIncome'})]
replDict = [('NumberOf','#'),
            ('Number','#'),
            ('DaysPastDueNotWorse','DPastDue'),
            ('LinesAndLoans','Lines'),
            ('LoansOrLines','Loans'),
            ('RevolvingUtilizationOf','Util')]
var_groups = VarRepl(dv,replDict)
traininig = pd.read_csv(path+'cs-training.csv',index_col=0)
test = pd.read_csv(path+'cs-test.csv',index_col=0)

todrop = [traininig[np.isnan(traininig['NumberOfDependents'])].index, \
          traininig[np.isnan(traininig['SeriousDlqin2yrs'])].index]
trainB = traininig.copy()
for d in todrop:
    if any(d):
        trainB = trainB.drop(d)

Muestreo y constrrucción del estimador

Para entrenar a los modelos, reduciendo el overfiting, se ha decidido usar cross validation. Sin embargo, desde mi experiencia personal, la mayoría de los modelos de machine learning no funcionan bien si las clases a predecir se encuentran muy sesgadas, aunque ciertos métodos de machine learning son más sensible a este fenómeno que otros.

Por este motivo se ha decido realizar un muestreo basado en cross validation pero un poco más complejo de lo habitual. Veamos el procedimiento con una imagen que intenta explicarlo mejor:

Se divide el total del conjunto de datos en 5 muestreos (representados por cajas verdes)
Cada muestreo se divide a su vez en un conjunto de training (cajas rojas) y uno de validación (cajas azules) Estos conjuntos tienen miembros de la clase de salida en la misma proporción que el conjunto total, es decir, cada caja roja/azul tiene aproximadamente un 7% de observaciones de clientes que cometieron impago y un 93% de clientes que no lo hicieron.
Como hemos dicho el tener conjuntos de entrenamiento tan sesgados no suele funcionar bien (de hecho se ha realizado el entrenamiento sobre estos conjuntos y se ha comprobado que era pésimo). Sin embargo, los conjuntos de validación son buenos, ya que son buenos representantes de la población total.
Para balancear los conjuntos de entrenamiento se construyen subgrupos (cajas amarillas dentro de las cajas rojas) de tal forma que la clase de salida esté equilibrada (50% de clientes cometen impago, frente a los que no), para eso se dividen los cliente que no comenten imapgo en grupos del mismo tamaño que el de clientes que lo cometen, y se crean grupos, repitiendo el conjunto de clientes que comenten impago por cada uno de los nuevos grupos de cliente sin impagar. </div>
</div>
Una vez se han creado los grupos, entrenamos un modelo diferente con cada uno de ellos (un modelo por cada caja amarilla). Pero para construir nuestro sistema final, debemos combinar las señales de estos modelos.

Existen varias formas de hacerlo, en general recurrimos a lo que se denominan métodos enssemble (métodos, que de forma similar a una red de neuronas, combinan cada modelo para obtener una señal final). De hecho la formulación puede realizarse de forma completamente análoga a la de una red neuronal.

Para simplificar este trabajo, hemos creado un método enssemble por cada caja verde, hemos medido el área bajo la curva roc en cada conjunto de validación y nos quedamos con el método que tiene mejor resultados, aunque se podrían combinar todos los modelos en un sólo enssemble.

Además en laimágen se e que hay una función de treeshold (muy similar a las funciones de activación de una red neuronal) y en este caso hemos escogido que sea de tipo step.

Este umbral se utiliza de la siguiente manera:

$$ ES = \left\{ \begin{array}{rl} 1 &\mbox{ $\sum{w_i*S_i} \geq treeshold$ } \\ 0 &\mbox{ $\sum{w_i*S_i} < treeshold$ } \end{array} \right. $$

Donde:

$ES$: señal del ensamblador (=1: el sliente comete impago, =0: el cliente no comete impago)
$w_i$: peso de cada modelo en el ensamblador
$S_i$: señal de cada modelo
$treeshold$: umbral de decisión

Veamos el código:

Muestreo estratificado:



In [85]:

    
import pickle as pkl
from sklearn.cross_validation import StratifiedKFold, KFold
from sklearn import metrics
KF = StratifiedKFold(np.r_[trainB[var_groups.out].values.flat], n_folds=5) # cada fold tiene un 10% del total

## balanceamos los test:
Fs = [f for f,k in KF]
Ks = [k for f,k in KF]
folds = {}
for i,k in enumerate(Ks):
    folds[i] = {'k':[],'f':Fs[i]}
    x = trainB[var_groups.out].iloc[k,:]
    xl = np.r_[(x==1).values.flat]
    subfolding_ratio = np.floor(1.*(x==0).sum()/(x==1).sum())[0]
    KF2 = KFold((x==0).sum()[0], subfolding_ratio)
    Ks2, Fs2 = [k2 for f2,k2 in KF2], [f2 for f2,k2 in KF2]
    for k2 in Ks2:
        folds[i]['k'] += [np.r_[k[k2],k[xl]]]

A continuación creamos las clases y funciones necesarias para crear el modelo ensamblador:



In [86]:

    
def build_classififers(Classifier, params, name):
    if ~os.path.exists(ROOT+name+'.pk'):   
        clf,P,mP,tP = {},{},{},{}
        INS = np.r_[var_groups.dues,var_groups.familiarB,var_groups.loans]
        OUT = var_groups.out
        for i,kf in folds.iteritems():        
            clf[i] = {x:Classifier(**params) for x in xrange(len(kf['k']))}
            for j,k in enumerate(kf['k']):
                clf[i][j].fit(trainB[INS].iloc[k].values, np.r_[trainB[OUT].iloc[k].values.flat])
            P[i] = np.c_[[forest.predict(trainB[INS].iloc[kf['f']]) for forest in clf[i].itervalues()]]   
            tP[i] = np.r_[trainB[OUT].iloc[kf['f']].values.flat]    
        roc = []
        # arboles equiponderados, cuando la imputación supera el treeshold, entonces se asigna clase =1, sino clase =0
        # máximo treeshold como óptimo    
        for treeshold in np.r_[0.:1.1:.1]:
            for i,kf in folds.iteritems():     
                mP[i] = P[i].sum(0)/(1.*P[i].shape[0])
                mP[i][mP[i]>treeshold],mP[i][mP[i]<=treeshold]=1,0
            roc += [np.mean([metrics.roc_auc_score(tp,p) for tp,p in zip(tP.itervalues(),mP.itervalues())])]    
        optimum = np.r_[roc].argmax()
        treeshold = np.r_[0.:1.1:.1][optimum]
        # Seleccionamos el conjunto de modelos óptimo entre todos los folds (con mejor roc) dado un treeshold óptimo
        for i,kf in folds.iteritems():     
            mP[i] = P[i].sum(0)/(1.*P[i].shape[0])
            mP[i][mP[i]>treeshold],mP[i][mP[i]<=treeshold]=1,0
        roc = np.r_[[metrics.roc_auc_score(tp,p) for tp,p in zip(tP.itervalues(),mP.itervalues())]]
        optimum = roc.argmax()
        pkl.dump(dict(clf=clf,P=P,tP=tP,mP=mP,treeshold=treeshold,optimum=optimum),open(ROOT+name+'.pk','w+'))    
    else:
        print 'Ya existe el archivo: '+ROOT+name+'.pk'

class MultiClassifier(object):
    def __init__(self,dictClassifiers):
        self.__dict__.update(dictClassifiers)
        self.classifiers = dictClassifiers
        self.weights = np.ones(len(self.classifiers))/len(self.classifiers)
        
    def predict_proba(self,test,weigths=None):
        if type(weigths)!=type(None):
            self.weights = weigths
        p = [self.weights[i]*clf.predict_proba(test)[:,1] for i,clf in self.classifiers.iteritems()]
        p = np.r_[p].sum(0)
        return np.c_[np.arange(len(p)),p]
    
# Los mejores son: Random Forest y SVM:
def predict_proba(classifier, test):
    probs = [[index + 1, x[1]] for index, x in enumerate(classifier.predict_proba(test[INS]))]
    np.savetxt(ROOT+'submission.csv', probs, delimiter=',', fmt='%d,%f', header='Id,Probability', comments = '')
    return probs

Los métodos anteriores se utilizarán de la siguiente forma:

build_classififers: esta función nos permitirá entrenar todos los modelos (uno por cada caja amarilla) de un clúster del cross validation (caja verde)
MultiClassifier: esta clase actúa como interfaz del método predict_proba, lo que hace es coger un conjunto de clasificadores y clacula las posibilidades de obtener la clase de salida del conjunto, utilizando el ensamblador como hemos descrito.
predict_proba: nos permite calcular las probabilidades de la variable target y crear el archivo para la submission a kaggle.

Veamos cómo lo hemos aplicado para construir el archivo de submission:



In [ ]:

    
# Configuramos los clasificadores a probar:
from sklearn.svm import SVC
from sklearn.naive_bayes import GaussianNB
from sklearn.ensemble import RandomForestClassifier
from sklearn.neighbors import KNeighborsClassifier

Classifiers = [RandomForestClassifier, KNeighborsClassifier, SVC, GaussianNB]
params = [dict(n_estimators=100, n_jobs=3), {}, {}, {}]
names = ['randomforest','knn','svm','naivebayes']

# Creamos los calsificadores
map(build_classififers,Classifiers,params,names)

# Calculamos AUC en función del cuál vamos a escoger el mejormétodo
# (en este caso Random Forest, porque tiene un AUC de 0.78, es el mayor)
auc = map(lambda name: (lambda d: metrics.roc_auc_score(d['tP'][d['optimum']],d['mP'][d['optimum']]))\
          (pkl.load(open(ROOT+name+'.pk','r'))),names) 
test.loc[np.isnan(test['NumberOfDependents']),'NumberOfDependents'] = \
    np.mean(test.loc[~np.isnan(test['NumberOfDependents']),'NumberOfDependents'])

d = pkl.load(open(ROOT+names[0]+'.pk','r'))
predict_proba(MultiClassifier(d['clf'][d['optimum']]),test)

clasif2 = {k:v for k,v in enumerate(list(d['clf'][d['optimum']].itervalues())+list(d1['clf'][d1['optimum']].itervalues()))}
predict_proba(MultiClassifier(clasif2),test)

Tras dar estos pasos y enviar el archivo de submission, se ha obtenido un AUC de la curva roc de probabilidades de 0.84, quedando en un puesto 594.

Regresión lineal de la variable Monthly Income

Ya se han estudiado las relaciones entre las variables de entrada y se ha observado que no había relaciones lineales entre ellas, por lo que sin crear nuevas features, es complicado obtener una buena regresión lineal que se ajuste bien a la variable Monthly Income, aún así dejamos una muestra de cómo se llevaría a cabo, utilizando por ejemplo, lo que hemos denominado "conjunto de variables relacionadas con la estructura familiar".

Veamos primero una descomposición en componentes principales, por si fuéramos capaces de inferir alguna relación entre las variables.



In [113]:

    
from pandas import DataFrame as df_
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler as ZScore

zscore = ZScore()
pca = PCA(n_components=2)
#INS2 = np.r_[INS,['MonthlyIncome']]
INS2 = np.r_[var_groups.familiar,var_groups.loans]
names_INS2 = np.r_[var_groups.repl.familiar,var_groups.repl.loans]

x = trainB[INS2][~np.isnan(trainB.MonthlyIncome)]
pca.fit(zscore.fit_transform(x))
plot_components()









    Out[113]:





u'age, MonthlyIncome, #Dependents, DebtRatio, UtilUnsecuredLines, #OpenCreditLines, #RealEstateLoans'

Como puede verse, si eliminamos la estructura de impagos, las componentes necesarias para representar de forma significativa las variables, se reducen de 3 a 2. Por tanto a la hora de ralizar la regresión lineal no tendremos en cuanta esas variables.

A continuación mostramos como haríamos una regresión lineal con esas variables:



In [143]:

    
# Regresión lineal:
from sklearn import metrics
from sklearn.linear_model import LinearRegression
INS = np.r_[var_groups.familiarB,var_groups.loans]    
lm = LinearRegression()
lma = lm.fit(trainB[INS][~np.isnan(trainB.MonthlyIncome)],trainB['MonthlyIncome'][~np.isnan(trainB.MonthlyIncome)])
y_true = trainB['MonthlyIncome'][~np.isnan(trainB.MonthlyIncome)]
y_pred = lma.predict(trainB[INS][~np.isnan(trainB.MonthlyIncome)])
accuracy = {'r2':metrics.r2_score(y_true,y_pred),\
            'MSE':metrics.mean_squared_error(y_true,y_pred)}
accuracy









    Out[143]:





{'MSE': 202373043.24858341, 'r2': 0.021960910265643174}

Cómo se puede ver, la estimación no es demasiado precisa, se podría mejorar esta estimación escogiendo mejor las variables explicativas o incluso aplicando el mismo método de cross validation que para el método ensamblador antes descrito.

	var 0	var 1	var 2	var 3	var 4	var 5	var 6	var 7	var 8	var 9
var 0	1.000000	0.985573	0.981770	-0.057981	-0.001235	-0.006064	-0.002680	-0.052329	-0.029059	-0.010217
var 1	0.985573	1.000000	0.991976	-0.052066	-0.000987	-0.007118	-0.010922	-0.068600	-0.038429	-0.011116
var 2	0.981770	0.991976	1.000000	-0.056152	-0.000995	-0.007972	-0.010176	-0.078079	-0.044271	-0.012743
var 3	-0.057981	-0.052066	-0.056152	1.000000	-0.005808	0.025773	-0.213303	0.157961	0.041663	0.037717
var 4	-0.001235	-0.000987	-0.000995	-0.005808	1.000000	0.004059	0.001557	-0.011021	0.006633	0.007124
var 5	-0.006064	-0.007118	-0.007972	0.025773	0.004059	1.000000	-0.040673	0.051044	0.116461	-0.028712
var 6	-0.002680	-0.010922	-0.010176	-0.213303	0.001557	-0.040673	1.000000	0.065322	0.124684	0.062647
var 7	-0.052329	-0.068600	-0.078079	0.157961	-0.011021	0.051044	0.065322	1.000000	0.430833	0.091455
var 8	-0.029059	-0.038429	-0.044271	0.041663	0.006633	0.116461	0.124684	0.430833	1.000000	0.124959
var 9	-0.010217	-0.011116	-0.012743	0.037717	0.007124	-0.028712	0.062647	0.091455	0.124959	1.000000