Mãos nas estatísticas Bayesianas com Python, PyMC3 e ArviZ

Inferências Gaussianas, Verificações Preditivas Posiores, Comparação de Grupos, Regressão Linear Hierárquica

Susan Li Seg. 17 de jul · 10 min ler

Se você acha que o teorema de Bayes é contra-intuitivo e as estatísticas bayesianas , que se baseiam no teorema de Baye, podem ser muito difíceis de entender. Estou com você.

Existem inúmeras razões pelas quais devemos aprender as estatísticas Bayesianas , em particular, as estatísticas Bayesianas estão surgindo como uma estrutura poderosa para expressar e entender as redes neurais profundas de próxima geração.

Acredito que para as coisas que temos que aprender antes de podermos fazer, aprendemos fazendo-as . E nada na vida é tão difícil que não possamos facilitar a maneira como a tomamos .

Então, esta é a minha maneira de tornar isso mais fácil: ao invés de muitas teorias ou terminologias no começo, vamos nos concentrar na mecânica da análise Bayesiana, em particular, como fazer análise Bayesiana e visualização com PyMC3 e ArviZ . Antes de memorizar as infinitas terminologias, codificaremos as soluções e visualizaremos os resultados, usando as terminologias e teorias para explicar os modelos ao longo do caminho.

PyMC3 é uma biblioteca Python para programação probabilística com uma sintaxe muito simples e intuitiva. ArviZ , uma biblioteca Python que trabalha de mãos dadas com o PyMC3 e pode nos ajudar a interpretar e visualizar distribuições posteriores.

E aplicaremos métodos bayesianos a um problema prático, para mostrar uma análise Bayesiana ponta-a-ponta, que passará do enquadramento da questão para a construção de modelos, para a obtenção de probabilidades prévias de implementação em Python da distribuição posterior posterior.

Antes de começarmos, vamos tirar algumas intuições básicas do caminho:

Os modelos bayesianos também são conhecidos como modelos probabilísticos porque são construídos usando probabilidades. E Bayesian usa as probabilidades como uma ferramenta para quantificar a incerteza. Portanto, as respostas que recebemos são distribuições e não estimativas pontuais.

Etapas da Abordagem Bayesiana

Etapa 1: Estabelecer uma crença sobre os dados, incluindo as funções Prior e Likelihood.

Passo 2, Use os dados e a probabilidade, de acordo com a nossa crença nos dados, para atualizar o nosso modelo, verifique se o nosso modelo está de acordo com os dados originais.

Etapa 3, Atualize nossa visão dos dados com base em nosso modelo.

Os dados

Como estou interessado em usar o aprendizado de máquina para otimização de preço, decido aplicar os métodos bayesianos a um conjunto de dados de preços de tickets do trem de alta velocidade espanhol que pode ser encontrado aqui . Aprecie a equipe de Gurus por raspagem do conjunto de dados.

 das estatísticas de importação do scipy 
importar arviz como az
import numpy como np
import matplotlib.pyplot como plt
import pymc3 como pm
importar no mar como sns
importar pandas como pd
a partir de theano import shared
from sklearn import preprocessing
print (formato 'Executando no PyMC3 v {}'. (pm .__ version__)) data = pd.read_csv ('renfe.csv')
data.drop ('Sem nome: 0', eixo = 1, local = Verdadeiro)
data = data.sample (frac = 0.01, random_state = 99)
data.head (3)

tabela 1

 data.isnull (). sum () / len (dados) 

figura 1

Existem 12% dos valores na coluna de preços em falta, resolvo preenchê-los com a média dos respectivos tipos de tarifas. Preencha também as outras duas colunas categóricas com os valores mais comuns.

 dados ['train_class'] = dados ['train_class']. fillna (dados ['train_class']. mode (). iloc [0]) 
data ['tarifa'] = data ['tarifa']. fillna (data ['tarifa']. mode (). iloc [0])
data ['price'] = data.groupby ('fare'). transform (lambda x: x.fillna (x.mean ()))

Inferências gaussianas

 az.plot_kde (data ['price']. valores, rug = True) 
plt.yticks ([0], alfa = 0);

Figura 2

O gráfico do KDE do preço do bilhete ferroviário mostra uma distribuição semelhante a Gauss , exceto por várias dezenas de pontos de dados que estão longe da média.

Vamos supor que uma distribuição gaussiana seja uma descrição adequada do preço do bilhete ferroviário. Como não sabemos a média ou o desvio padrão, devemos definir prioridades para os dois. Portanto, um modelo razoável poderia ser o seguinte.

Modelo

Vamos realizar inferências gaussianas sobre os dados do preço do bilhete. Aqui estão algumas das opções de modelagem que entram nisso.

Vamos instanciar os modelos no PyMC3 assim:

  • As especificações do modelo no PyMC3 são agrupadas em uma declaração com.

Escolhas de priors:

  • ?, média de uma população. Distribuição normal, muito larga. Eu não sei os possíveis valores de ?, eu posso definir os antecedentes refletindo minha ignorância. Pela experiência eu sei que o preço do bilhete de trem não pode ser menor que 0 ou maior que 300, então eu defino os limites da distribuição uniforme como 0 e 300. Você pode ter experiência diferente e definir os diferentes limites. Isso é totalmente bem. E se você tiver informações prévias mais confiáveis do que eu, por favor use!
  • ?, desvio padrão de uma população. Só pode ser positivo, portanto use a distribuição HalfNormal. Mais uma vez, muito largo.

Opções para a função de verossimilhança do preço do bilhete:

  • y é uma variável observada representando os dados provenientes de uma distribuição normal com os parâmetros ? e ?.
  • Desenhe 1000 amostras posteriores utilizando a amostragem NUTS.

Usando o PyMC3, podemos escrever o modelo da seguinte forma:

model_g.py

O y especifica a probabilidade. Esta é a maneira em que dizemos ao PyMC3 que queremos condicionar o desconhecido no know (data).

Nós traçamos o traçado do modelo gaussiano. Isso é executado em um gráfico Theano sob o capô.

 az.plot_trace (trace_g); 

Figura 3

  • À esquerda, temos um gráfico do KDE, – para cada valor de parâmetro no eixo x, obtemos uma probabilidade no eixo y que nos diz qual é a probabilidade desse valor de parâmetro.
  • À direita, obtemos os valores amostrados individuais em cada etapa durante a amostragem. A partir do gráfico de rastreamento, podemos obter visualmente os valores plausíveis do posterior.
  • O gráfico acima tem uma linha para cada parâmetro. Para este modelo, o posterior é bidimensional, e assim a figura acima está mostrando as distribuições marginais de cada parâmetro.

Há algumas coisas para notar aqui:

  • Nossas cadeias de amostragem para os parâmetros individuais (à esquerda) parecem bem convergentes e estacionárias (não há grandes variações ou outros padrões estranhos).
  • A estimativa posterior máxima de cada variável (o pico nas distribuições do lado esquerdo) é muito próxima dos parâmetros verdadeiros.

Podemos plotar uma distribuição conjunta de parâmetros.

 az.plot_joint (trace_g, tipo = 'kde', fill_last = False); 

Figura 4

Não vejo nenhuma correlação entre esses dois parâmetros. Isso significa que provavelmente não temos colinearidade no modelo. Isso é bom.

Também podemos ter um resumo detalhado da distribuição posterior de cada parâmetro.

 az.summary (trace_g) 

mesa 2

Também podemos ver visualmente o resumo acima, gerando um gráfico com a média e a Densidade Posterior Mais Alta (HPD) de uma distribuição, e para interpretar e relatar os resultados de uma inferência Bayesiana.

 az.plot_posterior (trace_g); 

Figura 5

  • Ao contrário da inferência freqüentista , na inferência bayesiana , obtemos toda a distribuição dos valores.
  • Toda vez que o ArviZ calcula e relata um HPD, ele usará, por padrão, um valor de 94%.
  • Por favor, note que os intervalos de HPD não são os mesmos que os intervalos de confiança.
  • Aqui podemos interpretar como tal que há 94% de probabilidade de que a crença esteja entre 63,8 e 64,4 euros para o preço médio do bilhete.

Podemos verificar a convergência das cadeias formalmente usando o teste Gelman Rubin. Valores próximos a 1,0 de convergência média.

 pm.gelman_rubin (trace_g) 

 bfmi = pm.bfmi (trace_g) 
max_gr = max (np.max (gr_stats) para gr_stats em pm.gelman_rubin (trace_g) .values ())
(pm.energyplot (trace_g, legenda = Falso, figsize = (6, 4)). set_title (formato "BFMI = {} nGelman-Rubin = {}". (bfmi, max_gr)));

Figura 6

Nosso modelo convergiu bem e a estatística Gelman-Rubin parece bem.

Verificações preditivas posteriores

  • As verificações preditivas posteriores (PPCs) são uma ótima maneira de validar um modelo. A ideia é gerar dados a partir do modelo usando parâmetros de desenhos do posterior.
  • Agora que calculamos o posterior, vamos ilustrar como usar os resultados da simulação para obter previsões.
  • A seguinte função irá sortear aleatoriamente 1000 amostras de parâmetros do rastreamento. Então, para cada amostra, ele irá desenhar 25798 números aleatórios de uma distribuição normal especificada pelos valores de ? e ? nessa amostra.
 ppc = pm.sample_posterior_predictive (trace_g, samples = 1000, modelo = model_g) np.asarray (ppc ['y']). 

Agora, ppc contém 1000 conjuntos de dados gerados (contendo 25798 amostras cada), cada um usando um parâmetro diferente da posterior.

 _, ax = plt.subplots (figsize = (10, 5)) 
ax.hist ([y.mean () para y em ppc ['y']], bins = 19, alfa = 0,5)
ax.axvline (data.price.mean ())
ax.set (título = 'Preversa posterior da média', xlabel = 'média (x)', ilabel = 'Frequência');

Figura 7

A média inferida é muito próxima da média real do preço do bilhete ferroviário.

Comparação de grupos

Podemos estar interessados em comparar o preço em diferentes tipos de tarifas. Vamos nos concentrar em estimar o tamanho do efeito, ou seja, quantificar a diferença entre duas categorias de tarifas. Para comparar as categorias de tarifa, vamos usar a média de cada tipo de tarifa. Como somos Bayesianos, trabalharemos para obter uma distribuição posterior das diferenças de médias entre as categorias de tarifa.

Nós criamos três variáveis:

  • A variável de preço, representando o preço do bilhete.
  • A variável idx, uma variável dummy categórica para codificar as categorias de tarifa com números.
  • E finalmente a variável grupos, com o número de categorias tarifárias (6)
 price = data ['price']. values 
idx = pd.Categorical (data ['tarifa'],
categories = ['Flexible', 'Promo', 'Promo +', 'Adulto ida', 'Mesa', 'Individual-Flexível']).
grupos = len (np.unique (idx))

O modelo para o problema de comparação de grupos é quase o mesmo que o modelo anterior. a única diferença é que ? e ? serão vetores ao invés de variáveis escalares. Isso significa que, para as prioridades, passamos um argumento de forma e, para a probabilidade, indexamos apropriadamente as médias e as variáveis sd usando a variável idx:

comparando_groups.py

Com 6 grupos (categorias de tarifas), é um pouco difícil traçar um gráfico de traços para ? e ? para cada grupo. Então, criamos uma tabela de resumo:

 flat_fares = az.from_pymc3 (trace = trace_groups) 
fares_gaussian = az.summary (flat_fares)
fares_gaussian

Tabela 3

É óbvio que existem diferenças significativas entre os grupos (ou seja, categorias de tarifa) na média.

Para torná-lo mais claro, plotamos a diferença entre cada categoria de tarifa sem repetir a comparação.

  • O d de Cohen é um tamanho de efeito apropriado para a comparação entre duas médias. Cohen's d introduz a variabilidade de cada grupo usando seus desvios padrão.
  • A probabilidade de superioridade (ps) é definida como a probabilidade de que um ponto de dados obtido aleatoriamente de um grupo tenha um valor maior do que um obtido aleatoriamente de outro grupo.

difference_plot.py Figura 8

Basicamente, o gráfico acima nos diz que nenhum dos casos de comparação acima, onde o 94% de HPD inclui o valor de referência de zero. Isso significa que para todos os exemplos, podemos descartar uma diferença de zero. As diferenças médias entre 6,1 e 63,5 euros são suficientemente grandes para justificar que os clientes comprem bilhetes de acordo com diferentes categorias de tarifas.

Regressão Linear Hierárquica Bayesiana

Queremos construir um modelo para estimar o preço do bilhete ferroviário de cada tipo de trem e, ao mesmo tempo, estimar o preço de todos os tipos de trem. Esse tipo de modelo é conhecido como modelo hierárquico ou modelo multinível.

  • Codificando a variável categórica.
  • A variável idx, uma variável dummy categórica para codificar os tipos de trem com números.
  • E finalmente a variável grupos, com o número de tipos de trem (16)

encode_cat.py Tabela 4

A parte relevante dos dados que modelaremos será semelhante à acima. E estamos interessados em saber se diferentes tipos de trem afetam o preço do ingresso.

Modelo Hierárquico

hierarchical_model.py Figura 9

Os marginais posteriores na coluna da esquerda são altamente informativos, “?_?_tmp” nos informa os níveis médios de preço do grupo, “?_?” nos diz que a compra da categoria de tarifa “Promo +” aumenta o preço significativamente comparado ao tipo de tarifa “Adulto ida” a categoria "Promo" aumenta o preço de forma significativa em comparação ao tipo de tarifa "Promo +" e assim por diante (sem massa abaixo de zero).

 pm.traceplot (hierarchical_trace, var_names = ['?_tmp'], coords = {'?_tmp_dim_0': intervalo (5)}); 

Figura 10

Entre os 16 tipos de trem, podemos observar como os cinco tipos de trem se comparam em termos de preço do ingresso. Podemos ver, olhando para os marginais, “?_tmp”, que há alguma diferença nos preços entre os tipos de trem; as diferentes larguras estão relacionadas com a confiança que temos em cada estimativa de parâmetro – quanto mais medições por tipo de trem, maior será a nossa confiança.

A quantificação da incerteza de algumas de nossas estimativas é uma das coisas poderosas da modelagem bayesiana. Temos um intervalo credível bayesiano para o preço de diferentes tipos de trem.

 az.plot_forest (hierarchical_trace, var_names = ['?_tmp', '?'], combinado = Verdadeiro); 

Figura 11

Por último, podemos querer calcular r ao quadrado:

 ppc = pm.sample_posterior_predictive (hierarchical_trace, samples = 2000, model = hierarchical_model) 
az.r2_score (data.price.values, ppc ['fare_like'])

O objetivo deste post é aprender, praticar e explicar bayesiano, não para produzir os melhores resultados possíveis a partir do conjunto de dados. Caso contrário, teríamos ido diretamente com o XGBoost .

O notebook Jupyter pode ser encontrado no Github , aproveite o resto da semana.

Referências:

GLM: Regressão Linear Hierárquica – Documentação do PyMC3 3.6

O conjunto de dados de radônio de Gelman et al. (2007) é um clássico para modelagem hierárquica. Neste conjunto de dados, a quantidade de…

docs.pymc.io

Texto original em inglês.