Regressão Bayesiana com PyMC3
Implementando a regressão Bayesiana em Python de forma simples e prática.
Índice
∘ Introdução ao PyMC3
∘ Monte Carlo
∘ Cadeias de Markov
∘ Markov Chain Monte Carlo
∘ Implementando o PyMC3
∘ Estudo de caso: Estimando a elasticidade-preço na demanda (EDP)
∘ Conclusão
Olá, olá, olá! Seja mais uma vez bem-vindo a mais um Turing Talks!
Hoje vamos introduzir a vocês, leitores, algumas práticas de análise e visualização de dados por meio da Regressão Bayesiana através do pacote computacional PyMC3. Nas ultimas décadas a Estatística Bayesiana teve um grande salto de popularidade. Tal ascensão se deve, em grande parte, a avanços científicos nos métodos de Monte Carlo via cadeias de Markov (MCMC) combinado ao barateamento dos recursos computacionais.
Introdução ao PyMC3
PyMC3 é uma biblioteca Python de programação probabilística com uma sintaxe muito simples e intuitiva para interpretar e visualizar distribuições a posteriori. O PyMC3 utiliza métodos de Monte Carlo via cadeias de Markov (MCMC). Uma grande vantagem do MCMC é sua flexibilidade — já que funciona com qualquer modelo e qualquer distribuição a priori. A seguir, veremos o essencial para entender os dois blocos de construção de um modelo de regressão bayesiano PyMC3: Cadeias de Monte Carlo e Markov.
Monte Carlo
De uma forma não rigorosa, podemos definir Monte Carlo como uma maneira de aproximar uma quantidade gerando-se números aleatórios. Por exemplo, de uma forma intuitiva, dado um circulo de raio r = 5:
A partir da formula A∼πr² sabemos que a área do circulo é A∼78.5. Mas como podemos aproximar a área sem a fórmula?
- Para começar, podemos definir um quadrado 10 x 10 sobre o círculo.
- Em seguida, geramos pontos aleatórios dentro do quadrado. Quanto mais pontos, mais precisa a aproximação.
- Considerando apenas 25 pontos, onde 19/25 pontos = 76% caíram dentro do círculo, isto é, significa que a área do círculo é aproximadamente 76% da área do quadrado, ou : 76×10×10/100 ≈ 76. Nada mal usando números aleatórios!
Cadeias de Markov
As Cadeias de Markov, (explicadas com mais detalhes no Turing Talks sobre POS Tagging), são modelos de uma sequência de estados onde cada transição entre esses estados ocorre com uma certa probabilidade. Para exemplificar imagine um urso que só caça, come e dorme. A tabela a seguir mostra as probabilidades de transição entre esses três estados.
Da primeira linha, vemos que, se o urso estiver caçando agora, a probabilidade dele caçar novamente ou dormir num estado futuro é 10% e a probabilidade de comer é 80%.
Algumas Cadeias de Markov têm a propriedade de que, após passarem de estados muitas vezes, elas atingirão o chamado estado estacionário. Isso significa que não importa onde o urso tenha começado, as probabilidades de ele estar em determinados estados em um futuro distante são as mesmas (vide a tabela abaixo).
Markov Chain Monte Carlo
Agora juntaremos ambos os conceitos, para extrair um conjunto de diferentes valores plausíveis para os parâmetros (ou seja, os coeficientes) de um modelo probabilístico. Começamos gerando um ponto aleatório. Depois outro, próximo ao primeiro — essa é a parte de Monte Carlo, geração aleatória. Então, verificamos o quão bem esse novo ponto explica nossos dados, ou: qual é a verossimilhança com esse valor do parâmetro. Então, aceitamos ou rejeitamos esse novo ponto. Quanto melhor explicar os dados, maior a probabilidade de aceitá-los.
Aqui, os pontos que explicam bem nossos dados — e, portanto, são aceitos— , estão em cor verde. Eventualmente, há muitos pontos aceitos. Isso gera uma Cadeia de Markov, e as probabilidades de amostragem de valores específicos convergem para o estado estacionário, que é nossa distribuição a posteriori. Finalmente, descartamos as observações iniciais, que correspondem a amostragem que antecede a convergência da Cadeia de Markov, pois são simplesmente aleatórios. Os restantes (pontos verdes mais escuros) são as melhores estimativas dos valores de parâmetros da distribuição a posteriori.
Sem mais delongas, agora vamos colocar todos esses conceitos em prática, solucionando (ao estilo Bayesiano) um problema de negócio com dados reais através do pacote PyMC3. Uma vez interessados em otimização de preços, decidimos aplicar esses métodos Bayesianos a um conjunto de dados históricos sobre preços de abacate e volume de vendas em mercados dos EUA (O dataset pode ser encontrado aqui )
Implementando o PyMC3
Bibliotecas
Primeiro temos que importar as bibliotecas necessárias:
# Importando o PyMC3 e também o arviz para visualização
import pymc3 as pm
import arviz as az# Importando outros pacotes importantes de data science
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
Estudo de caso: Estimando a elasticidade-preço na demanda (EDP)
Após anos de aprendizado no grupo “Turing USP” você acabou de conseguir o emprego dos sonhos: um cargo de analista de dados na “Incorporação Abacate”, uma empresa que vende, sim, você adivinhou certo, abacates.
A empresa preparou alguns dados para você, incluindo o preço do abacate por unidade da fruta, o volume vendas para um determinado preço e um indicador para classificar se a fruta é orgânica. A partir desse conjunto de dados você obteve um dataframe.
Objetivo: Estimar a elasticidade-preço na demanda (EDP) de abacates e otimizar os preços.
(elasticidade-preço = impacto da mudança de preço no volume de vendas)
Etapas da Abordagem Bayesiana:
1 — Fitar um modelo de regressão Bayesiana.
2 — Inspecionar a confiabilidade do modelo.
3 — Predizer o volume de vendas para diferentes preços e propor o preço que maximiza o lucro e mensurar a incerteza associada.
Etapa 1
Fitar um modelo de regressão Bayesiana.
Como primeiro passo é necessário definir o modelo de regressão linear que consiste na variável resposta, que é o volume de vendas e o parâmetros da regressão:
Onde:
N ➝ denota uma distribuição normal. Vale ressaltar a diferença entre a abordagem Bayesiana em relação a abordagem frequentista: aqui tratamos a resposta como uma variável aleatória e assumimos que ela tem uma distribuição normal com a média definida pela equação de regressão e algum desvio padrão (σ).
β₀ ➝ Corresponde ao volume de vendas sem o impacto de qualquer feature;
β₁ ➝ A EDP, ou seja, o impacto da feature “price” no volume de vendas;
β₂ ➝ O impacto da feature “TypeOrganic” no volume de vendas.
Além disso, você assume que a EDP seja a mesma tanto para os abacates orgânicos e não orgânicos. E você também espera que a EDP seja negativa, pois quanto mais alto for o preço, menores serão as vendas (esse é o caso da maioria dos bens de consumo). Para incorporar este conhecimento a priori no modelo, uma vez familiarizado com o dataset você decide utilizar uma distribuição normal com média -80 representando o termo a priori para o preço.
Como construir um modelo deste tipo no PyMC3?
De uma forma bem simples o modelo é construído usando o dataframe e dois módulos do PyMC3:
formula = "volume ~ price + type_organic"
with pm.Model() as model:
priori = {"price": pm.Normal.dist(mu = -80)}
pm.GLM.from_formula(formula, data=avocado, priors=priori)
# Fitando model utilizando 1000 draws
traço = pm.sample(draws = 1000, tune = 500, chains = 4, init="adapt_diag")Explicando o código acima, definiu-se a instancia pm.Model() como model e nas instruções determinou-se:
- O modelo a partir da função
GLM.from_formula, passando como argumentos da função a formula "volume ~ price + type_organic", o dataframe e a distribuição a priori (vale ressaltar que, se fosse o caso, poderíamos definir um conhecimento a priori para todas as features do modelo, assim como poderíamos definir também uma função de verossimilhança; nesse caso, deixamos o "default"). - Por fim, com a função
sample, especificamos os draws (observações) que correspondem à geração de números aleatórios de onde será extraído a verdadeira distribuição a posteriori. No processo de amostragem, definimos 1000 draws válidos que correspondem ao estado estacionário e 500 draws iniciais especificados pelo parâmetro tune para fins de convergência. Estes serão descartados. No parâmetro chains (cadeias), especificamos o número de repetições independentes do processo de amostragem para o caso de alguma cadeia não convergir. Portanto, no total tem-se 4000 draws para cada um dos 4 parâmetros do modelo.
Etapa 2
Inspecionar a confiabilidade do modelo.
Prezando-se por boas práticas de análise de dados, agora os draws podem ser inspecionados por meio das funções plot_trace e summary.
# Gráficando as observações
pm.plot_trace(traço)# Printando um resumo estatístico
resumo = pm.summary(traço)
print(resumo)
Observaremos os resultados acima, começando pela função summary, pela qual podemos obter uma tabela com um resumo estatístico dos dados. Nas duas primeiras colunas, podemos ver a média e o desvio padrão dos draws para cada parâmetro. Em seguida, temos as extremidades dos intervalo de credibilidade, e vale destacar a última coluna chamada r_hat. Esse número só é calculado se executarmos mais de uma cadeia. Valores de r_hat maiores que 1 são um indicativo que cadeias de Markov não convergiram. Em relação à função plot_trace, esta produz dois subplots por parâmetro:
- À esquerda, temos o gráfico da função de densidade a posteriori. Há quatro linhas por parâmetro em cada subplot, uma para cada cadeia de Markov. O fato de serem todos semelhantes é um indicativo de confiabilidade do modelo. Apesar de ocorrer uma pequena instabilidade em uma das cadeias da feature
pricea densidade ainda está bem próxima das outras. - À direita, temos os gráficos dos 4000 draws por cadeia. Percebe-se que estes oscilam com amplitude em torno de alguma média constante, o que prova uma boa convergência. Assim, em suma, não precisamos nos preocupar com isso e podemos usar o modelo com segurança para otimizar o preço!
Etapa 3
Otimizar o preço e mensurar a incerteza associada.
Uma vez efetuado a inspeção do modelo, vamos ao que interessa: Seu chefe pediu que você forneça o preço do abacate que maximiza o lucro e qual lucro pode ser esperado. Além disso, ele quer que o preço seja divisível por US$ 0.25 para que os clientes possam pagar facilmente com quarto de dólar.
Neste caso você deve utilizar o modelo desenvolvido para prever o volume e o lucro considerando alguns preços razoáveis. Em seguida, você visualizará as distribuições preditivas para escolher o preço ideal. Por fim, você ira calcular o intervalo de credibilidade de previsão de lucro.
Para isso é necessário calcular a média da distribuição a posteriori dos 4 parâmetros obtidos anteriormente. Assim sendo:
# Obtendo a média a posteriori de cada parâmetro do modelo
intercepto_media = np.mean(traço.get_values("Intercept"))
price_media = np.mean(traço.get_values("price"))
organic_media = np.mean(traço.get_values("type_organic"))
sd_media = np.mean(traço.get_values("sd"))Para cada preço provável que maximiza o lucro, calcula-se a média preditiva para obter a amostragem da distribuição do volume de vendas e, assim, usar o volume predito para prever o lucro. Veja:
lucro_pred_por_preço = {}
for preço in [0.5, 0.75, 1, 1.25]:
pred_media = (intercepto_media + price_media*preço + organic_media)
volume_pred = np.random.normal(pred_mean, sd_mean, size=1000)
lucro_pred = preço * volume_pred
lucro_pred_por_preço.update({preço: lucro_pred})Por fim, plotamos um gráfico em floresta do lucro predito:
pm.plot_forest(lucro_pred_por_preço)
plt.show()
Vale ressaltar que optar por esse tipo de gráfico em floresta é uma escolha conveniente que, nesse caso, supera os gráficos de densidade usando seaborn, em que tudo ficaria muito confuso. Para implementar o gráfico em floresta, conforme o código acima, só foi necessário criar um dicionário que denota o lucro predito conforme o preço estipulado. Conforme os resultados acima, cada densidade é representada com uma linha azul no gráfico. A linha mais grossa no meio de cada densidade denota o intervalo interquartil (grau de espalhamento de dados), enquanto que as linhas mais finas mostram o intervalo de credibilidade de 94% (que é o default do PyMC3).
Baseado nesses resultados, a escolha que maximiza o lucro é US$ 0.75. Agora, conhecendo o preço ideal, pode-se melhorar o cálculo de intervalo de credibilidade para 99%, onde a chance de previsão para o caso de lucro (valores positivos) ou prejuízo (valores negativos) não cair nesse intervalo é de 1%, portanto sendo muito improvável um erro de estimativa fora desse intervalo.
# função hpd para obter um novo intervalo de credibilidade
opt_hpd = pm.hpd(lucro_pred_por_preço[0.75],hdi_prob = 0.99)
print("Região de credibilidade (99%):", opt_hpd)
Good joob! Com um preço de abacate mais alto (US$ 1.25) ou mais baixo (US$ 0.5), sua empresa perderia lucro, mas graças às suas habilidades com Estatística Bayesiana, eles conseguiram definir o melhor preço possível. Mais do que isso, sabendo da incerteza na previsão do lucro, eles podem se preparar até para o pior cenário (no qual o lucro é negativo!).
Código
Para ver na integra o código utilizado para fazer este post, confira o notebook no GitHub, aqui:
Conclusão
Neste artigo, colocamos em prática fundamentos da Estatística Bayesiana adotando uma abordagem diferente da Regressão Linear em comparação as Estatísticas Frequentistas, no caso, utilizando Monte Carlo e Cadeias de Markov. Em seguida, aplicamos esses conceitos solucionando um estudo de caso através de um modelo linear bayesiano usando o pacote PyMC3.
Vimos bastante sobre Estatística Bayesiana neste texto. Ressalta-se, porém, que nós mostramos apenas a ponta do iceberg. Por exemplo: dentro da estrutura bayesiana, podemos combinar várias equações de regressão (que chamamos de modelos hierárquicos); existem também versões bayesianas de praticamente todos os modelos estatísticos comuns, como a regressão logística ou de Poisson; além disso, há diversas técnicas populares de machine learning que dependem da inferência bayesiana— inclusive, temos as redes neurais bayesianas!
Por isso, para mais informações, encorajo você, leitor, a verificar a documentação do PyMC3 e o livro Think Bayes de Allen Downey, que é uma excelente leitura introdutória — e está disponível gratuitamente no site do autor!
Além disso, Se quiser conhecer um pouco mais sobre o que fazemos no Turing USP, não deixe de nos seguir nas redes sociais: Facebook, Instagram, LinkedIn e, claro, ler nossos posts no Medium. E se preferir acompanhar ainda mais de perto e participar de nossas discussões e eventos, você também pode entrar no nosso servidor do Discord.
Até a próxima!
