Funções com Multiprocessing para processamento de textos

Quem acompanhou o post A small journey in the valley of Natural Language Processing and Text Pre-Processing for German language acompanhou um pouco dos desafios de modelar um classificador de textos em alemão.

No entanto uma coisa que me salvou na parte de pre-processing foi que eu praticamente usei o multiprocessing para paralelizar o pré-processamento na coluna de texto e isso me salvou um tempo incrível (relembrando: eu tinha 1+ milhão de registros de texto, com 250 palavras média por registro (com um desvio padrão de 700, tudo isso usando biblioteca interna).

import time
import numpy as np
import pandas as pd
import nlp_pre_processing
# An internal NLP lib to process text
nlp = nlp_pre_processing_library.NLPPreprocessor()
# Multiprocessing library that uses pool
# to distribute the task for all processors
from multiprocessing import Pool
print(f'Start processing…: {(time.time() start_time)}')
# Tracking the time
start_time = time.time()
# Number of partitions that
# the Pandas Dataframe will be
# splited to parallel processing
num_partitions = 20
# Number of cores that will be used
# more it's better
num_cores = 16
print(f'Partition Number: {num_partitions} – Number of Cores: {num_cores}…')
def main_process_pipeline(df, func):
Function that will split the dataframe
and process all those parts in a n number
of processors
df (Pandas dataframe): Dataframe that will be splited
func (function): Python function that will be executed in parallel
df: Dataframe with all parts concatenated after the function be applied
df_split = np.array_split(df, num_partitions)
pool = Pool(num_cores)
df = pd.concat(, df_split))
return df
def pre_process_wrapper(df):
""" Will take the Dataframe and apply a function using lambda"""
df['text'] = df['text'].apply(lambda text: nlp.pre_processing_pipeline(text))
return df
# Unite the Dataframe and the Wrapper
processed_df = main_process_pipeline(df, pre_process_wrapper)
print(f'Processing finished in seconds: {(time.time() start_time)}')
view raw hosted with ❤ by GitHub

É isso: Simples e tranquilo.

Funções com Multiprocessing para processamento de textos

Security in Machine Learning Engineering: A white-box attack and simple countermeasures

Some weeks ago during a security training for developers provided by Marcus from Hackmanit (by the way, it’s a very good course that goes in some topics since web development until vulnerabilities of NoSQL and some defensive coding) we discussed about some white box attacks in web applications (e.g.attacks where the offender has internal access in the object) I got a bit curious to check if there’s some similar vulnerabilities in ML models. 

After running a simple script based in [1],[2],[3] using Scikit-Learn, I noticed there’s some latent vulnerabilities not only in terms of objects but also in regarding to have a proper security mindset when we’re developing ML models. 

But first let’s check a simple example.

A white-box attack in a Scikit-Learn random forest object

I have a dataset called Layman Brothers that consists in a database of loans that I did grab from internet (if someone knows the authors let me know to give the credit) that contains records regarding consumers of a bank that according some variables indicates whether the consumer defaulted or not. This is a plain vanilla case of classification and for a matter of simplicity I used a Random Forest to generate a classification model. 

The main points in this post it’s check what kind of information the Scikit-Learn object (model) reveals in a white-box attack and raises some simple countermeasures to reduce the attack surface in those models.

After ran the classifier, I serialized the Random Forest model using Pickle. The model has the following performance against the test set:

# Accuracy: 0.81
# status
# 0 8071
# 1 929
view raw hosted with ❤ by GitHub

Keep attention in those numbers because we’re going to talk about them later on in this post. 

In a quick search in internet the majority of applications that uses Scikit-Learn for production environments deals with a pickled (or serialized model in Joblib) object that it’s hosted in a machine or S3 bucket and an API REST take care to do the servicing of this model. The API receives some parameters in the request and bring back some response (the prediction result). 

In our case, the response value based on the independent variables (loan features) will be defaulted {1}or not {0}. Quite straightforward.  

Having access in the Scikit-Learn object I noticed that the object discloses valuable pieces of information that in the hands of a white-box attacker could be potentially very damaging for a company. 

Loading the pickled object, we can check all classes contained in a model:

So, we have a model with 2 possible outcomes, {0} and {1}. From the perspective of an attacker we can infer that this model has some binary decision containing a yes {1}or no {0}decision. 

I need to confess that I expected to have only a read accessin this object (because the Scikit-Learn documentation gives it for grant), but I got surprised when I discovered that I can write in the objecti.e. overriding the training object completely. I made that using the following snippet:

# Load model from Pickle
model_rf_reload_pkl = pickle.load(open(‘model_rf.pkl’, ‘rb’))
# Displays prediction classes
# >>> array([0, 1])

One can noticed that with this command I changed all the possible classes of the model only using a single numpy array and hereafter this model will contain only the outcome {1}.

Just for a matter of exemplification I ran the same function against the test dataset to get the results and I got the following surprise:

# Actual against test set
# Accuracy: 0.2238888888888889
# status
# 1 9000
# Previous against test set
# Accuracy: 0.8153333333333334
# status
# 0 8071
# 1 929
view raw hosted with ❤ by GitHub

In this simple example we moved more than 8k records to the wrong class. It’s unnecessary to say how damaging this could be in production in some critical domain like that. 

If we do a simple mental exercise, where this object could be a credit score application, or some classifier for in a medical domain, or some pre-order of some market stocks; we can see that it brings a very cold reality that we’re not close to be safe doing the traditional ML using the popular tools. 

In the exact moment that we Machine Learning Engineers or Data Scientists just run some scripts without even think in terms of vulnerabilities and security, we’re exposing our employers, business and exposing ourselves in such liability risk that can cause a high and unnecessary damage because of the lack of a better security thinking in ML models/objects. 

After that, I opened an issue/question in Scikit-Learn project to check the main reason why this type of modification it’s possible. Maybe I was missing something that was thought by the developers during the implementation phase. My issue in the project can be seeing below:

And I got the following response:

Until the day when this post was published there’s no answer for my last question about this potential vulnerability in a parameter that should not be changed after model training.

This is a simple white-box attack that can interfere directly in the model object itself. Now let’s pretend that we’re not an attacker in the object, but we want to explore other attack surfaces and check which valuable information those models can give for us. 

Models revealing more than expected

Using the same object, I’ll explore the same attributes that is given in the docs to check if we’re able to fetch more information from the model and how this information can be potentially useful.

First, I’ll try to see the number of estimators:

print(fNumber of Estimators: {len(model_rf_reload_pkl.estimators_)}’)
# >>> Number of Estimators: 10
view raw hosted with ❤ by GitHub

This Random Forest training should be not so complex, because we have only 10 estimators (i.e. 10 different trees) and grab all the complexities of this Random Forest won’t be so hard to a mildly motivated attacker. 

I’ll grab a single estimator to perform a quick assessment in the tree (estimator) complexity:

# >>> DecisionTreeClassifier(class_weight=None, criterion='gini',
# max_depth=5,max_features='auto',
# max_leaf_nodes=5, min_impurity_decrease=0.0,
# min_impurity_split=None, min_samples_leaf=100,
# min_samples_split=2, min_weight_fraction_leaf=0.0,
# presort=False, random_state=1201263687,
# splitter='best')
view raw hosted with ❤ by GitHub

Then this tree it’s not using a class_weight to perform training adjustments if there’s some unbalance in the dataset. As an attacker with this piece of information, I know that if I want to perform attacks in this model, I need to be aware to alternate the classes during my requests. 

It means that if I get a single positive result, I can explore in alternate ways without being detected as the requests are following by a non-weighted distribution.

Moving forward we can see that this tree has only 5 levels of depth (max_depth) with a maximum 5 leaf nodes (max_leaf_nodes) with a minimum of 100 records per leaf (min_samples_leaf).

It means that even with such depth I can see that this model can concentrate a huge amount of cases in some leaf nodes (i.e. low depth with limited number of leaf nodes). As an attacker maybe I don’t could not have access in the number of transactions that Layman Brothers used in the training, but I know that the tree it’s simple and it’s not so deep. 

In other words, it means that my search space in terms of parameters won’t be so hard because with a small number of combinations I can easily get a single path from the root until the leaf node and explore it.

As an attacker I would like to know how many features one estimator contains. The point here is if I get the features and their importance, I can prune my search space and concentrate attack efforts only in the meaningful features. To know how many features one estimator contains, I need to run the follow snippet:

# Extract single tree
estimator = model_rf_reload_pkl.estimators_[5]
print(fNumber of Features: {estimator.n_features_}’)
# >> Number of Features: 9
view raw hosted with ❤ by GitHub

As we can see, we have only 9 parameters that were used in this model. Then, my job as an attacker could not be better. This is a model of dreams for an attacker. 

But 9 parameters can be big enough in terms of a search space. To prune out some non-fruitful attacks of my search space, I’ll try to check which variables are relevant to the model. With this information I can reduce my search space and go directly in the relevant part of the attack surface. For that let’s run the following snippet:

features_list = [str(x + 0) for x in range(estimator.n_features_)]
# >>> ['0', '1', '2', '3', '4', '5', '6', '7', '8']
importances = estimator.feature_importances_
indices = np.argsort(importances)
plt.title(‘Feature Importances’)
plt.barh(range(len(indices)), importances[indices], color=b’, align=center’)
plt.yticks(range(len(indices)), indices)
plt.xlabel(‘Relative Importance’)

Let me explain: I did grab the number of the features and put an id for each one and after that I checked the relative importance using np.argsort()to assess the importance of those variables. 

As we can see I need only to concentrate my attack in the features in the position [4][3]and [5]. This will reduce my work in tons because I can discard other 6 variables and the only thing that I need to do it’s just tweak 3 knobs in the model. 

But trying something without a proper warmup could be costly for me as an attacker and I can lose the momentum to perform this attack (e.g. the ML team can update the model, someone can identify the threat and remove the model and rollback an old artifact, etc).

To solve my problem, I’ll check the values of one of trees and use it as a warmup before to do the attack. To check those values, I’ll run the following code that will generate the complete tree:

from sklearn.tree import export_graphviz
# Export as dot file
feature_names = features_list,
rounded = True, proportion = False,
precision = 2, filled = True)
# Convert to png using system command (requires Graphviz)
from subprocess import call
call([‘dot’, ‘Tpng’, ‘’, ‘o’, ‘tree.png’, ‘Gdpi=600’])
# Display in jupyter notebook
from IPython.display import Image
Image(filename =tree.png’)
view raw hosted with ❤ by GitHub

Looking the tree graph, it’s clearer that to have our Loan in Layman Brothers always with false {0}in the defaultvariable we need to tweak the values in {Feature 3<=1.5 && Feature 5<=-0.5 && Feature 4<=1.5}.

Doing a small recap in terms of what we discovered: (i) we have the complete model structure, (ii) we know which features it’s important or not, (iii) we know the complexity of the tree, (iv) which features to tweak and (v) the complete path how to get our Loan in Layman Brothers bank always approved and game the system. 

With this information until the moment that the ML Team changes the model, as an attacker I can explore the model using all information contained in a single object. 

As discussed before this is a simple white-box approach that takes in consideration the access in the object. 

The point that I want to make here it’s how a single object can disclose a lot for a potential attacker and ML Engineers must be aware about it.

Some practical model security countermeasures

There are some practical countermeasures that can take place to reduce the surface attack or at least make it harder for an attacker. This list it’s not exhaustive but the idea here is give some practical advice for ML Engineers of what can be done in terms of security and how to incorporate that in ML Production Pipelines. Some of those countermeasures can be:

  • Consistency Checks on the object/model: If using a model/object it’s unavoidable, one thing that can be done is load this object and using some routine check (i) the value of some attributes (e.g. number of classes in the model, specific values of the features, tree structure, etc), (ii) get the model accuracy against some holdout dataset (e.g. using a holdout dataset that has 85% of accuracyand raise an error in any different value), (iii) object size and last modification.  

  • Use “false features: False features will be some information that can be solicited in the request but in reality, won’t be used in the model. The objective here it’s to increase the complexity for an attacker in terms of search space (e.g.a model can use only 9 features, but the API will request 25 features (14 false features)). 

  • Model requests monitoring: Some tactics can be since monitoring IP requests in API, cross-check requests based in some patterns in values, time intervals between requests.

  • Incorporate consistency checks in CI/CD/CT/CC: CI/CD it’s a common term in ML Engineering but I would like to barely scratch two concepts that is Continuous Training (CT) and Continuous Consistency (CC). In Continuous Training the model will have some constant routine of training during some period of time in the way that using the same data and model building parameters the model will always produce the same results.  In Continuous Consistency it’s an additional checking layer on top of CI/CD to assess the consistency of all parameters and all data contained in ML objects/models. In CC if any attribute value got different from the values provided by the CT, the pipeline will break, and someone should need to check which attribute it’s inconsistent and investigate the root cause of the difference.

  • Avoid expose pickled models in any filesystem (e.g. S3) where someone can have access:  As we saw before if someone got access in the ML model/objects, it’s quite easy to perform some white-box attacks, object access will reduce the exposure/attack surface.

  • If possible, encapsulates the coefficients in the application code and make it private: The heart of those vulnerabilities it’s in the ML object access. Remove those objects and incorporates only model coefficients and/or rules in the code (using private classes) can be a good way out to disclose less information about the model.

  • Incorporate the concept of Continuous Training in ML Pipeline: The trick here it’s to change the model frequently to confuse potential attackers (e.g.different positions model features between the API and ML object, check the reproducibility of results (e.g.accuracy, recall, F1 Score, etc) in the pipeline).

  • Use heuristics and rules to prune edge cases: Attackers likes to start test their search space using some edge (absurd) cases and see if the model gives some exploitable results and fine tuning on top of that. Some heuristics and/or rules in the API side can catch those cases and throw a cryptic error to the attacker make their job quite harder.

  • Talk its silver, silence its gold: One of the things that I learned in security it’s less you talk about what you’re doing in production, less you give free information to attackers. This can be harsh but its a basic countermeasure regarding social engineering. I saw in several conferences people giving details about the training set in images like image sizing in the training, augmentation strategies, pre-checks in API side, and even disclosing the lack of strategies to deal with adversarial examples. This information itself can be very useful to attackers in order to give a clear perspective of model limitations. If you want to talk about your solution talk more about reasons (why) and less in terms of implementation (how). Telling in Social Media that I keep all my money underneath my pillow, my door has only a single lock and I’ll arrive from work only after 8PM do not make my house safer. Remember: Less information = less exposure. 


I’ll try to post some countermeasures in a future post, but I hope that as least this post can shed some light in the ML Security. 

There’s a saying in aviation culture (a good example of industry security-oriented) that means “the price of safety it’s the eternal vigilance” and I hope that hereafter more ML Engineers and Data Scientists can be more vigilant about ML and security. 

As usual, all codes and notebooks are in Github.

Security in Machine Learning Engineering: A white-box attack and simple countermeasures

Instituto Mises Brasil – Uma análise editorial usando Natural Language Processing

Este post será o primeiro de uma pequena série onde eu fiz algumas analises sobre a mudança editorial do Instituto Mises Brasil usando um pouco de Natural Language Processing (NLP) e Latent Dirichlet Allocation (LDA). Este texto inicial é apenas a descrição do repositório em que os dados e os scripts estão armazenados. Tudo é livre e pode ser copiado e contestado sem nenhum tipo de restrição.

O endereço do repositório é:

Sabe quando algo muda na redação de algum jornal, revista ou algum meio de comunicação, mas você não sabe o que é? Pois bem, eu fiquei com a mesma dúvida e resolvi usar algumas ferramentas para validar se houve uma mudança ou não.


Tem algum tempo que eu venho acompanhando a política brasileira da perspectiva de ensaístas de vertentes ligadas ao libertarianismo, anarcocapitalismo, secessão, autopropriedade e assuntos correlatos; e um fato que me chamou bastante a atenção foi a mudança editorial que está acontecendo de forma lenta em um dos principais Think Tanks liberais do Brasil que é o Instituto Mises Brasil (IMB).

Para quem não sabe, em meados de 2015 houve uma ruptura no núcleo do IMB em que de um lado ficou o Presidente do IMB (Hélio Beltrão) e do outro ficaram o Chiocca Brothers (Fernando, Cristiano, Roberto) que na sequência criaram o Instituto Rothbard. O motivo dessa ruptura foi devido a divergências relativas a artigos ligados a secessão.

E devido a esse processo de ruptura que eu penso que houve essa transição do IMB para uma linha editorial mais leve no que diz respeito a assuntos ligados à liberdade o que contrário as ideias do próprio Ludwig von Mises.

Qual o motivo desse repositório?

O principal motivo é fazer uma análise de dados simples usando Natural Processing Language (NLP) em todos textos do Mises Brasil para validar uma hipótese que é:

  • Hipótese [0]: Houve uma mudança editorial no Instituto Mises Brasil em que assuntos ligados ao austrolibertarianismo, liberdade, ética, e secessão e outras questões relacionadas deram espaço para temas efêmeros como financismo, burocracia e principalmente política.

Caso a resposta da H0 seja positiva, eu vou tentar chegar nas repostas das seguintes perguntas que são:

  • Pergunta [1]: Caso H0 for verdadeira, assuntos ligados ao austrolibertarianismo como praxeologia, fim do estatismo, ética argumentativa, e secessão estão sendo deixados de lado em termos editoriais?
  • Pergunta [2]: O Instituto Mises Brasil está se tornando em termos editoriais mais liberal-mainstream do que libertário?
  • Pergunta [3]: Houve uma mudança em relação ao grupo de assuntos que são tratados ao longo do tempo como também a mudança do espectro de assuntos dos articulistas presentes?


Tudo isso foi gerado em um MacMini com Python 3.6, mas também pode ser executado em computadores com Linux com a pré-instalação das seguintes bibliotecas:

$ pip install numpy==1.17.2
$ pip install pandas==0.25.1
$ pip install requests==2.22.0
$ pip install spacy==2.2.1
$ pip install beautifulsoup4==4.8.1
$ pip install bs4==0.0.1
$ python -m spacy download pt_core_news_sm

Sinceramente: Usem o R para geração dos seus próprios gráficos. Eu adoro o Seaborn e o Matplotlib para geração de gráficos, mas nesse sentido o R é muito mais flexível e precisa de bem menos “hacking” para fazer as coisas ficarem legais.

Extração de dados

A base que está no repositório foi gerada em 16.10.2019 para fins de congelar a análise e deixar a mesma com um grau maior de replicabilidade.

A extração busca todos os textos, independente de ser artigo do blog ou post da página principal. Isso ocorre devido ao fato de que não há uma divisão das URLs que faça essa distinção e algumas vezes temos artigos do blog que viram posts na página principal.

Outro ponto é que deve ser mencionado é que o Leandro Roque é o principal tradutor/ ensaísta do site e alguns posts de tradução são assinados por ele (o que é correto). Isso leva a dois efeitos que são 1) ele e muito profícuo com o fluxo de artigos no site e definitivamente isso distorce as estatísticas individuais dele como ensaísta e 2) por causa das traduções ele tem um espectro de assuntos bem mais diversos do que os outros autores, e isso tem que ser considerado quando analisarmos os assuntos os quais ele mais escreve. Pessoalmente eu desconsideraria ele de todas as análises dado esses dois pontos colocados. Mas aí vai de cada um.

Para quem quiser gerar uma base de dados nova com os dados até a presente data basta executar o comando abaixo no terminal:

$ python3

Ao final da execução vão aparecer as seguintes informações:

Fetching Time: 00:16:52
Articles fetched: 2855

Avisos gerais

Essa análise é apenas para fins educacionais. É obvio que uma análise editorial que envolva questões linguísticas/semânticas é algo muito complexo até mesmo para nós seres humanos, e colocar que uma máquina consiga fazer isso é algo que não tem muito sentido dado a natureza da complexidade da linguagem e as suas nuances.

Esse repositório como também a análise não tem a menor pretensão de ser algo “cientifico”. Isso significa que não haverá elementos de linguística cognitiva, linguística computacional, análise do discurso ou ciências similares. Esse repositório traz muitas visões pessoais e observações que porventura usa alguns dados e alguns scripts.

Distribuição e usos

Todos os dados, scripts, gráficos podem ser usados livremente sem nenhum tipo de restrição. Quem puder ajudar faz um hyperlink para o meu site/blog ou pode citar academicamente que vai ajudar bastante.

Eu não sou dono dos direitos dos textos do Instituto Mises Brasil e aqui tem apenas uma compilação dos dados extraídos do site, dados estes públicos e que podem ser extraídos por qualquer pessoa.

Garantias, erros e afins

Não existe garantia nestas análises, gráficos, dados e scripts e o uso está por conta de quem usar. Vão ter muitos erros (principalmente gramaticais, sintáticos e semânticos) e na medida que forem acontecendo podem abrir um Pull Request ou me mandar um e-mail que eu vou ajustando. Porém, como eu vou escrevendo na velocidade dos meus pensamentos nem sempre o sistema responsável pela correção sintática vai funcionar bem.

Instituto Mises Brasil – Uma análise editorial usando Natural Language Processing

O Campeonato Brasileiro está ficando mais injusto?

Uma análise exploratória usando o Coeficiente de Gini

Todos os dados e a análise completa pode ser encontrada no repositório brasileirao-gini. Lá tem todas as instruções para executar passo a passo a análise.


Uma das coisas boas da era digital é que com o uso do WhatsApp podemos ter as nossas mesas redondas virtuais com os nossos amigos não importa o quanto estamos longe estamos deles. E dentro de um dessas resenhas virtuais estavamos discutindo o ótimo trabalho do treinador do Flamengo na atual temporada (i.e. aspectos relativos de como ele mudou positivamente o time).

Porém, apareceu um tópico relevante durante o debate que foi a hipótese que na medida em que entra mais dinheiro no campeonato e nos times maiores o Brasileirão vai ficando cada vez mais injusto com um conjunto de times ganhando muitos jogos e os times pequenos ficarem relegados a meros coadjuvantes dentro da liga.

Em outras palavras, isso significa que o campeonato sempre é disputado pelos mesmos clubes com poderio econômico e por causa dessa disparidade a competitividade poderia estar diminuindo ao longo do tempo.

Isso levantou a questão dentro do nosso grupo que foi: O brasileirão está ficando mais injusto ao longo do tempo?

E é essa pergunta que eu vou tentar responder no final desse post.

Como checar se há uma desigualdade estrutural no Campeonato Brasileiro?

E para responder inicialmente essa pregunta eu vou usar o Coeficiente de Gini que é uma métrica para medir a dispersão estatística que inicialmente foi criada para medir a distribuição de renda e riqueza entre países e é amplamente usada em economia como um importante indicador de monitoramento dessas questões.

O Coeficiente de Gini é usado como medida relativa para benchmark e monitoramento de desigualdade e pobreza através de diversos países e serve como base de analise e desenvolvimento de politicas publicas como você pode conferir nos trabalhos de SobottkaMoreira e autores e pelo Instituto de Pesquisa Econômica Aplicada o IPEA no trabalho de Barros e autores.

Como este post não é para falar em profundidade sobre este indicador eu sugiro a leitura do trabalho original de Conrrado Gini chamado Variabilità e Mutabilità ou este trabalho sobre como esse indicador foi usado de parâmetro para alguns algumas analises de bem estar social no Brasil.

Isto é, eu vou medir de forma simples a variância da distribuição dos pontos dentro de cada uma das edições do Campeonato Brasileiro (que eu vou chamar daqui em diante de Brasileirão) para verificar se há uma desigualdade latente estrutural dentro dos campeonatos e ao longo dos anos.

Algumas considerações prévias sobre algumas limitações sobre a forma de medir essa desigualdade

De acordo com o que coloquei anteriormente podemos fazer de maneira bem simples a mensuração de fatores importantes como desigualdade econômica através da renda e/ou riqueza usando o Coeficiente de Gini, e aqui eu vou usar uma alegoria simples de como eu vou aplicar isso nos dados do Brasileirão.

Se esse índice serve de alguma maneira medir a desigualdade entre países considerando a sua riqueza ou mesmo a sua renda, trazendo para o mundo do futebol podemos usar como representação dos dados os pontos ganhos de um time dentro de uma temporada como se fosse a renda e aplicar o Coeficiente de Gini regular não somente dentro de um ano especifico da liga mas também monitorar essa desigualdade ao longo do tempo.

Para isso eu vou usar a base de dados de todos os resultados do Brasileirão desde 2003 até 2018 extraídos da Wikipédia. E aqui eu tenho que fazer duas considerações que são:

São algumas limitações importantes que devem ser consideradas, dado que o número de pontos em disputa foi modificado e um ajuste é necessário.

Cabe ressaltar que o Coeficiente de Gini e da análise em si possuem inúmeras limitações que devem ser entendidas como:

  • Não levar o “Fator Tradição” em consideração a um efeito de produtividade dos clubes ao longo do tempo, i.e., um clube grande e antigo carrega uma estrutura financeira/econômica/institucional maior do que os clubes mais novos, e isso pode ser visto na distribuição dos campeões na era dos pontos corridos (Em economia isso seria similar ao efeito de transição de produtividade de uma base instalada produtiva ao longo dos anos);
  • Diferenças econômicas das regiões em que os clubes têm as suas bases;
  • O Coeficiente de Gini olha somente o resultado final sem levar em considerações fatores conjunturais que podem influenciar estes mesmos resultados como gestão, e momento econômico do time, e outros eventos como Olimpíada e Copa do Mundo;
  • O conceito da natureza da geração dos pontos (o que em economia seria a renda) podem ter dinâmicas muito diferentes dado que o índice não captura se os pontos foram gerados através de 3 empates (1 ponto multiplicado por 3) ou por uma vitória e duas derrotas (3 pontos);
  • Analisar a geração de pontos em si ao longo do tempo pode ser bem complicado e não representaria uma comparação plausível. Por exemplo: O Flamengo campeão brasileiro de 2009 seria na melhor das possibilidades um mero quinto colocado em 2014;
  • O índice em si por tratar somente do resultado final, não mostra a transitividade desses pontos ao longo do campeonato. Explico: Se um time nas rodadas finais do campeonato não tem mais nenhum tipo de chance de ir para uma boa competição(Copa Libertadores), de ser rebaixado para divisões inferiores, ou mesmo já conquistou o titulo pode acontecer desses times entrarem com menos disposição para ganhar os jogos. Incentivos financeiros entre os times também podem ocorrerEste artigo da Investopedia fala um pouco sobre esse efeito;

Algumas outras limitações do Coeficiente de Gini podem ser encontradas no trabalho de Tsai, no Working Paper de Osberg, no site do HSRC ou para quem quiser uma critica mais vocal tem esse ensaio do Craig Wright.

Para calcular o Coeficiente de Gini eu usei o código da Olivia Guest apenas para fins de simplicidade, mas qualquer software pode ser usado uma vez que os dados estão disponíveis no repositório.

Como o código têm muito boilerplate code de coisas que eu já fiz no passado e o meu foco é a analise em si, eu não vou comentar o código inteiro.

Vamos dar uma olhada inicial apenas para ver se os dados foram carregados corretamente.

# Imports, 538 theme in the charts and load data
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
%pylab inline'fivethirtyeight')
df_brasileirao = pd.read_csv('dataset-2003-2018.csv', delimiter=';')
print(f'Number of Records: {df_brasileirao.shape[0]} – Number of Columns: {df_brasileirao.shape[1]}')
# Loading Check
view raw hosted with ❤ by GitHub

Aparentemente tudo OK com os dados, então eu vou realizar o calculo do Coeficiente de Gini usando todas as edições do Brasileirão.

Ranking de desigualdade entre todas as edições do Brasileirão usando o Coeficiente de Gini

Para quem quiser calcular, basta usar essas duas funções:

# Gini function as PyGini package
def gini(arr, eps=1e-8):
Reference: PyGini (I owe you a beer @o_guest)
Calculate the Gini coefficient of a numpy array.
Based on bottom eq on [2]_.
.. [2]_
# All values are treated equally, arrays must be 1d and > 0:
arr = np.abs(arr).flatten() + eps
# Values must be sorted:
arr = np.sort(arr)
# Index per array element:
index = np.arange(1, arr.shape[0]+1)
# Number of array elements:
N = arr.shape[0]
# Gini coefficient:
return(np.sum((2*index N 1)*arr))/(N*np.sum(arr))
def get_gini_df(df):
"""Generate DF with Gini Index
df : Pandas Dataframe
Dataframe with Brasileirão data
gini_df : Pandas Dataframe
Returns a Pandas Dataframe with the year, team and gini index
gini_per_year = []
for year in df['year'].unique():
championship_index = gini(np.array(df[df['year'] == year]['points']))
champion = (df[(df['year'] == year) & (df['position'] == 1)]['team'])
gini_per_year.append((year, champion.values[0], round(championship_index, 4)))
gini_df = pd.DataFrame(gini_per_year)
gini_df.columns = ['year', 'team', 'gini']
# Indexing the date field for graph it smoothly
gini_df.set_index('year', inplace=True)
return gini_df
gini_df = get_gini_df(df_brasileirao)
gini_df.sort_values(by=['gini'], ascending=True)
view raw hosted with ❤ by GitHub

Usando os dados que foram carregados e o código acima, eu realizei o cálculo do Coeficiente de Gini e obtive a seguinte tabela:

Coeficiente de Gini Calculado com os respectivos campeões de cada edição

Considerando o Coeficiente de Gini como principal métrica de ordenação, podemos ver que a edição do Brasileirão de 2017, 2005 (ambas com o Corinthians campeão) e de 2009 (Flamengo campeão) foram as que tiveram mais igualdade dentro da era dos pontos corridos.

Por outro lado, as edições de 2018 (Palmeiras campeão), 2014 (Cruzeiro campeão) e 2012 (Fluminense campeão) foram as mais desiguais no que se refere a distribuição final dos pontos. Um fato curioso é que se pegarmos as 5 temporadas mais desiguais, veremos o Palmeiras e o Fluminense com 2 títulos cada (respectivamente 2018, 2016 e 2012, 2010).

Essas informações apontam que Corinthians e Flamengo tendem a ganhar as temporadas com a menor distribuição de pontos finais, e quando Palmeiras e Fluminense ganham são geralmente temporadas mais desiguais da perspectiva de distribuição dos pontos no final. Aliaá isso seria uma boa hipótese inicial para ser testada com mais dados.

Vamos olhar os dois extremos que são as temporadas de 2017 (mais igual) e 2018 (mais desigual).

Brasileirão 2017
Brasileirão 2018

Olhando a distribuição dos pontos aqui, podemos ver que a diferença de pontos entre os campeões foi de 8 pontos (80-72) e considerando a distancia do campeão com o quinto colocado nos dois campeonatos, se em 2017 temos uma distância de 15 pontos (72-57) em 2018 temos 17 pontos (80-63), ou seja muito parecida.

Entretanto, sabendo que o campeão sempre têm um pouco de margem considerando apenas esses extremos, se considerarmos apenas a distância entre o vice-campeão e o quinto colocado em 2017 temos apenas 6 pontos (63-57) enquanto em 2018 essa distância vai para 9 pontos (72-63).

Realizando o mesmo exercício entre o campeão e o pior time do campeonato em 2017 temos uma diferença de 36 pontos (72-36) enquanto em 2018 chegamos a expressiva marca de 57 pontos de diferença (80-23). Isso mostra que mesmo entre os piores times ao longo desses dois campeonatos temos a expressiva diferença de 13 pontos (36 (Atlético Goianiense/2017 – 23 (Paraná). Vamos ficar atentos a essas informações, pois eu vou voltar aqui posteriormente.

Agora eu vou gerar o gráfico de como ficou esse Coeficiente de Gini ao longo do tempo.

def get_graph_ts(df, column, title, label):
"""Generate graph of a Time Series
df : Pandas Dataframe
Dataframe with Brasileirão data
column : string
Column with the metric to be ploted
title : string
Graph title to be displayed
label : string
Name of the series that will be placed
as legend
plt.ylabel('Gini Index')
plt.plot(df[column], label=label)
'Gini Index in Brasileirão',
'Gini Index',
view raw hosted with ❤ by GitHub
Coeficiente de Gini no Brasileirão ao longo do tempo

Em uma primeira análise podemos perceber algumas curiosidades:

  • Ao que parece temos sim uma tendência não tão clara de crescimento da desigualdade, mas com vales e picos bem distintos e considerando que o formato de 20 times têm apenas 12 temporadas completas isso tem que ser levado com cautela;
  • Algo surpreendente é que os vales costumam acontecer nos anos ímpares e os picos nos anos pares. Isso talvez seja explicado por algum efeito externo, tal como as Olimpíadas e Copa do Mundo que ocorrem em anos pares. E uma hipótese bem fraca, mas ainda sim é uma hipótese. (Nota do Autor: Se algueém tiver uma explicação razoável pode me mandar um comentário que eu coloco aqui e dou o crédito);
  • Ao que parece depois de 2011 houve uma subida mais consistente e manutenção desse aumento da desigualdade com o Coeficiente de Gini voltando para baixo de 0.115 apenas em 2017, ou seja, uma janela de 6 anos acima do patamar de 0.115;
  • E falando em 2017, essa temporada ao que parece foi uma total quebra em relação a essa desigualdade, dado que considerando o segundo, terceiro e quarto colocados terminarem com a diferença de apenas 1 ponto3 times ficaram com 43 pontos, com dois deles sendo rebaixados para a segunda divisão.

Entretanto, uma coisa que eu considerei foi que talvez tenha algum efeito que eu chamo de “Última temporada na Série A” ou “Já estou rebaixado mesmo, não tem mais o que fazer” em relação a essa (des)igualdade, dado que sempre saem/entram 4 times por ano (i.e. 20% dos times são trocados todos os anos). Para remover esse potencial efeito, eu vou considerar uma média móvel considerando os 3 últimos campeonatos. Ou seja, sempre haverá a combinação de a) dois anos desiguais e um ano igual e b) dois anos iguais com um ano desigual.

# Some graphs with rolling average
date_range = [2003, 2004, 2005, 2006, 2007,
2008, 2009, 2010, 2011, 2012,
2013, 2014, 2015, 2016, 2017,
def get_graph_ts_rolling_average(ts, title, window, date_range=date_range):
"""Generate graph of a Time Series with a simple rolling average
ts : Pandas Dataframe column
Dataframe column with a metric to be ploted
title : Pandas Dataframe
Graph title to be displayed
window : int
Rolling back window to be considered in the average
date_range : Array
Array to be used in the ploting. Matplotlib has a
very bad way to deal with that, so I need to use this
workaround to place all years properly
plt.plot(date_range, ts.rolling(window=window, center=False).mean(), label='gini');
plt.ylabel('Gini Index')
'Rolling Average in Gini Index in Brasileirão – Rolling Average window=',
view raw hosted with ❤ by GitHub
Coeficiente de Gini quando consideramos uma janela de 3 edições para computar a média

Realizando essa suavização a grosso modo, podemos ver o efeito desse aumento da desigualdade mais claro, quase que de forma linear de 2009 até 2016 sendo quebrado apenas pela famigerada temporada de 2017, e com a temporada de 2018 sofrendo um efeito claro com esse ajuste.

Anteriormente eu falei sobre essa disparidade dos pontos entre os primeiros colocados, entre os últimos colocados e os campeões em relação aos piores times da liga. Podemos perceber que sempre nestes casos extremos temos o campeão vencendo com uma pequena folga, e com os times subsequentes com algum tipo de disputa em ordens diferentes de magnitude.

Para validar esse ponto de que essa desigualdade está aumentando de uma maneira um pouco mais robusta, eu vou remover os nossos outliers desses campeonatos, ou resumindo: Eu vou tirar o campeão e o pior time da temporada da análise.

(Nota do Autor: Eu sei que existem inúmeras abordagens de como se fazer isso da maneira correta, com recursos quase que infinitos na internet de técnicas (1234 e inclusive com testes estatísticos muito robustos de como se fazer isso da maneira correta, caso seja necessário. Por questões de simplicidade e para reforçar o meu ponto, eu estou removendo dado que dentro dessas situações extremas como eu tenho como resultados a) o campeão com uma leve vantagem no final e b) o pior time da temporada sendo muito pior mesmo, a ideia é ver a igualdade dos outros times no campeonato. Lembrem-se que eu quero ver de maneira geral (des)igualdade da liga e tirar o efeito dos supercampeões e dos vergonhosos sacos de pancada da temporada).)

Dito isto, vou remover esses outliers e calcular novamente o Coeficiente de Gini.

def get_brasileirao_no_outliers(df):
"""Generate a DF removing the champion and the worst team of the championship
df : Pandas Dataframe
Dataframe with Brasileirão data
df_concat : Pandas Dataframe
Returns a Pandas Dataframe without the outliers
df_concat = pd.DataFrame()
for year in df['year'].unique():
pos_min = df[df['year'] == year]['position'].min()
pos_max = df[df['year'] == year]['position'].max()
df_filtered = df[(df['year'] == year) \
& (~df['position'].isin([pos_min, pos_max]))]
df_concat = df_concat.append(df_filtered)
return df_concat
def get_gini(df):
"""Generate a DF with the year and the following Gini Index calculated
df : Pandas Dataframe
Dataframe with Brasileirão data
gini_df : Pandas Dataframe
Returns a Pandas Dataframe with the year, and gini index
gini_per_year = []
for year in df['year'].unique():
championship_index = gini(np.array(df[df['year'] == year]['points']))
gini_per_year.append((year, round(championship_index, 4)))
gini_df = pd.DataFrame(gini_per_year)
gini_df.columns = ['year', 'gini']
# Indexing the date field for graph it smoothly
gini_df.set_index('year', inplace=True)
return gini_df
# Outlier removal
df_brasileirao_no_outliers = get_brasileirao_no_outliers(df_brasileirao)
df_brasileirao_no_outliers_gini = get_gini(df_brasileirao_no_outliers)
df_brasileirao_no_outliers_gini.sort_values(by=['gini'], ascending=True)
Gini Index com o campeão e o último colocado do campeonato removidos

Agora temos algumas mudanças significativas no nosso painel que são:

  • Se antes tínhamos as temporadas de 2017, 2005 e 2009 como as mais iguais, agora tiramos a temporada de 2009 vencida pelo Flamengo, e colocamos a de 2007 vencida pelo São Paulo;
  • Pelo lado das temporadas mais desiguais, no caso tínhamos a ordem de desigualdade pelas temporadas de 2018, 2014, e 2012, agora temos a nova ordem em relação as edições de 2014 (Cruzeiro), 2012 (Fluminense) e 2018 (Palmeiras)

Vamos checar como ficou a tabela das edições mais extremas sem os campeões em 2007 e 2014.

Brasileirão 2007
Brasileirão 2014

Olhando as duas tabelas finais em que removemos o campeão e o pior time, logo de cara podemos reparar que enquanto na temporada de 2014 temos dois vales de diferenças de 7 pontos (4° e 5° e 7° e 8°), no campeonato de 2007 a maior diferença foi de 4 pontos (15° e 16°) e fazendo um pequeno exercício de imaginação podemos dizer que o Paraná que teve 41 pontos em 2007 não teria perigo nenhum de ir para a segunda divisão, se tivesse o mesmo número de pontos no campeonato de 2014. Como eu coloquei anteriormente essa análise de transpor um time ao longo do tempo não é muito valida, porém, mencionei isso apenas para ressaltar a importância da distribuição dos pontos ao invés do número absoluto de pontos no resultado final.

Vamos gerar o gráfico apenas para verificar se a tendência do aumento da desigualdade permanece ou não.

Gini Index ao longo do tempo, removendo o campeão e o pior time da temporada
Gini Index ao longo do tempo, removendo o campeão e o pior time da temporada, considerando a média das 3 últimas edições

Olhando o Coeficiente de Gini removendo o campeão e o pior time, podemos ver que ainda temos a tendência de aumento da desigualdade dentro da liga, seja usando a métrica de forma regular ou aplicando uma média móvel com uma janela de 3 temporadas ao longo do tempo.


Ao longo dos dados apresentados nessa análise eu cheguei aos seguintes fatos:

  • Temos uma tendência crescente na desigualdade em relação ao número de pontos entre os times dentro do campeonato brasileiro;
  • Esse crescimento começa de forma mais substancial em 2010;
  • As temporadas mais recentes (2018 e 2017) são respectivamente as temporadas com a maior e a menor desigualdade;
  • Essa tendência do aumento da desigualdade ocorre mesmo removendo o campeão e o pior time da temporada, dado que o campeão nos casos extremos possui uma pequena vantagem em relação ao segundo colocado, e o pior time termina com uma pontuação virtualmente impossível de reverter em uma ocasião de última rodada;
  • Nas temporadas que há uma desigualdade maior existe vales de pontos entre blocos de times, no caso que vimos esses vales foram de 7 pontos (2 vitorias + 1 empate);
  • Nas temporadas mais iguais esses vales de blocos times provavelmente são menores ou em alguns casos inexistentes, no caso observado havia somente um bloco de 4 pontos (1 vitória + 1 empate) e excluindo este fato não havia nenhuma distância maior do que 1 vitória por toda a tabela entre posições imediatamente superiores.

Conclusão e considerações para o futuro

Se tivermos que responder a nossa pergunta principal que foi “O brasileirão está ficando mais injusto ao longo do tempo?” a resposta correta seria:

“Sim. Com o uso do Coeficiente de Gini como métrica para mensurar se há uma desigualdade estrutural mostrou que existem sim elementos latentes dessa desigualdade”.

O leitor mais atento pode repaar que uma coisa que eu tomei muito cuidado aqui foi para não realizar afirmações relativas à competitividade, afirmações relativas às condições financeiras dos clubes, incremento de premiações ou ausência de incentivos para os piores colocados, etc.; dado que estes aspectos são difíceis de mensurar e há pouquíssimos dados disponíveis de maneira confiável para a análise, mas empiricamente talvez podemos realizar algumas afirmações nessas direções.

Em relação a análises futuras existem muitas hipóteses que podem ser testadas como um potencial problema fundamental de competitividade devido ao fato de muitos clubes que não representam uma elite (i.e. muitos times fracos na liga principal) e inclusive existem pautas relativas ao aumento do número dos clubes rebaixados, hipóteses que elencam fatores importantes como a disparidade financeira como um dos potenciais fatores de existir poucos supertimes, há uma hipótese que ganha tração também que é a respeito das cotas financeiras de direitos de exibição de jogos na televisão que compõe grande parte da receita desses clubes que a contar do ano que vem vai punir pesadamente os times que forem rebaixados.

São muitas hipóteses em discussão e muitos aspectos que poderiam ser a razão para o aumento dessa desigualdade. Há um grande número de aspectos que podem ser estudados e no final desse post tem alguns links para referências em outras ligas.

Todos os dados e a análise completa pode ser encontrada no repositório brasileirao-gini. Lá tem todas as instruções para executar passo a passo a análise.

Referências e links úteis

Inequality in the Premier League – Çınar Baymul

An Analysis Of Parity Levels In Soccer – Harvard Sports

Which Sports League has the Most Parity? – Harvard Sports

Major League Soccer and the Effect of Egalitarianism – Harvard Sports

The Gini Coefficient as a Measure of League Competitiveness and Title Uncertainty – Australia Sports Betting

Mourão, P. R., & Teixeira, J. S. (2015). Gini playing soccer. Applied Economics, 47(49), 5229-5246

How “fair” are European soccer leagues? Gini index applied to points distribution of 5 soccer leagues between 2000 and 2015 – r/soccer

Footballomics: Estimating League Disparity Performance with a Point-Rank Gini Index – Christoforos Nikolaou

O Campeonato Brasileiro está ficando mais injusto?

Reproducibility in FastText

A few days ago I wrote about FastText and one thing that is not clear in docs it’s about how to make the experiments reproducible in a deterministic day.

In default settings of train_supervised() method i’m using the thread parameter with multiprocessing.cpu_count() - 1 as value.

This means that we’re using all the CPUs available for training. As a result, this implies a shorter training time if we’re using multicore servers or machines.

However, this implicates in a totally non-deterministic result because of the optimization algorithm used by fastText (asynchronous stochastic gradient descent, or Hogwildpaper here), the obtained vectors will be different, even if initialized identically.

This very gentle guide of FastText with Gensim states that:

for a fully deterministically-reproducible run, you must also limit the model to a single worker thread (workers=1), to eliminate ordering jitter from OS thread scheduling. (In Python 3, reproducibility between interpreter launches also requires use of the PYTHONHASHSEED environment variable to control hash randomization).

Radim Řehůřek in FastText Model

So for that particular reason the main assumption here it’s even playing in a very stocastic environment of experimentation we’ll consider only the impact of data volume itself and abstract this information from the results, for the reason that this stocastic issue can play for both experiments.

To make reproducible experiments the only thing that it’s needed it’s to change the value of thread parameter from multiprocessing.cpu_count() - 1to 1.

So for the sake of reproducibility the training time will take longer (in my experiments I’m facing an increase of 8000% in the training time.

Reproducibility in FastText

Dica de Python: Dask

Para quem não aguenta mais sofrer com o Pandas e não quer lidar com as inúmeras limitações do Scala o Dask é uma ótima biblioteca para manipulação de dados e computação em Python.

Direto da documentação:

Familiar: Provides parallelized NumPy array and Pandas DataFrame objects

Flexible: Provides a task scheduling interface for more custom workloads and integration with other projects.

Native: Enables distributed computing in pure Python with access to the PyData stack.

Fast: Operates with low overhead, low latency, and minimal serialization necessary for fast numerical algorithms

Scales up: Runs resiliently on clusters with 1000s of cores

Scales down: Trivial to set up and run on a laptop in a single process

Responsive: Designed with interactive computing in mind, it provides rapid feedback and diagnostics to aid humans

Dica de Python: Dask

FastText – A great tool for Text Classification

In some point of time, I’ll post a field report about FastText in a project for Text Classification. My opinion until this moment (16.03.19): For a fast Alpha version of a text classification with robust use of Bag-of-Tricks and WordNGrams it’s amazing in terms of practical results (especially Recall) and speed of development.


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

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

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

Principais soluções do

Agora que já sabemos um pouco sobre essa solução, vamos entender um pouco do ecossistema de soluções do H2O, e ver as principais características e aplicações de cada uma.


Essa plataforma é o carro chefe da empresa, o qual eles apostam tanto na versão Desktop para aplicação de Machine Learning quando também na versão para processamento distribuído para altos volumes de dados.

Essa versão tem alguns algoritmos prontos on the shelf como boosting, regressão linear e logística, algoritmos baseados em árvores, e alguns algoritmos que utilizam gradiente como método de otimização. Nada muito complexo, mas bem funcional.

Essa versão é ideal para usar se você quer conhecer um pouco mais da ferramenta e não quer gastar muito tempo instalando ou configurando coisas antes de sair aplicando os algoritmos, ou mesmo para um teste inicial das funções de processamento distribuído em cluster. 

Abaixo, um pouco da arquitetura da solução:


Sparkling Water

Essa solução tem como principal característica utilizar os próprios algoritmos, mas com a vantagem de usar todas as features de processamento distribuído e integração do Spark. Nesta solução todas as tarefas de computação também podem ser feitas dentro do Spark usando Scala e com uma interface via Web.

Essa solução é a mais recomendada para construção de aplicações de Machine Learning seja para microserviços ou até mesmo para embutir dentro de uma plataforma/sistema toda a parte algorítmica e computacional.

h2ospark h2ospark2

Deep Water

O Deep Water é a solução voltada para implementação de Deep Learning usando otimização computacional com GPUs com frameworks como o Tensor Flow, Theano, Caffe entre outros.

Neste caso, a plataforma do H2O será a interface onde serão incorporados todos os parâmetros de treinamento do modelo (cross validation, amostragem, critério de parada, hiperparametrização, etc) e o backend com o Tensor Flow, Theano etc. faz o processamento utilizando GPUs.


O Steam é uma plataforma que realiza todo o link entre os modelos de machine learning usando o H2O e também propriedades de desenvolvimento para incorporar modelos de Machine Learning em aplicações, tudo isso de forma colaborativa, algo muito similar ao Domino.

A principal vantagem do Steam é que ele abstraí toda a parte de engenharia por trás da tarefa de incorporar modelos e machine learning em produção como infraestrutura, auto-scale de infra estrutura de acordo com a carga de requisições, bem com algumas tarefas de Data Science como retreinamento de modelos; além de reduzir e muito os custos/investimentos de TI.


Agora que sabemos quais são os principais produtos do H2O, em breve teremos alguns posts com alguns tutoriais explorando um pouco mais essa ferramenta.

Links úteis

Documentação Técnica

Principais soluções do

Porque você REALMENTE deveria considerar o como uma das ferramentas do seu stack de Machine Learning?

Quando falamos de ferramentas de machine learning logo vem a cabeça a tríade Tensor Flow,  Scikit-Learn e Spark MLLib.
Contudo, há uma ferramenta que vem discretamente ganhando espaço que é o
Essa ferramenta nasceu originalmente em 2011 em que o seu time teve como principal objetivo democratizar e tornar escalável machine learning através de uma plataforma mais visual e que tivesse uma boa experiência para os usuários, independente do seu nível técnico.
Algumas características do
Abaixo alguns vídeos sobre o em ação: para detecção de fraudes
Customer Churn usando

Em alguns posts futuros falaremos um pouco sobre questões de arquitetura, e mergulharemos em alguns tutoriais.
Porque você REALMENTE deveria considerar o como uma das ferramentas do seu stack de Machine Learning?

Comparação entre R e Python utilizando Florestas Aleatórias e Classificação

Neste post do blog do Yhat tem o código, os dados e os resultados.

Pessoalmente gosto muito da abordagem dos autores em comparação de classificadores usando as métricas de Acurácia, Erro Quadrático Médio e  para regressão e tempo de treinamento.

Para projetos curtos de avaliação de uma série de classificadores essas medidas são suficientes para dar uma linha de base. Essas medidas podem auxiliar na escolha de quais modelos estão com melhor convergência e podem indicar um melhor tratamento dos dados em termos de quais variáveis são pertinentes ao modelo escolhido.



Comparação entre R e Python utilizando Florestas Aleatórias e Classificação