PlaidML: An open source portable deep learning engine


We’re pleased to announce the next step towards deep learning for every device and platform. Today Vertex.AI is releasing PlaidML, our open source portable deep learning engine. Our mission is to make deep learning accessible to every person on every device, and we’re building PlaidML to help make that a reality. We’re starting by supporting the most popular hardware and software already in the hands of developers, researchers, and students. The initial version of PlaidML runs on most existing PC hardware with OpenCL-capable GPUs from NVIDIA, AMD, or Intel. Additionally, we’re including support for running the widely popular Keras framework on top of Plaid to allow existing code and tutorials to run unchanged.

PlaidML: An open source portable deep learning engine

Baidu are bringing HPC Techniques to Deep Learning

Via Baidu Research Blog.

The Ring all-reduce approach came to save a lot of work when training deep neural networks. The approach to propagate and update the gradients (and control the convergence of the model) are well explained below:

The Ring Allreduce

The main issue with the simplistic communication strategy described above was that the communication cost grew linearly with the number of GPUs in the system. In contrast, a ring allreduce is an algorithm for which the communication cost is constant and independent of the number of GPUs in the system, and is determined solely by the slowest connection between GPUs in the system; in fact, if you only consider bandwidth as a factor in your communication cost (and ignore latency), the ring allreduce is an optimal communication algorithm [4]. (This is a good estimate for communication cost when your model is large, and you need to send large amounts of data a small number of times.)

The GPUs in a ring allreduce are arranged in a logical ring. Each GPU should have a left neighbor and a right neighbor; it will only ever send data to its right neighbor, and receive data from its left neighbor.

The algorithm proceeds in two steps: first, a scatter-reduce, and then, an allgather. In the scatter-reduce step, the GPUs will exchange data such that every GPU ends up with a chunk of the final result. In the allgather step, the GPUs will exchange those chunks such that all GPUs end up with the complete final result.

More can be found here. To implement there’s a Github project called Tensor All-reduce that can be used for distributed deep learning.

Baidu are bringing HPC Techniques to Deep Learning

Horovod: Uber’s Open Source Distributed Deep Learning Framework for TensorFlow

Via Uber Engineering Blog.

Yet, this is another tool for Deep Learning, but I think that those guys hit the nail exposing and fixing one of the major concerns about Tensor Flow that is distributed training.

When Uber needed to use Deep Learning they found some endeavors to use the conventional Data Parallelism architecture.  Using Data Parallelism arch they can distribute the training using several instances in parallel and when the gradients for every batch are calculated in each instance (node/worker) these gradients are propagated for all nodes and averaged to control the convergence (update) of the model in the training phase. The following image explains better than words.

But using this architecture Uber faced two problems that were a) the right ratio of worker to parameter servers (to avoid/deal with network and processing bottleneck) and b) the complexity of TensorFlow code (more details here).

To avoid these problems they used an idea of a 2009 paper  “Bandwidth Optimal All-reduce Algorithms for Clusters of Workstations” called Ring all-reduce.

They explain the workflow of this approach:

In the ring-allreduce algorithm, each of N nodes communicates with two of its peers 2*(N-1) times. During this communication, a node sends and receives chunks of the data buffer. In the first N-1 iterations, received values are added to the values in the node’s buffer. In the second N-1 iterations, received values replace the values held in the node’s buffer. Baidu’s paper suggests that this algorithm is bandwidth-optimal, meaning that if the buffer is large enough, it will optimally utilize the available network.

The implementations details can be found here.



Horovod: Uber’s Open Source Distributed Deep Learning Framework for TensorFlow

Tensorflow sucks (?)

This post of Nico’s blog has good points about why Pytorch even without all Google support and money is taking out users from Tensor Flow.

[…]The most interesting question to me is why Google chose a purely declarative paradigm for Tensorflow in spite of the obvious downsides of this approach. Did they feel that encapsulating all the computation in a single computation graph would simplify executing models on their TPU’s so they can cut Nvidia out of the millions of dollars to be made from cloud hosting of deep learning powered applications? It’s difficult to say. Overall, Tensorflow does not feel like a pure open source project for the common good. Which I would have no problem with, had their design been sound. In comparison with beautiful Google open source projects out there such as Protobuf, Golang, and Kubernetes, Tensorflow falls dramatically short.

While declarative paradigms are great for UI programming, there are many reasons why it is a problematic choice for deep learning.

Take the React Javascript library as an example, the standard choice today for interactive web applications. In React, the complexity of how data flows through the application makes sense to be hidden from the developer, since Javascript execution is generally orders of magnitudes faster than updates to the DOM. React developers don’t want to worry about the mechanics of how state is propagated, so long as the end user experience is “good enough”.

On the other hand, in deep learning, a single layer can literally execute billions of FLOP’s! And deep learning researchers care very much about the mechanics of how computation is done and want fine control because they are constantly pushing the edge of what’s possible (e.g. dynamic networks) and want easy access to intermediate results.[…]

Tensorflow sucks (?)

How to pass environment variables in Jupyter Notebook

(Sharing some personal suffering)

One thing that gets me mad it’s to make several .csv/.txt files in my computer to perform some analysis. I personally prefer to connect directly in some RDBMS (Redshift) and get the data in some straightforward way and store the query inside the Jupiter Notebook.

The main problem with this approach is: a high number of people put their passwords inside the notebooks/scripts and this is very unsafe. (You don’t need to believe me, check it by yourself)

I was trying to pass the environment variables in a traditional way using export VARIABLE_NAME=xptoSomeValue  but after starting the Jupyter Notebook I get the following error:



KeyError                                  Traceback (most recent call last)
<ipython-input-13-2288aa3f6b7a> in <module>()
      2 import os
----> 4 HOST = os.environ['REDSHIFT_HOST']
      5 PORT = os.environ['REDSHIFT_PORT']
      6 USER = os.environ['REDSHIFT_USER']

/usr/local/Cellar/python/2.7.13/Frameworks/Python.framework/Versions/2.7/lib/python2.7/UserDict.pyc in __getitem__(self, key)
     38         if hasattr(self.__class__, "__missing__"):
     39             return self.__class__.__missing__(self, key)
---> 40         raise KeyError(key)
     41     def __setitem__(self, key, item):[key] = item
     42     def __delitem__(self, key): del[key]


For some reason, this approach didn’t work. I make a small workaround to start using some environmental variables when I call of jupyter notebook command in that way:

env REDSHIFT_HOST='myRedshiftHost' REDSHIFT_USER='flavio.clesio' REDSHIFT_PORT='5439' REDSHIFT_DATA='myDatabase' REDSHIFT_PASS='myVeryHardPass' jupyter notebook

I hope it helps!

How to pass environment variables in Jupyter Notebook

Six reasons your boss must send you to Spark Summit Europe 2017

It’s redundant to say that Apache Spark is becoming the most prominent open-source big data cluster-computing framework in the last 2 years, where this technology not only shattered old paradigms of general purpose distributed data processing, but also built a very vibrant, innovation-driven, and receptive community.

This is my first time at Spark Summit, and for me personally, it’s a great time as Machine Learning professional to be part of such event that has grown dramatically in the last 2 years only.

Here in Brazil we do not have such tradition to invest in conferences (that are some cultural reasons involved that needed to break down in another blog post), but this is the six reasons that your boss must send you to Spark Summit Europe 2017:

  1. Accomplish more than the rest: While some your company competitors are heavily busy making re-work in old frameworks, your company can stay focused to solve real problems that permit scalability for your business using bleeding edge technologies.
  2. Stay ahead of the game: You can choose one of these two sentences to put in your resumé: 1) “Worked with Apache Spark, the most prominent open-source cluster-computing framework for Big Data Projects“; or 2) “Worked with <<Put some obsolete framework the needs a couple USD millions to be deployed and have 70% fewer features than Apache Spark and the most stable version was written 9 years ago and the whole marketing are migrating>>”. It’s up to you.
  3. Connect with Apache Spark experts: In Spark Summit you’ll meet some real dealers of Apache Spark, not someone with marketing pitch (no offense) offering difficulties (e.g. closed-source, buggy platform) to sell facilities (e.g. never-ending-consulting-until-drain-your-entire-budget style, sell (buggy) plugins, add-ons, etc… ). Some of Spark experts are Tim Hunter, Tathagata Das,  Sue Ann Hong, Holden Karau, to name a few.
  4. Network that matters: I mean people with shared interest in enthusiasm over an open-source framework Apache Spark and technology, headhunters of good companies that understand that data plays a strong role at business; not some B.S. artist or pseudo-tech-cloaked-sellers someone else.
  5. Applied knowledge produce innovation, and innovation produce results: Some cases using Apache Spark to innovate and help business – Saving more than US$ 3 million using Apache Spark and Machine Learning, managing 300TB data workload using Apache Spark, real-time anomaly detection in some systems, changing the game of digital marketing using Apache Spark,  and predicting traffic using weather data.
  6. Opting out will destroy your business and your career: Refuse to get knowledge and apply that it’s the fast way to destroy your career with stagnation in old methods/process/platforms and become obsolete in a few months. For your company, opting out of innovation or learning new methods and technologies that can help to scale the business or enhance productivity, it’s a good way to get out of business in a few years.

To register and learn more about the event, please visit Spark Summit 2017 and follow spark_summit on Twitter.

Six reasons your boss must send you to Spark Summit Europe 2017

See you at Spark Summit Europe 2017

In October 26, my friend Eiti Kimura and I will provide a talk called Preventing leakage and monitoring distributed systems with Machine Learning at Spark Summit Europe 2017 where we’ll show our solution to monitoring a highly complex distributed system using Apache Spark as a tool for Machine Learning.

We’re very excited to share our experience in this journey, and how we solved a complex problem using a simple solution that saved more than US$ 3 million in the last 19 months.

See you at Spark Summit at Dublin.

See you at Spark Summit Europe 2017

Como criar um Virtualenv no Python sem bullshit

Via Eiti Kimura.

Direto e reto:

1) Realize a instalação do virtualenv pelo pip

$ pip install virtualenv

2) Faça a definição do seu diretório

$ mkdir deep-learning-virtual-env

3) Após a definição, entre no diretório

$ cd deep-learning-virtual-env

4) Faça a inicialização do seu virtualenv

$ virtualenv .

5) Com isso realize a ativação do seu virtualenv

$ source bin/activate

6) Para facilitar o seu trabalho, criamos até mesmo um arquivo de requirements com o Theano, Keras, Jupyter Notebook, Scikit-Learn. Para fazer isso basta rodar o seguinte comando:

$ pip install -r requirements.key




Como criar um Virtualenv no Python sem bullshit

RLScore: Regularized Least-Squares Learners

Uma boa alternativa para ensemble quando a dimensionalidade dos datasets for alta, ou as alternativas com Elastic Net, Lasso e Ridge não derem a convergência desejada.

RLScore: Regularized Least-Squares Learners

RLScore is a Python open source module for kernel based machine learning. The library provides implementations of several regularized least-squares (RLS) type of learners. RLS methods for regression and classification, ranking, greedy feature selection, multi-task and zero-shot learning, and unsupervised classification are included. Matrix algebra based computational short-cuts are used to ensure efficiency of both training and cross-validation. A simple API and extensive tutorials allow for easy use of RLScore.

Regularized least squares (RLS) is a family of methods for solving the least-squares problem while using regularization to further constrain the resulting solution.

RLS is used for two main reasons. The first comes up when the number of variables in the linear system exceeds the number of observations. In such settings, the ordinary least-squares problem is ill-posed and is therefore impossible to fit because the associated optimization problem has infinitely many solutions. RLS allows the introduction of further constraints that uniquely determine the solution.

The second reason that RLS is used occurs when the number of variables does not exceed the number of observations, but the learned model suffers from poor generalization. RLS can be used in such cases to improve the generalizability of the model by constraining it at training time. This constraint can either force the solution to be “sparse” in some way or to reflect other prior knowledge about the problem such as information about correlations between features. A Bayesian understanding of this can be reached by showing that RLS methods are often equivalent to priors on the solution to the least-squares problem.

To sse in Depth

1) $ pip install rlscore
2) $ export CFLAGS="-I /usr/local/lib/python2.7/site-packages/numpy/core/include $CFLAGS"

Original post

In [1]:
# Import libraries
import numpy as np
from rlscore.learner import RLS
from rlscore.measure import sqerror
from rlscore.learner import LeaveOneOutRLS
In [2]:
# Function to load dataset and split in train and test sets
def load_housing():
    D = np.loadtxt("/Volumes/PANZER/Github/learning-space/Datasets/02 - Classification/housing_data.txt")
    X = D[:,:-1] # Independent variables
    Y = D[:,-1]  # Dependent variable
    X_train = X[:250]
    Y_train = Y[:250]
    X_test = X[250:]
    Y_test = Y[250:]
    return X_train, Y_train, X_test, Y_test
In [3]:
def print_stats():
    X_train, Y_train, X_test, Y_test = load_housing()
    print("Housing data set characteristics")
    print("Training set: %d instances, %d features" %X_train.shape)
    print("Test set: %d instances, %d features" %X_test.shape)

if __name__ == "__main__":
Housing data set characteristics
Training set: 250 instances, 13 features
Test set: 256 instances, 13 features

Linear regression with default parameters

In [4]:
# Function to train RLS method
def train_rls():
    #Trains RLS with default parameters (regparam=1.0, kernel='LinearKernel')
    X_train, Y_train, X_test, Y_test = load_housing()
    learner = RLS(X_train, Y_train)
    #Leave-one-out cross-validation predictions, this is fast due to
    #computational short-cut
    P_loo = learner.leave_one_out()
    #Test set predictions
    P_test = learner.predict(X_test)
    # Stats
    print("leave-one-out error %f" %sqerror(Y_train, P_loo))
    print("test error %f" %sqerror(Y_test, P_test))
    #Sanity check, can we do better than predicting mean of training labels?
    print("mean predictor %f" %sqerror(Y_test, np.ones(Y_test.shape)*np.mean(Y_train)))

if __name__=="__main__":
leave-one-out error 25.959399
test error 25.497222
mean predictor 81.458770

Choosing regularization parameter with leave-one-out

Regularization parameter with grid search in exponential grid to catch the lowest LOO-CV error.

In [5]:
def train_rls():
    #Select regparam with leave-one-out cross-validation
    X_train, Y_train, X_test, Y_test = load_housing()
    learner = RLS(X_train, Y_train)
    best_regparam = None
    best_error = float("inf")
    #exponential grid of possible regparam values
    log_regparams = range(-15, 16)
    for log_regparam in log_regparams:
        regparam = 2.**log_regparam
        #RLS is re-trained with the new regparam, this
        #is very fast due to computational short-cut
        #Leave-one-out cross-validation predictions, this is fast due to
        #computational short-cut
        P_loo = learner.leave_one_out()
        e = sqerror(Y_train, P_loo)
        print("regparam 2**%d, loo-error %f" %(log_regparam, e))
        if e < best_error:
            best_error = e
            best_regparam = regparam
    P_test = learner.predict(X_test)
    print("best regparam %f with loo-error %f" %(best_regparam, best_error)) 
    print("test error %f" %sqerror(Y_test, P_test))

if __name__=="__main__":
regparam 2**-15, loo-error 24.745479
regparam 2**-14, loo-error 24.745463
regparam 2**-13, loo-error 24.745431
regparam 2**-12, loo-error 24.745369
regparam 2**-11, loo-error 24.745246
regparam 2**-10, loo-error 24.745010
regparam 2**-9, loo-error 24.744576
regparam 2**-8, loo-error 24.743856
regparam 2**-7, loo-error 24.742982
regparam 2**-6, loo-error 24.743309
regparam 2**-5, loo-error 24.750966
regparam 2**-4, loo-error 24.786243
regparam 2**-3, loo-error 24.896991
regparam 2**-2, loo-error 25.146493
regparam 2**-1, loo-error 25.537315
regparam 2**0, loo-error 25.959399
regparam 2**1, loo-error 26.285436
regparam 2**2, loo-error 26.479254
regparam 2**3, loo-error 26.603001
regparam 2**4, loo-error 26.801196
regparam 2**5, loo-error 27.352322
regparam 2**6, loo-error 28.837002
regparam 2**7, loo-error 32.113350
regparam 2**8, loo-error 37.480625
regparam 2**9, loo-error 43.843555
regparam 2**10, loo-error 49.748687
regparam 2**11, loo-error 54.912297
regparam 2**12, loo-error 59.936226
regparam 2**13, loo-error 65.137825
regparam 2**14, loo-error 70.126118
regparam 2**15, loo-error 74.336978
best regparam 0.007812 with loo-error 24.742982
test error 24.509981

Training with RLS and simultaneously selecting the regularization parameter with leave-one-out using LeaveOneOutRLS

In [6]:
def train_rls():
    #Trains RLS with automatically selected regularization parameter
    X_train, Y_train, X_test, Y_test = load_housing()
    # Grid search
    regparams = [2.**i for i in range(-15, 16)]
    learner = LeaveOneOutRLS(X_train, Y_train, regparams = regparams)
    loo_errors = learner.cv_performances
    P_test = learner.predict(X_test)
    print("leave-one-out errors " +str(loo_errors))
    print("chosen regparam %f" %learner.regparam)
    print("test error %f" %sqerror(Y_test, P_test))

if __name__=="__main__":
leave-one-out errors [ 24.74547881  24.74546295  24.74543138  24.74536884  24.74524616
  24.74501033  24.7445764   24.74385625  24.74298177  24.74330936
  24.75096639  24.78624255  24.89699067  25.14649266  25.53731465
  25.95939943  26.28543584  26.47925431  26.6030015   26.80119588
  27.35232186  28.83700156  32.11334986  37.48062503  43.84355496
  49.7486873   54.91229746  59.93622566  65.1378248   70.12611801
chosen regparam 0.007812
test error 24.509981

Learning nonlinear predictors using kernels

RLS using a non-linear kernel function.

In [7]:
def train_rls():
    #Selects both the gamma parameter for Gaussian kernel, and regparam with loocv
    X_train, Y_train, X_test, Y_test = load_housing()
    regparams = [2.**i for i in range(-15, 16)]
    gammas = regparams
    best_regparam = None
    best_gamma = None
    best_error = float("inf")
    for gamma in gammas:
        #New RLS is initialized for each kernel parameter
        learner = RLS(X_train, Y_train, kernel="GaussianKernel", gamma=gamma)
        for regparam in regparams:
            #RLS is re-trained with the new regparam, this
            #is very fast due to computational short-cut
            #Leave-one-out cross-validation predictions, this is fast due to
            #computational short-cut
            P_loo = learner.leave_one_out()
            e = sqerror(Y_train, P_loo)
            #print "regparam", regparam, "gamma", gamma, "loo-error", e
            if e < best_error:
                best_error = e
                best_regparam = regparam
                best_gamma = gamma
    learner = RLS(X_train, Y_train, regparam = best_regparam, kernel="GaussianKernel", gamma=best_gamma)
    P_test = learner.predict(X_test)
    print("best parameters gamma %f regparam %f" %(best_gamma, best_regparam))
    print("best leave-one-out error %f" %best_error)
    print("test error %f" %sqerror(Y_test, P_test))
if __name__=="__main__":
best parameters gamma 0.000031 regparam 0.000244
best leave-one-out error 21.910837
test error 16.340877

Binary classification and Area under ROC curve

In [8]:
from rlscore.utilities.reader import read_svmlight

# Load dataset and stats
def print_stats():
    X_train, Y_train, foo = read_svmlight("/Volumes/PANZER/Github/learning-space/Datasets/02 - Classification/a1a.t")
    X_test, Y_test, foo = read_svmlight("/Volumes/PANZER/Github/learning-space/Datasets/02 - Classification/a1a")
    print("Adult data set characteristics")
    print("Training set: %d instances, %d features" %X_train.shape)
    print("Test set: %d instances, %d features" %X_test.shape)

if __name__=="__main__":
Adult data set characteristics
Training set: 30956 instances, 123 features
Test set: 1605 instances, 119 features
In [ ]:
from rlscore.learner import RLS
from rlscore.measure import accuracy
from rlscore.utilities.reader import read_svmlight

def train_rls():
    # Train ans test datasets    
    X_train, Y_train, foo = read_svmlight("/Volumes/PANZER/Github/learning-space/Datasets/02 - Classification/a1a.t")
    X_test, Y_test, foo = read_svmlight("/Volumes/PANZER/Github/learning-space/Datasets/02 - Classification/a1a", X_train.shape[1])
    learner = RLS(X_train, Y_train)
    best_regparam = None
    best_accuracy = 0.
    #exponential grid of possible regparam values
    log_regparams = range(-15, 16)
    for log_regparam in log_regparams:
        regparam = 2.**log_regparam
        #RLS is re-trained with the new regparam, this
        #is very fast due to computational short-cut
        #Leave-one-out cross-validation predictions, this is fast due to
        #computational short-cut
        P_loo = learner.leave_one_out()
        acc = accuracy(Y_train, P_loo)
        print("regparam 2**%d, loo-accuracy %f" %(log_regparam, acc))
        if acc > best_accuracy:
            best_accuracy = acc
            best_regparam = regparam
    P_test = learner.predict(X_test)
    print("best regparam %f with loo-accuracy %f" %(best_regparam, best_accuracy)) 
    print("test set accuracy %f" %accuracy(Y_test, P_test))

if __name__=="__main__":
RLScore: Regularized Least-Squares Learners

Accelerating the XGBoost algorithm using GPU computing

A fronteira final em relação ao uso com GPU de um dos mais poderosos algoritmos de todos os tempos está aqui.

Abstract: We present a CUDA based implementation of a decision tree construction algorithm within the gradient boosting library XGBoost. The tree construction algorithm is executed entirely on the GPU and shows high performance with a variety of datasets and settings, including sparse input matrices. Individual boosting iterations are parallelized, combining two approaches. An interleaved approach is used for shallow trees, switching to a more conventional radix sort based approach for larger depths. We show speedups of between 3-6x using a Titan X compared to a 4 core i7 CPU, and 1.2x using a Titan X compared to 2x Xeon CPUs (24 cores). We show that it is possible to process the Higgs dataset (10 million instances, 28 features) entirely within GPU memory. The algorithm is made available as a plug-in within the XGBoost library and fully supports all XGBoost features including classification, regression and ranking tasks. 

Accelerating the XGBoost algorithm using GPU computing

Análise de Múltipla Correspondência no R para o problema de Churn

Via Data Science Plus

Analytical challenges in multivariate data analysis and predictive modeling include identifying redundant and irrelevant variables. A recommended analytics approach is to first address the redundancy; which can be achieved by identifying groups of variables that are as correlated as possible among themselves and as uncorrelated as possible with other variable groups in the same data set. On the other hand, relevancy is about potential predictor variables and involves understanding the relationship between the target variable and input variables.
Multiple correspondence analysis (MCA) is a multivariate data analysis and data mining tool for finding and constructing a low-dimensional visual representation of variable associations among groups of categorical variables. Variable clustering as a tool for identifying redundancy is often applied to get a first impression of variable associations and multivariate data structure.
The motivations of this post are to illustrate the applications of: 1) preparing input variables for analysis and predictive modeling, 2) MCA as a multivariate exploratory data analysis and categorical data mining tool for business insights of customer churn data, and 3) variable clustering of categorical variables for the identification of redundant variables.

Análise de Múltipla Correspondência no R para o problema de Churn

Deep Learning AMI Amazon Web Services

Para quem quer escalar processamento em Machine Learning e não tem grana para comprar GPUs, o Deep Learning AMI da Amazon é uma ótima alternativa em termos de custos.

The Deep Learning AMI is an Amazon Linux image supported and maintained by Amazon Web Services for use on Amazon Elastic Compute Cloud (Amazon EC2). It is designed to provide a stable, secure, and high performance execution environment for deep learning applications running on Amazon EC2. It includes popular deep learning frameworks, including MXNet, Caffe, Tensorflow, Theano, CNTK and Torch as well as packages that enable easy integration with AWS, including launch configuration tools and many popular AWS libraries and tools. It also includes the Anaconda Data Science Platform for Python2 and Python3. Amazon Web Services provides ongoing security and maintenance updates to all instances running the Amazon Linux AMI. The Deep Learning AMI is provided at no additional charge to Amazon EC2 users.

The AMI Ids for the Deep Learning Amazon Linux AMI are the following:
us-east-1 : ami-e7c96af1
us-west-2: ami-dfb13ebf
eu-west-1: ami-6e5d6808

Release tags/Branches used for the DW Frameworks:
MXNet : v0.9.3 tag
Tensorflow : v1.0.0 tag
Theano : rel-0.8.2 tag
Caffe : rc5 tag
CNTK : v2.0beta12.0 tag
Torch : master branch
Keras : 1.2.2 tag

Deep Learning AMI Amazon Web Services

Ferramenta para Machine Learning – MLJAR

Para quem busca uma alternativa paga para Machine Learning em ambientes fora da própria infraestrutura o MLJAR pode ser a resposta.


MLJAR is a human-first platform for machine learning.
It provides a service for prototyping, development and deploying pattern recognition algorithms.
It makes algorithm search and tuning painless!


You pay for computational time used for models training, predictions and data analysis. 1 credit is 1 computation hour on machine with 8 CPU and 15GB RAM. Computational time is aggregated per second basis.

Ferramenta para Machine Learning – MLJAR

Akid: Uma biblioteca de Redes Neurais para pesquisa e produção

Finalmente começaram a pensar em eliminar esse vale entre ciência/academia e indústria.

Akid: A Library for Neural Network Research and Production from a Dataism Approach – Shuai Li
Abstract: Neural networks are a revolutionary but immature technique that is fast evolving and heavily relies on data. To benefit from the newest development and newly available data, we want the gap between research and production as small as possibly. On the other hand, differing from traditional machine learning models, neural network is not just yet another statistic model, but a model for the natural processing engine — the brain. In this work, we describe a neural network library named {\texttt akid}. It provides higher level of abstraction for entities (abstracted as blocks) in nature upon the abstraction done on signals (abstracted as tensors) by Tensorflow, characterizing the dataism observation that all entities in nature processes input and emit out in some ways. It includes a full stack of software that provides abstraction to let researchers focus on research instead of implementation, while at the same time the developed program can also be put into production seamlessly in a distributed environment, and be production ready. At the top application stack, it provides out-of-box tools for neural network applications. Lower down, akid provides a programming paradigm that lets user easily build customized models. The distributed computing stack handles the concurrency and communication, thus letting models be trained or deployed to a single GPU, multiple GPUs, or a distributed environment without affecting how a model is specified in the programming paradigm stack. Lastly, the distributed deployment stack handles how the distributed computing is deployed, thus decoupling the research prototype environment with the actual production environment, and is able to dynamically allocate computing resources, so development (Devs) and operations (Ops) could be separated. 

Akid: Uma biblioteca de Redes Neurais para pesquisa e produção

Churn-at-Risk: Aplicação de Survival Analysis no controle de churn de assinaturas em Telecom


Um dos assuntos mais recorrentes em qualquer tipo de serviço de assinatura é como reduzir o Churn (saída de clientes), dado que conquistar novos clientes é bem mais difícil (e caro) do que manter os antigos.

Cerca de 70% das empresas sabem que é mais barato manter um cliente do que ter que ir atrás de um novo.

Fazendo uma analogia simples, o lucro dos serviços de assinatura são como uma espécie de sangue na corrente sanguínea de uma empresa e uma interrupção de qualquer natureza prejudica todo o negócio, dado que esse é um modelo de receita que se baseia na recorrência de tarifação e não no desenvolvimento, ou mesmo venda de outros produtos.

Em modelos de negócios baseados no volume de pessoas que estão dispostas a terem uma cobrança recorrente o negócio fica bem mais complicado, dado que diferentemente de produtos que tem uma elasticidade maior o fluxo de receita é extremamente sujeito aos sabores do mercado e dos clientes.

Dentro desse cenário, para todas as empresas que tem o seu fluxo de receita baseado nesse tipo de business, saber quando um cliente entrará em uma situação de saída através do cancelamento do serviço (Churn) é fundamental para criar mecanismos de retenção mais efetivos, ou mesmo criação de réguas de contato com os clientes para evitar ou minimizar a chance de um cliente sair da base de dados.

Sendo assim, qualquer mecanismo ou mesmo esforço para minimizar esse efeito é de grande valia. Nos baseamos na teoria estatística buscar respostas para as seguintes perguntas:

  • Como diminuir o Churn?
  • Como identificar um potencial cliente que irá entrar em uma situação de Churn? Quais estratégias seguir para minimizar esse Churn?
  • Quais réguas de comunicação com os clientes devemos ter para entender os motivos que estão fazendo um assinante cancelar o serviço e quais são as estratégias de customer winback possíveis nesse cenário?

E pra responder essa pergunta, fomos buscar as respostas na análise de sobrevivência dado que essa área da estatística é uma das que lidam melhor em termos de probabilidade de tempo de vida com dados censurados, seja de materiais (e.g. tempo de falha de algum sistema mecânico) ou no tempo de vida de pessoas propriamente ditas (e.g. dado uma determinada posologia qual é a estimativa de um paciente sobreviver a um câncer), e no nosso caso quanto tempo de vida um assinante tem até deixar cancelar a sua assinatura.

Análise de Sobrevivência

A análise de sobreviência é uma técnica estatístisca que foi desenvolvida na medicina e tem como principal finalidade estimar o tempo de sobrevivência ou tempo de morte de um determinado paciente dentro de um horizonte do tempo.

O estimador de Kaplan-Meier (1958) utiliza uma função de sobrevivência que leva em consideração uma divisão entre o número de observações que não falharam no tempo t pelo número total de observações no estudo em que cada intervalo de tempo tem-se o número de falhas/mortes/churn distintos bem como é calculado o risco de acordo com o número de indivíduos restantes no tempo subsequente.

Já o estimador Nelson-Aalen (1978) é um estimador que tem as mesmas características do Kaplan-Meier, com a diferença que esse estimador trabalha com uma função de sobrevivência que é a cumulative hazard rate function.

Os elementos fundamentais para caracterização de um estudo que envolve análise de sobrevivência são, o (a) tempo inicial, (b)escala de medida do intervalo de tempo e (c) se o evento de churn ocorreu.

Os principais artigos são de Aalen (1978), Kaplan-Meier (1958) e Cox (1972).

Esse post não tem como principal objetivo dar algum tipo de introdução sobre survival analysis, dado que tem muitas referências na internet sobre o assunto e não há nada a ser acrescentado nesse sentido por este pobre blogueiro.

Assim como a análise de cohort, a análise de sobrevivência tem como principal característica ser um estudo de natureza longitudinal, isto é, os seus resultados tem uma característica de temporalidade seja em aspectos de retrospecção, quanto em termos de perspectivas, isso é, tem uma resposta tipicamente temporal para um determinado evento de interesse.

O que vamos usar como forma de comparação amostral é o comportamento longitudinal, de acordo com determinadas características de amostragens diferentes ao longo do tempo, e os fatores que influenciam no churn.

Devido a questões óbvias de NDA não vamos postar aqui características que possam indicar qualquer estratégia de negócios ou mesmo caracterização de alguma informação de qualquer natureza.

Podemos dizer que a análise de sobrevivência aplicada em um caso de telecom, pode ajudar ter uma estimativa em forma de probabilidade em relação ao tempo em que uma assinatura vai durar até o evento de churn (cancelamento) e dessa forma elaborar estratégias para evitar esse evento, dado que adquirir um novo cliente é mais caro do que manter um novo e entra totalmente dentro de uma estratégia de Customer Winback (Nota: Esse livro Customer Winback do Jill Griffin e do Michael Lowenstein é obrigatório para todos que trabalham com serviços de assinaturas ou negócios que dependam de uma recorrência muito grande como comércio).

No nosso caso o tempo de falha ou tempo de morte, como estamos falando de serviços de assinaturas, o nosso evento de interesse seria o churn, ou cancelamento da assinatura. Em outras palavras teríamos algo do tipo Time-to-Churn ou um Churn-at-Risk. Guardem esse termo.


Usamos dados de dois produtos antigos em que os dados foram anonimizados e aplicados um hash de embaralhamento uniforme (que obedece uma distribuição específica) nos atributos (por questões de privacidade) que são:

  • id = Identificador do registro;
  • product = produto;
  • channel = canal no qual o cliente entrou na base de dados;
  • free_user = flag que indica se o cliente entrou na base em gratuidade ou não;
  • user_plan = se o usuário é pré-pago ou pós-pago;
  • t = tempo que o assinante está na base de dados; e
  • c = informa se o evento de interesse (no caso o churn (cancelamento da assinatura) ocorreu ou não.

Eliminamos o efeito de censura à esquerda retirando os casos de reativações, dado que queríamos entender a jornada do assinante como um todo sem nenhum tipo de viés relativo a questões de customer winback. Em relação à censura à direita temos alguns casos bem específicos que já se passaram alguns meses desde que essa base de dados foi extraída.

Um aspecto técnico importante a ser considerado é que esses dois produtos estão em categorias de comparabilidade, dado que sem isso nenhum tipo de caractericação seria nula.

No fim dessa implementação teremos uma tabela de vida em relação a esses produtos.


Primeiramente vamos importar as bibliotecas: Pandas (para manipulação de dados), matplotlib (para a geração de gráficos), e lifelines para aplicação da análise de sobrevivência:

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import lifelines

Após realizar a importação das bibliotecas, vamos ajustar o tamanho das imagens para uma melhor visualização:

%pylab inline
pylab.rcParams['figure.figsize'] = (14, 9)

Vamos realizar o upload da nossa base de dados criando um objecto chamado df e usando a classe read_csv do Pandas:

df = pd.read_csv('')

Vamos checar a nossa base de dados:

id product channel free_user user_plan t c
0 3315 B HH 1 0 22 0
1 2372 A FF 1 1 16 0
2 1098 B HH 1 1 22 0
3 2758 B HH 1 1 4 1
4 2377 A FF 1 1 29 0

Então como podemos ver temos as 7 variáveis na nossa base de dados.

Na sequência vamos importar a biblioteca do Lifelines, em especial o estimador de KaplanMaier:

from lifelines import KaplanMeierFitter

kmf = KaplanMeierFitter()

Após realizar a importação da classe relativa ao estimador de Kaplan Meier no objeto kmf, vamos atribuir as nossas variáveis de tempo (T) e evento de interesse (C)

T = df["t"]

C = df["c"]

O que foi feito anteriormente é que buscamos no dataframe df o array t e atribuímos no objeto T, e buscamos o array da coluna c no dataframe e atribuímos no objeto C.

Agora vamos chamar o método fit usando esses dois objetos no snippet abaixo:, event_observed=C )
<lifelines.KaplanMeierFitter: fitted with 10000 observations, 6000 censored>
Objeto ajustado, vamos agora ver o gráfico relativo a esse objeto usando o estimador de Kaplan Meier.
plt.title('Survival function of Service Valued Add Products');
plt.ylabel('Probability of Living (%)')
plt.xlabel('Lifespan of the subscription (in days)')
<matplotlib.text.Text at 0x101b24a90>

Como podemos ver no gráfico, temos algumas observações pertinentes, quando tratamos a probabilidade de sobrevivência desses dois produtos no agregado que são:

  • Logo no primeiro dia há uma redução substancial do tempo de sobrevivência da assinatura em aproximadamente 22%;
  • Há um decaimento quase que linear depois do quinto dia de assinatura; e
  • Depois do dia número 30, a probabilidade de sobrevivência de uma assinatura é de aproximadamente de 50%. Em outras palavras: depois de 30 dias, metade dos novos assinantes já estarão fora da base de assinantes.

No entanto, vamos plotar a mesma função de sobrevivência considerando os intervalos de confiança estatística.

plt.title('Survival function of Service Valued Add Products - Confidence Interval in 85-95%');
plt.ylabel('Probability of Living (%)')
plt.xlabel('Lifespan of the subscription')
<matplotlib.text.Text at 0x10ad8e0f0>

Contudo nesse modelo inicial temos duas limitações claras que são:

  • Os dados no agregado não dizem muito em relação à dinâmicas que podem estar na especificidade de alguns atributos/dimensões;
  • Não são exploradas as dimensões (ou quebras) de acordo com os atributos que vieram na base de dados; e
  • Não há a divisão por produto.

Para isso, vamos começar a entrar no detalhe em relação a cada uma das dimensões e ver o que cada uma tem de influência em relação à função de sobrevivência.

Vamos começar realizando a quebra pela dimensão que determina se o cliente entrou via gratuidade ou não (free_user).

ax = plt.subplot(111)

free = (df["free_user"] == 1)[free], event_observed=C[free], label="Free Users")
kmf.plot(ax=ax, ci_force_lines=True)[~free], event_observed=C[~free], label="Non-Free Users")
kmf.plot(ax=ax, ci_force_lines=True)
plt.title("Lifespans of different subscription types");
plt.ylabel('Probability of Living (%)')
<matplotlib.text.Text at 0x10ad8e908>

Este gráfico apresenta algumas informações importantes para os primeiros insights em relação a cada uma das curvas de sobrevivência em relação ao tipo de gratuidade oferecida como fator de influência para o churn que são:

  • Os assinantes que entram como não gratuitos (i.e. não tem nenhum tipo de gratuidade inicial) após o 15o dia apresenta um decaimento brutal de mais de 40% da chance de sobrevivência (tratando-se do intervalo de confiança);
  • Após o 15o dia os assinantes que não desfrutam de gratuidade tem a sua curva de sobrevivência em uma relativa estabilidade em torno de 60% na probabilidade de sobrevivência até o período censurado;
  • Ainda nos usuários sem gratuidade, dado o grau de variabilidade do intervalo de confiança podemos tirar como conclusão que muitos cancelamentos estão ocorrendo de forma muito acelerada, o que deve ser investigado com mais calma pelo time de produtos; e
  • Já os usuários que entram via gratuidade (i.e. ganham alguns dias grátis antes de serem tarifados) apresenta um nível de decaimento do nível de sobrevivência maior seja no período inicial, quando ao longo do tempo, contudo uma estabilidade é encontrada ao longo de toda a série sem maiores sobressaltos.

Dado essa análise inicial das curvas de sobrevivência, vamos avaliar agora as probabilidades de sobrevivência de acordo com o produto.

ax = plt.subplot(111)

product = (df["product"] == "A")[product], event_observed=C[product], label="Product A")
kmf.plot(ax=ax, ci_force_lines=True)[~product], event_observed=C[~product], label="Product B")
kmf.plot(ax=ax, ci_force_lines=True)

plt.title("Survival Curves of different Products");
plt.ylabel('Probability of Living (%)')
<matplotlib.text.Text at 0x10aeaabe0>

Este gráfico apresenta a primeira distinção entre os dois produtos de uma forma mais clara.

Mesmo com os intervalos de confiança com uma variação de 5%, podemos ver que o produto A (linha azul) tem uma maior probabilidade de sobrevivência com uma diferença percentual de mais de 15%; diferença essa amplificada depois do vigésimo dia.

Em outras palavras: Dado um determinada safra de usuários, caso o usuário entre no produto A o mesmo tem uma probabilidade de retenção de cerca de 15% em relação a um usuário que por ventura entre no produto B, ou o produto A apresenta uma cauda de retenção superior ao produto B.

Empiricamente é sabido que um dos principais fatores de influência de produtos SVA são os canais de mídia os quais esses produtos são oferecidos.

O canal de mídia é o termômetro em que podemos saber se estamos oferencendo os nossos produtos para o público alvo correto.

No entanto para um melhor entendimento, vamos analisar os canais nos quais as assinaturas são originadas.

A priori vamos normalizar a variável channel para realizar a segmentação dos canais de acordo com o conjunto de dados.

df['channel'] = df['channel'].astype('category');
channels = df['channel'].unique()

Após normalização e transformação da variável para o tipo categórico, vamos ver como está o array.

[HH, FF, CC, AA, GG, ..., BB, EE, DD, JJ, ZZ]
Length: 11
Categories (11, object): [HH, FF, CC, AA, ..., EE, DD, JJ, ZZ]

Aqui temos a representação de 11 canais de mídia os quais os clientes entraram no serviço.

Com esses canais, vamos identificar a probabilidade de sobrevivência de acordo com o canal.

for i,channel_type in enumerate(channels):
    ax = plt.subplot(3,4,i+1)
    ix = df['channel'] == channel_type T[ix], C[ix], label=channel_type )
    kmf.plot(ax=ax, legend=True)
    if i==0:
        plt.ylabel('Probability of Survival by Channel (%)')
Fazendo uma análise sobre cada um desses gráficos temos algumas considerações sobre cada um dos canais:
  • HH, DD: Uma alta taxa de mortalidade (churn) logo antes dos primeiros 5 dias, o que indica uma característica de efemeridade ou atratividade no produto para o público desse canal de mídia.
  • FF: Apresenta menos de 10% de taxa de mortalidade nos primeiros 20 dias, e tem um padrão muito particular depois do 25o dia em que praticamente não tem uma mortalidade tão alta. Contém um intervalo de confiança com uma oscilação muito forte.
  • CC: Junto com o HH apesar de ter uma taxa de mortalidade alta antes do 10o dia, apresenta um grau de previsibilidade muito bom, o que pode ser utilizado em estratégias de incentivos de mídia que tenham que ter uma segurança maior em termos de retenção a médio prazo.
  • GG, BB: Apresentam uma boa taxa de sobrevivência no inicio do período, contudo possuem oscilações severas em seus respectivos intervalos de confiança. Essa variável deve ser considerada no momento de elaboração de uma estratégia de investimento nesses canais.
  • JJ: Se houvesse uma definição de incerteza em termos de sobrevivência, esse canal seria o seu melhor representante. Com os seus intervalos de confiança oscilando em mais de 40% em relação ao limite inferior e superior, esse canal de mídia mostra-se extremamente arriscado para os investimentos, dado que não há nenhum tipo de regularidade/previsibilidade de acordo com esses dados.
  • II: Apesar de ter um bom grau de previsibilidade em relação à taxa de sobrevivência nos primeiros 10 dias, após esse período tem uma curva de hazard muito severa, o que indica que esse tipo de canal pode ser usado em uma estratégia de curto prazo.
  • AA, EE, ZZ: Por haver alguma forma de censura nos dados, necessitam de mais análise nesse primeiro momento. (Entrar no detalhe dos dados e ver se é censura à direita ou algum tipo de truncamento).

Agora que já sabemos um pouco da dinâmica de cada canal, vamos criar uma tabela de vida para esses dados.

A tabela de vida nada mais é do que uma representação da função de sobrevivência de forma tabular em relação aos dias de sobrevivência.

Para isso vamos usar a biblioteca utils do lifelines para chegarmos nesse valor.

from lifelines.utils import survival_table_from_events

Biblioteca importada, vamos usar agora as nossas variáveis T e C novamente para realizar o ajuste da tabela de vida.

lifetable = survival_table_from_events(T, C)

Tabela importada, vamos dar uma olhada no conjunto de dados.

print (lifetable)
          removed  observed  censored  entrance  at_risk
0            2250      2247         3     10000    10000
1             676       531       145         0     7750
2             482       337       145         0     7074
3             185       129        56         0     6592
4             232        94       138         0     6407
5             299        85       214         0     6175
6             191        73       118         0     5876
7             127        76        51         0     5685
8             211        75       136         0     5558
9            2924        21      2903         0     5347
10            121        27        94         0     2423
11             46        27        19         0     2302
12             78        26        52         0     2256
13            111        16        95         0     2178
14             55        35        20         0     2067
15            107        29        78         0     2012
16            286        30       256         0     1905
17            156        23       133         0     1619
18            108        18        90         0     1463
19             49        11        38         0     1355
20             50        17        33         0     1306
21             61        13        48         0     1256
22            236        23       213         0     1195
23             99         6        93         0      959
24            168         9       159         0      860
25            171         7       164         0      692
26             58         6        52         0      521
27             77         2        75         0      463
28             29         6        23         0      386
29            105         1       104         0      357
30             69         0        69         0      252
31            183         0       183         0      183

Diferentemente do R que possuí a tabela de vida com a porcentagem relativa à probabilidade de sobrevivência, nesse caso vamos ter que fazer um pequeno ajuste para obter a porcentagem de acordo com o atributo entrance e at_risk.

O ajuste se dará da seguinte forma:

survivaltable = lifetable.at_risk/np.amax(lifetable.entrance)

Ajustes efetuados, vamos ver como está a nossa tabela de vida.

0     1.0000
1     0.7750
2     0.7074
3     0.6592
4     0.6407
5     0.6175
6     0.5876
7     0.5685
8     0.5558
9     0.5347
10    0.2423
11    0.2302
12    0.2256
13    0.2178
14    0.2067
15    0.2012
16    0.1905
17    0.1619
18    0.1463
19    0.1355
20    0.1306
21    0.1256
22    0.1195
23    0.0959
24    0.0860
25    0.0692
26    0.0521
27    0.0463
28    0.0386
29    0.0357
30    0.0252
31    0.0183
Name: at_risk, dtype: float64

Vamos transformar a nossa tabela de vida em um objeto do pandas para melhor manipulação do conjunto de dados.

survtable = pd.DataFrame(survivaltable)

Para casos de atualização de Churn-at-Risk podemos definir uma função que já terá a tabela de vida e poderá fazer a atribuição da probabilidade de sobrevivência de acordo com os dias de sobrevivência.

Para isso vamos fazer uma função simples usando o próprio python.

def survival_probability( int ):
   print ("The probability of Survival after", int, "days is", survtable["at_risk"].iloc[int]*100, "%") 

Nesse caso vamos ver a chance de sobrevivência usando o nosso modelo Kaplan-Meier já ajustado para uma assinatura que tenha 22 dias de vida.

In [22]:
The probability of Survival after 22 days is 11.95 %

Ou seja, essa assinatura tem apenas 11.95% de probabilidade de estar ativa, o que significa que em algum momento muito próximo ela pode vir a ser cancelada.


Como podemos ver acima, usando análise de sobrevivência podemos tirar insights interessantes em relação ao nosso conjunto de dados, em especial para descobrirmos a duração das assinaturas em nossa base de dados, e estimar um tempo até o evento de churn.

Os dados utilizados refletem o comportamento de dois produtos reais, porém, que foram anonimizados por questões óbvias de NDA. Contudo nada impede a utilização e a adaptação desse código para outros experimentos. Um ponto importante em relação a essa base de dados é que como pode ser observado temos uma censura à direita muito acentuada o que limita um pouco a visão dos dados a longo prazo, principalmente se houver algum tipo de cauda longa no evento de churn.

Como coloquei no São Paulo Big Data Meetup de Março há uma série de arquiteturas que podem ser combinadas com esse tipo de análise, em especial métodos de Deep Learning que podem ser um endpoint de um pipeline de predição.

Espero que tenham gostado e quaisquer dúvidas mandem uma mensagem para flavioclesio at

PS: Agradecimentos especiais aos meus colegas e revisores Eiti Kimura, Gabriel Franco e Fernanda Eleuterio.

Churn-at-Risk: Aplicação de Survival Analysis no controle de churn de assinaturas em Telecom

Deep Dive com Gradient Boosting Machine com H2O + R (Mais Grid Search!)

Dando sequência a alguns tutoriais sobre o uso do R como linguagem de programação junto H2O como backend de processamento e memória (duas principais limitações do R) vamos falar um pouco de Gradient Boosting Machine e usar uma base de dados de crédito de um banco fictício chamado “Layman Brothers”.

Gradient Boosting Machine é um meta-algoritmo de aprendizado supervisionado que é geralmente utilizado em problemas de classificação e regressão. O principio algorítmico por trás do GBM é a produção de previsões/classificações derivadas de modelos preditivos fracos (Weak Learners), em especial árvores de decisão essas que por sua vez combinadas via ensemble learning para redução de vieses dos algoritmos.

Essas previsões são geradas através da combinação da meta-heurística de gradiente descendente para otimização paramétrica face a minimização de uma função de custo (loss function), e do Boosting que é combinação de diversos classificadores fracos (Weak Learners) em série para (ou meta-classificador) para combinação de resultados desses algoritmos.

Como podemos supor, com essa combinação heurística de algoritmos, em especial dos weak learners (que dão uma robustez substancial ao modelo) é de se esperar uma determinada insensibilidade á distribuição de cauda longa que pode ser espessa e detonar as suas previsões (e.g. distribuição da renda mundial em que poucos (20%) tem muito dinheiro e muitos (80%) tem pouco) , outliers (i.e. eventos extremos, também conhecidos como cisnes negros), além de uma boa resposta a não-linearidade. (Nota: Se você não entendeu nada do que está aqui, uma boa pedida são dois livros do Nassim Taleb que são Black Swan (A lógica do cisne negro) e Antifragile (Antifrágil)).

Como dito anteriormente, a base de dados que será usada aqui é de um banco fictício chamado “Layman Brothers”, que é uma alusão simpática ao Lehman Brothers; e o nosso objetivo é ter um sistema de crédito um pouco mais confiável do que o deles o que não é uma tarefa que demande muita inteligência ou stamina intelectual. (Nota: Essa base é originalmente do repositório do UCI, mas estou rebatizando para dar um tom cênico mais descontraído aqui no post).

A nossa base de dados de créditos tem as seguintes colunas:

  • ID: Número da transação
  • LIMIT_BAL: Crédito concedido em dólares
  • SEX: Sexo (1 = masculino; 2 = feminino).
  • EDUCATION: Nível escolar d@ cliente (1 = ensino médio; 2 = universidade; 3 = ensino superior completo; 4 = outros)
  • MARRIAGE: Estado civil (1 = casad@; 2 = solteir@; 3 = outros).
  • AGE: Idade d@ cliente
  • PAY_X: Histórico do pagamento passado. Foi rastreado o pagamento passado mensal (de abril até setembro de 2005) da seguinte forma: PAY_1 o status de repagamento do mês de setembro de 2005, PAY_2: o status do repagamento mês de agosto de 2005, etc. A escala de medida do repagamento é :-1 = Pago em dia, 1 = pago com um mês de atraso, 2 = pagamento atrasado por 2 meses, 8 = pagamento atrasado por 8 meses, etc.
  • BILL_AMTX: Montante do saldo ainda não amortizado dos meses anteriores. BILL_AMT1 = Saldo ainda não amortizado em setembro de 2005, BILL_AMT2 = saldo ainda não amortizado em agosto de 2005, etc.
  • PAY_AMTX: Montante pago anteriormente (em dólares) relativos ao mês anterior. PAY_AMT1 = valor pago em setembro de 2005, PAY_AMT2 = valor pago em agosto de 2005, etc.
  • DEFAULT: Se @ cliente deixou de pagar o empréstimo no mês seguinte.

Base de dados apresentada, vamos ao código.

Primeiramente, se você não instalou o H2O via R ou está com a versão desatualizada, é só executar esse código abaixo que ele vai remover a versão antiga, instalar todas as dependências, e instalar o H2O:

# The following two commands remove any previously installed H2O packages for R.
if ("package:h2o" %in% search()) { detach("package:h2o", unload=TRUE) }
if ("h2o" %in% rownames(installed.packages())) { remove.packages("h2o") }

# Next, we download packages that H2O depends on.
if (! ("methods" %in% rownames(installed.packages()))) { install.packages("methods") }
if (! ("statmod" %in% rownames(installed.packages()))) { install.packages("statmod") }
if (! ("stats" %in% rownames(installed.packages()))) { install.packages("stats") }
if (! ("graphics" %in% rownames(installed.packages()))) { install.packages("graphics") }
if (! ("RCurl" %in% rownames(installed.packages()))) { install.packages("RCurl") }
if (! ("jsonlite" %in% rownames(installed.packages()))) { install.packages("jsonlite") }
if (! ("tools" %in% rownames(installed.packages()))) { install.packages("tools") }
if (! ("utils" %in% rownames(installed.packages()))) { install.packages("utils") }

# Now we download, install and initialize the H2O package for R.
install.packages("h2o", type="source", repos=(c("")))

Agora vamos carregar a biblioteca e iniciar o nosso cluster (que nesse caso ainda estará no meu notebook) com o tamanho máximo de memória de 8 gigas, e vai usar todos os processadores (-1):

# Load library

# Start instance with all cores
h2o.init(nthreads = -1, max_mem_size = "8G")

# Info about cluster

# Production Cluster (Not applicable because we're using in the same machine)
#localH2O <- h2o.init(ip = '', port =54321, nthreads=-1) # Server 1
#localH2O <- h2o.init(ip = '', port =54321, nthreads=-1) # Server 2

Cluster iniciado, vamos buscar os nossos dados que estão no repositório remoto do Github e na sequência vamos carregar no nosso objeto .hex (extensão do H2O):

# URL with data
LaymanBrothersURL = ""

# Load data 
creditcard.hex = h2o.importFile(path = LaymanBrothersURL, destination_frame = "creditcard.hex")

Com os dados carregados, vamos realizar a transformação das variáveis categóricas, e em seguida vamos ver o sumário dessas variáveis:

# Convert DEFAULT, SEX, EDUCATION, MARRIAGE variables to categorical
creditcard.hex[,25] <- as.factor(creditcard.hex[,25]) # DEFAULT
creditcard.hex[,3] <- as.factor(creditcard.hex[,3]) # SEX
creditcard.hex[,4] <- as.factor(creditcard.hex[,4]) # EDUCATION
creditcard.hex[,5] <- as.factor(creditcard.hex[,5]) # MARRIAGE

# Let's see the summary

Como podemos ver pelo summary() temos algumas estatísticas descritivas básicas interessantes sobre essa base de dados, como:


  • A maioria dos empréstimos foram feitos por pessoas que se declararam do sexo feminino (60%);
  • 63% de todos os empréstimos foram feitos para a população classificada como universitária ou que tem curso superior completo;
  • Há um equilíbrio entre o estado civil em relação aos empréstimos concedidos;
  • Com um terceiro quartil de 41 e uma média e medianas bem próximas (35 e 34), podemos ver que grande parte dos empréstimos foram feitos por pessoas na idade adulta que estão na meia idade; e
  • Temos muitas pessoas que pegaram empréstimos altos (acima de 239 mil dólares), porém, a média do valor concedido é de 167 mil dólares.

Óbvio que caberiam mais algumas análises de perfil, correlações, e até mesmo alguns gráficos para exemplificar melhor a composição demográfica dessa base, mas como esse não é o objetivo desse post, fica aberto para que algum dos 5 leitores desse site blog faça isso e compartilhe.

Com essas análises feitas, vamos dividir a nossa base nos conjuntos de treinamento, teste e validação usando o comando splitFrame:

# We'll get 3 dataframes Train (60%), Test (20%) and Validation (20%)
creditcard.split = h2o.splitFrame(data = creditcard.hex
                                  ,ratios = c(0.6,0.2)
                                  ,destination_frames = c("creditcard.train.hex", "creditcard.test.hex", "creditcard.validation.hex")
                                  ,seed = 12345)

# Get the train dataframe(1st split object)
creditcard.train = creditcard.split[[1]]

# Get the test dataframe(2nd split object)
creditcard.test = creditcard.split[[2]]

# Get the validation dataframe(3rd split object)
creditcard.validation = creditcard.split[[3]]

Para checarmos a real proporção de cada base, podemos usar o comando table para ver a composição de cada base de dados (e principalmente ver se elas estão balanceadas):

# See datatables from each dataframe

# 1       0 14047
# 2       1  4030


# 1       0  4697
# 2       1  1285


# 1       0  4620
# 2       1  1321

Agora vamos criar dois objetos para passar ao nosso algoritmo: um objeto para definir quem será a nossa variável dependente (Y) e outro para definir as nossas variáveis independentes (X):

# Set dependent variable

# Set independent variables

# I intentionally removed sex variable from the model, to avoid put any gender bias inside the model. Ethics first guys! 😉

Os mais atentos podem verificar que eu removi a variável SEX. Fiz isso intencionalmente dado que não vamos colocar nenhum tipo de viés discriminatório no modelo (Atenção amigos: esse é um bom tempo para considerar seriamente essas questões de discriminação/ética em modelos de Machine Learning como etnia, gênero, etc).

Agora com esses objetos prontos, vamos treinar o nosso modelo:

# Train model
creditcard.gbm <- h2o.gbm(y = Y
                          ,x = X
                          ,training_frame = creditcard.train
                          ,validation_frame = creditcard.validation                      
                          ,ntrees = 100
                          ,seed = 12345
                          ,max_depth = 100
                          ,min_rows = 10
                          ,learn_rate = 0.2
                          ,distribution= "bernoulli"
                          ,model_id = 'gbm_layman_brothers_model'
                          ,build_tree_one_node = TRUE
                          ,balance_classes = TRUE
                          ,score_each_iteration = TRUE
                          ,ignore_const_cols = TRUE

Explicando alguns desses parâmetros:

  • x: Vetor que contém os nomes das variáveis independentes do modelo;
  • y: índice ou objeto que representa a variável dependente do modelo;
  • training frame: Um objeto de dados do H2O (H2OFrame) que contém as variáveis do modelo;
  • validation frame: Um objeto de dados do H2O (H2OFrame) que contém as variáveis do modelo para validação do modelo. Se estiver vazia os dados de treinamento são usados por padrão;
  • ntrees: Um inteiro não negativo que define o número de árvores. O valor default é 50;
  • seed: Semente dos números aleatórios a serem gerados. É usado para reprodutibilidade amostral;
  • max depth: Valor definido pelo usuário do número máximo da profundidade das árvores. O valor default é 5;
  • min rows: O número mínimo de linhas a serem designadas para cada nó terminal. O padrão é 10;
  • learn rate: Um inteiro que define a taxa de aprendizado do modelo. Vai de 0.1 até 1.0;
  • distribution: Escolhe uma distribuição de probabilidade entre AUTO, bernoulli, multinomial, gaussian, poisson, gamma ou tweedie. O default é AUTO;
  • model id: ID único que identifica o modelo. Se não especificado é gerado automaticamente;
  • build tree one node: Especifica se o modelo será processado em um nó apenas. Isso serve para evitar overhead de rede e com isso menos CPUs são usadas no processo. É ideal para pequenos datasets, e o default é FALSE;
  • balance classes: Faz o balanceamento de classes do conjunto de treinamento, caso os dados estejam com subamostragem ou desbalanceados. O default é falso;
  • score each iteration: Um binário que indica se haverá o processo de scoring durante cada interação do modelo. O default é falso; e
  • ignore const cols: Um binário que indica se colunas com constantes serão ignoradas. O Default é TRUE.

Alguns conselhos práticos de quem já sofreu (muito) na pele para parametrizar GBM que você não vai ter do seu professor na faculdade:

a) O H2O oferece e a opção validation_frame, porém, se você for mais purista o ideal é checar na etapa de prediction e ver o bias do modelo através da análise dos erros (sim gente, vai ter que rolar estatística aqui, ok?). Isso além de dar um ajuste mais fino, te dá o maior entendimento dos erros modelo. Se fosse em minas, o pessoal lá diria que isso faz bem pra saúde e forma o caráter. Faça o mesmo.;

b) Tenha bastante parcimônia para ajustar o número ideal de árvores (ntrees) dado que isso eleva demais o custo computacional (processamento + memória) do modelo. Via de regra, eu gosto de usar intervalos de 50 árvores para cada step até o limite de 300; e assim que eu chego em um meio termo eu vou ajustando na unha via grid search até chegar em uma árvore que eu tenha um bom desempenho sem overfitting. Isso é necessário pois grande parte das vezes você tem uma elevação ridícula de até 8 horas no tempo de treinamento pra ganhar no máximo 0.01 no AUC, ou uma redução de 0.005% nos falsos positivos. Em resumo: Vai com calma no ajuste. Faz bem pra saúde e forma o caráter; e além do mais economiza mais de 20 dólares na Amazon pra treinar um modelo caso você esteja usando máquinas on-demand fora da sua infra;

c) É o seed que vai garantir que os seus números estão corretos quando você for passar para alguém fazer o code review ou mesmo antes do deployment. Então use sempre que puder por questões óbvias de reprodutibilidade;

d) O parâmetro max depth costuma ser o que eu chamo de cemitério do malandro em Machine Learning. Isso devido ao fato de que qualquer iniciante em seu primeiro contato com esse parâmetro vai colocar o maior número possível em geral quase o mesmo número de instâncias da base de treinamento (isso é quando o malandro não coloca cross-validation pra coisa ficar ainda mais bonita) o que deixa a árvore mais específica e leva na maioria das vezes aquele overfittingTem iniciantes que conseguem a proeza de fazer overfitting mesmo usando max depth com leave-one-out cross validation. (Pequena dica empírica: pessoalmente eu nunca consegui resultados bacanas com uma profundidade de níveis que excedam 0.005% do número de registros no conjunto de treinamento (100/((30000/100)*70 =0.005%). Ainda estou tentando saber se isso está correto ou não, mas ao menos pra mim funciona bem;

e)  Quanto menor o valor do min rows, mais específica será a árvore e pode ocorrer que ela generalize menos. Por isso muita parcimônia com esse parâmetro;

f) Desnecessário dizer que um número muito pequeno pode influenciar no tempo de processamento e convergência do modelo, e um número alto pode cair em um mínimo local e estragar todo o trabalho. Dica prática: tá com pouco tempo? Vai de 0.35 até 0.75 com incremento de 0.1. Tá com tempo de sobra? Vai de 0.1 até 0.5 com incremento de 0.03;

g) Realmente vale a pena gastar um pouco de neurônios para entender melhor as distribuições de probabilidade (distribution) para escolher a correta. Se você não tiver tempo, escolha a AUTO e seja feliz;

h) A não ser que você esteja enfrentando uma situação de concorrência de rede e de processamento, o parâmetro build tree one node sempre deve estar desligado;

i) Se você está usando o parâmetro balance classes significa que o seu trabalho de amostragem está um lixo e você precisa da ferramenta pra fazer algo básico pode não ser o mais correto. Eu recomendo fortemente uma seriedade no processo de amostragem que é o coração de qualquer treinamento de machine learning. Caso sejam situações amostrais muito esquisitas (e.g. modelagem de sistemas de combate á fraudes, classificador de reclamações em Call Center, et cetera) ou por falta de tempo, vale a pena usar esse parâmetro (Dica prática: caso haja uma situação de desbalanceamento muito grave de classes (algo na proporção 9:1) o ideal é esquecer as outras métricas de avaliação de modelos e ir direto para o coeficiente de matthews que é bem mais consistente para lidar com esse tipo de caso);

j) Se você está usando o parâmetro ignore const cols é porque o seu trabalho de pré-processamento (Feature Extraction e Feature Engineering) está um lixo pode não estar sendo o melhor.

Modelo treinado e parâmetros explicados, vamos ver a performance do modelo usando os dados de validação:

# See algo performance
h2o.performance(creditcard.gbm, newdata = creditcard.validation)

# H2OBinomialMetrics: gbm

# MSE:  0.1648487
# RMSE:  0.4060157
# LogLoss:  0.8160863
# Mean Per-Class Error:  0.3155595
# AUC:  0.7484422
# Gini:  0.4968843

# Confusion Matrix for F1-optimal threshold:
#   0    1    Error        Rate
# 0      3988  632 0.136797   =632/4620
# 1       653  668 0.494322   =653/1321
# Totals 4641 1300 0.216294  =1285/5941

# We have an AUC of 74,84%, not so bad!

Com esse modelo tivemos um AUC de 74,84%. Razoável, considerando que usamos um conjunto de parametrizações simples.

A seguir, vamos conferir a importância de cada uma de nossas variáveis:

# Variable importance
imp <- h2o.varimp(creditcard.gbm)

head(imp, 20)

# Variable Importances: 
#   variable relative_importance scaled_importance percentage
# 1  EDUCATION        17617.437500          1.000000   0.380798
# 2   MARRIAGE         9897.513672          0.561802   0.213933
# 3      PAY_0         3634.417480          0.206297   0.078557
# 4        AGE         2100.291992          0.119217   0.045397
# 5  LIMIT_BAL         1852.831787          0.105170   0.040049
# 6  BILL_AMT1         1236.516602          0.070187   0.026727
# 7   PAY_AMT5         1018.286499          0.057800   0.022010
# 8  BILL_AMT3          984.673889          0.055892   0.021284
# 9  BILL_AMT2          860.909119          0.048867   0.018608
# 10  PAY_AMT6          856.006531          0.048589   0.018502
# 11  PAY_AMT1          828.846252          0.047047   0.017915
# 12 BILL_AMT6          823.107605          0.046721   0.017791
# 13 BILL_AMT4          809.641785          0.045957   0.017500
# 14  PAY_AMT4          771.504272          0.043792   0.016676
# 15  PAY_AMT3          746.101196          0.042350   0.016127
# 16 BILL_AMT5          723.759521          0.041082   0.015644
# 17     PAY_3          457.857758          0.025989   0.009897
# 18     PAY_5          298.554657          0.016947   0.006453
# 19     PAY_4          268.133453          0.015220   0.005796
# 20     PAY_2          249.107925          0.014140   0.005384

Nesse modelo podemos ver que o nível educacional tem um papel essencial na definição de quem vai entrar em default (38%), seguindo do estado civil (21%) e fechando com o pagamento anterior relativo ao mês de setembro de 2008 (7%) e da idade do tomador de crédito e o saldo emprestado (4%).

Em outras palavras: essas variáveis acima respondem por 74% do comportamento de crédito.

Com isso algumas questões hipóteses (Hx) e ações (Ax) podem ser tomadas pelo Layman Brothers:

H1: O nível educacional está muito relacionado com o default,  isso acontece de forma positiva ou não em relação à inadimplência?

H2: Será que universitários que tradicionalmente são pessoas com menos poder aquisitivo tem maiores dificuldades (ou facilidades) para o pagamento?

H3: De que forma o estado civil influencia na capacidade de pagamento do crédito emprestado?

H4: Porque o saldo não amortizado exerce efeito tão grande em relação às outras variáveis financeiras?

H5: Porque a pontualidade no pagamento não é tão determinante, com exceção da primeira parcela?

H6: O perfil educacional influencia o quanto em relação à capacidade de pagamento?

A1: De acordo com a escolaridade, ter diferentes taxas de juros para empréstimos.

A2: Ter ações de cobrança efetivas/intensas já no primeiro mês de atraso.

A3: Ter linhas de crédito mais específicas para cada perfil educacional com taxas e saldos correspondentes ao risco de default.

A4: Entender e criar linhas de financiamento de acordo com cada objetivo de acordo com o estado civil (e.g. entender se o gasto é para investimento (voltado para a geração de mais receita como cursos, maquinário, ou outros fatores que aumentem a produtividade; ou para despesas como consumo, contas de inúmeras naturezas, outros empréstimos, et cetera) .

Adiante, podemos agora usar o nosso modelo treinado para fazer previsões:

# Predict using GLM model
pred = h2o.predict(object = creditcard.gbm, newdata = creditcard.test)

# See predictions
head(pred, 5)

# predict        p0           p1
# 1       0 0.9990856 0.0009144487
# 2       0 0.9945627 0.0054373206
# 3       0 0.9997726 0.0002273775
# 4       0 0.9968271 0.0031728833
# 5       0 0.9991758 0.0008242144

Agora, vamos para um ajuste mais fino no nosso modelo com o objetivo de melhorar o nosso AUC (que é atualmente de 74,84%), e para isso vamos usar Grid Search.

Primeiramente vamos gerar uma lista de valores para os nossos hiper-parâmetros (hyper parameters) do modelo GBM. Os parâmetros que vamos usar serão ntrees (número de árvores), max_depth (profundidade das árvores) e learn_rate (taxa de aprendizado). Após isso vamos jogar dentro de uma meta lista que vamos usar para ajustar o nosso objeto de grid.

# Set hyparameters (Did not work using sequence. :o( )
ntrees_list <- list(50,100,150,200)

max_depth_list <- list(1,2,3,4,5,6,7,8,9,10)

learnrate_list <- list(.10,.20,.30,.40,.50,.60,.70,.80,.90)
# Full list of hyper parameters that will be used
hyper_parameters <- list(ntrees = ntrees_list
                         ,max_depth = max_depth_list
                         ,learn_rate = learnrate_list)

# See hyparameters lists

Ou seja, teremos uma combinação com 50, 100, 150 e 200 árvores, níveis de profundidade da árvore indo de 1 até 10 e taxa de aprendizado indo de 0.10 até 0.90.

Uma pequena experiência da trincheira deste escriba que não foi muito inteligente é ter uma boa combinação de números de parâmetros na meta lista em relação com a capacidade de processamento disponível para fazer o treinamento.

Isso se faz necessário pois como abaixo vamos usar a estratégia cartesiana para o nosso critério de busca (i.e. vamos usar todas as combinações paramétricas possíveis) vamos ter o seguinte cenário:

ntrees = 4
max_depth = 10
learn_rate = 9

Logo teremos 4 * 10 * 9 = 360 modelos/combinações!

Ou seja: Pode levar bastante tempo para processar (no meu caso levou 11m34min pra acabar, e houve uma porção de erros do H2O por incapacidade de processamento).

Após o processamento do grid vamos ordenar os modelos do melhor para o pior usando o AUC:

# sort the grid models by decreasing AUC
sortedGrid <- h2o.getGrid("depth_grid", sort_by="auc", decreasing = TRUE)    
# Let's see our models

# H2O Grid Details
# ================
# Grid ID: depth_grid 
# Used hyper parameters: 
# -  learn_rate 
# -  max_depth 
# -  ntrees 
# Number of models: 380 
# Number of failed models: 2940 

# Hyper-Parameter Search Summary: ordered by decreasing auc
# learn_rate max_depth ntrees            model_ids                auc
# 1        0.1         6    100 depth_grid_model_200 0.7811807105334736
# 2        0.1         6     50   depth_grid_model_5 0.7811440893197138
# 3        0.2         3    150 depth_grid_model_264 0.7809025695475355
# 4        0.2         3    100 depth_grid_model_174  0.780834324645831
# 5        0.1         6    200 depth_grid_model_380 0.7808292451933633

Agora, vamos pegar o melhor modelo (com menor AUC) e vamos ver algumas das suas características:

# Summary

# Model Details:
# ==============
# H2OBinomialModel: gbm
# Model Key:  depth_grid_model_200 
# Model Summary: 
#   number_of_trees number_of_internal_trees model_size_in_bytes min_depth max_depth mean_depth
# 1             100                      100               52783         6         6    6.00000
# min_leaves max_leaves mean_leaves
# 1         12         56    36.93000

# H2OBinomialMetrics: gbm
# ** Reported on training data. **
# MSE:  0.1189855
# RMSE:  0.3449427
# LogLoss:  0.3860698
# Mean Per-Class Error:  0.2593832
# AUC:  0.8371354
# Gini:  0.6742709

# Confusion Matrix for F1-optimal threshold:
# 0    1    Error         Rate
# 0      12424 1623 0.115541  =1623/14047
# 1       1625 2405 0.403226   =1625/4030
# Totals 14049 4028 0.179676  =3248/18077

Esse nosso modelo tem 100 árvores, uma profundidade de 6 níveis, e em média 37 instâncias em cada nó folha.

Como podemos ver tivemos um AUC de 83,71%, ou 11% de melhoria em comparação com o antigo AUC que foi de 74,84% em menos de 12 minutos.

Um fato curioso é que olhando a importância das variáveis novamente com esse modelo temos os seguintes resultados:

# Variable importance (again...)
imp2 <- h2o.varimp(best_glm)

head(imp2, 20)

# Variable Importances: 
#   variable relative_importance scaled_importance percentage
# 1      PAY_0         2040.270508          1.000000   0.358878
# 2      PAY_2          902.637390          0.442411   0.158772
# 3  LIMIT_BAL          385.425659          0.188909   0.067795
# 4        AGE          274.609589          0.134595   0.048303
# 5  BILL_AMT1          209.715469          0.102788   0.036888
# 6      PAY_3          168.518372          0.082596   0.029642
# 7  EDUCATION          150.365280          0.073699   0.026449
# 8  BILL_AMT2          146.754837          0.071929   0.025814
# 9      PAY_5          139.303482          0.068277   0.024503
# 10  PAY_AMT5          139.206543          0.068229   0.024486
# 11 BILL_AMT5          133.963348          0.065660   0.023564
# 12     PAY_4          124.926552          0.061230   0.021974
# 13  PAY_AMT6          123.267151          0.060417   0.021682
# 14 BILL_AMT6          114.012253          0.055881   0.020054
# 15  PAY_AMT1          112.402290          0.055092   0.019771
# 16     PAY_6          108.483795          0.053171   0.019082
# 17 BILL_AMT3          103.207893          0.050585   0.018154
# 18  PAY_AMT3           97.335411          0.047707   0.017121
# 19 BILL_AMT4           90.403320          0.044309   0.015902
# 20  MARRIAGE           61.917801          0.030348   0.010891

Ou seja, se antigamente o nível educacional e o estado civil tiveram uma participação importante, nesse modelo (com melhor AUC) a pontualidade, o montante de crédito concedido e a idade exercem mais influência.

Com esse melhor modelo, podemos fazer as nossas previsões e salvar em um arquivo .csv para upload em algum sistema ou isso pode ser feito via API via requisição.

# Get model and put inside a object
model = best_glm

# Prediction using the best model
pred2 = h2o.predict(object = model, newdata = creditcard.validation)

# Frame with predictions
dataset_pred =

# Write a csv file
write.csv(dataset_pred, file = "predictions.csv", row.names=TRUE)

 E após finalizado todo o trabalho, podemos desligar o nosso cluster:

# Shutdown the cluster 

# Are you sure you want to shutdown the H2O instance running at http://localhost:54321/ (Y/N)? Y
# [1] TRUE

Bem pessoal como vocês podem ver, usar um modelo usando Gradient Boosting Machine no R não é nenhum bicho de 7 cabeças no H2O, basta um pouquinho de parcimônia na parametrização que tudo dá certo.

Se tiverem dúvidas deixem o seu comentário inteligente e educado aqui nos comentários ou me mandem por e-mail.

Forte abraço!


Deep Dive com Gradient Boosting Machine com H2O + R (Mais Grid Search!)