• Português
  • English
  • Español
  • Introdução a Programação Probabilística com PyMC3

    Por: em 1 de dezembro de 2017

    Quando lidamos com análise de dados, geralmente lidamos com dados com alguma incerteza. Essas incertezas podem ser desde incertezas intrínsecas aos dados, assim como erros de medida ou de armazenamento dos dados. Para modelar o nosso problema, portanto, é natural que usemos técnicas probabilísticas de inferência. Essas técnicas podem ser simplificações de modelos probabilísticos (como boa parte das técnicas de machine learning), mas também podemos usar técnicas bayesianas para resolver os mesmos problemas. A vantagem de se fazer isso, como veremos a diante, é o fato que nunca deixamos de considerar nossa incerteza sobre nossas inferências, e também somos capazes muitas vezes de adicionar conhecimento a priori sobre o problema para delimitar sua solução. As técnicas de programação probabilística que introduzimos nesse post são capazes de modelar e resolver problemas desta maneira de forma simples e prática. Usaremos para isso a biblioteca PyMC3, que usa o Theano para fazer inferência neste tipo de modelo.

     

    Regressão Logística

    Um dos problemas mais simples que podemos querer resolver é descobrir a relação linear entre duas variáveis $x_1$ e $x_2$ e um target $y$. Em geral, essa relação pode ser escrita como:

    $$\mu = \alpha + \beta_1 x_1 + \beta_2 x_2$$

    No entanto, em geral, temos ainda ruído sobre essa relação, o que leva a distribuição de $y$ a ser escrita como:

    $$y = \mathcal{N}(\mu, \sigma)$$

    Onde sigma é uma variância que está ligada a este ruído. Resolver este problema é equivalente a estimar $\alpha$, $beta_1$, $beta_2$ e $sigma$, de forma que consigamos, dados $x_1$ e $x_2$, estimar o valor correspondente de $y$.

    Esta estimativa pode ser feita a partir de um conjunto de dados de algumas formas. Podemos estimar os parâmetros usando máximo a posteriori, ou seja, escolher os parâmetros de forma que a maximizar a probabilidade dos dados. Podemos também ter como resultado a distribuição dos parâmetros a posteriori, o que nos permite falar da incerteza que temos sobre os parâmetros encontrados. Esta distribuição é $P(\lambda\|D)$, e pode ser estimada usando técnicas de Monte Carlo ou resolvendo o problema analiticamente. Neste post mostraremos como estimar o máximo a posteriori e a distribuição a posteriori usando Monte Carlo, fazendo uso do PyMC3, uma ferramenta para descrever e resolver estes problemas usando a linguagem de programação Python.

     

    Dados sintéticos

    Podemos gerar dados sintéticos, fazendo uso do `numpy`:

    from numpy import random
    from matplotlib import pyplot as plt
    from mpl_toolkits.mplot3d import Axes3D

    # Valores reais dos parametros
    alpha = 1
    sigma = 0.1
    beta = [3, 2.5]

    # Tamanho do dataset
    size = 1000

    # Dados de entrada
    X1 = random.randn(size)
    X2 = random.randn(size) * 0.2

    # Simular variável de target
    Y = alpha + beta[0]*X1 + beta[1]*X2 + random.randn(size)*sigma

    # Vizualization
    fig = plt.figure()
    ax = Axes3D(fig)
    ax.scatter(X1, X2, Y)

     

    Modelando o problema com PyMC3

    Suponha que tenhamos um conhecimento a prior sobre os parâmetros que desejamos estimar:

    $$\alpha \tilde \mathcal{N}(0, 100)$$
    $$\beta_1 \tilde \mathcal{N}(0, 100)$$
    $$\beta_2 \tilde \mathcal{N}(0, 100)$$
    $$\sigma \tilde \|\mathcal{N}(0, 1)\|$$

    Podemos específicar nosso modelo usando PyMC3 da seguinte maneira:

    import pymc3

    with pymc3.Model() as linear_model:
    # Priors para os parâmetros
    alpha = pymc3.Normal('alpha', mu=0, sd=10)
    beta = pymc3.Normal('beta', mu=0, sd=10, shape=2)
    sigma = pymc3.HalfNormal('sigma', sd=1)

    # Valor esperado do target
    mu = alpha + beta[0]*X1 + beta[1]*X2

    # Likelihood das observações
    Y_obs = pymc3.Normal('Y_obs', mu=mu, sd=sigma, observed=Y)

    Este código é o suficiente para especificar nosso modelo: `pymc3.Model()` define um container para as variáveis do modelo, enquanto as outras linhas definem o modelo de maneira quase idêntica a notação matemática que usamos até agora. Note que cada variável possui um nome e pode ser tanto variáveis estocásticas (definidas em função de uma distribuição, como nossos priors) quanto determinísticas (cujo valor pode ser determinado completamente como uma função de outras variáveis, como `mu` em nosso caso).

    A definição do target também pode receber como parâmetro um vetor de observações. Este vetor será útil ao estimarmos os parâmetros como MAP por exemplo, o que faremos na próxima seção.

     

    Máximo a posteriori

    No PyMC3 é quase trivial obter o MAP dos parâmetros de um modelo dadas algumas observações. Após definido o modelo, podemos fazer:

    map_estimate = pymc3.find_MAP(model=linear_model)

    E encontramos algo como

    {
    'alpha': array(1.0063472447160637),
    'beta': array([ 3.00424673, 2.51779307]),
    'sigma': array(0.09715281616720224),
    'sigma_log__': array(-2.331470115753791)
    }

    A estimativa é feita usando algoritmos de otimização como o BFGS para encontrar o máximo do log da distribuição posterior. No entanto, encontrar o MAP pode nem sempre ser a melhor ideia, já que podemos ter múltiplos picos na distribuição, ou regiões com maior probabilidade bem longe do pico. Outro problema é que perdemos a noção sobre a incerteza desta estimativa. Métodos de Monte Carlo podem ser úteis para resolver este tipo de problema, como veremos na seção a seguir.

     

    Métodos de sampling

    Podemos simular a distribuição dos parâmetros dadas as observações usando método de Monte Carlo. Apesar de a explicação do Monte Carlo estar fora do escopo deste artigo, podemos entender esta técnica como um sequência de amostras dos parâmetros que, aplicando algumas regras envolvendo a probabilidade das observações, converge para a distribuição a posteriori. Para calcular isso, basta fazer:

    with linear_model:
    # 1000 amostras do posterior
    trace = pymc3.sample(draws=1000)

    # Plotar o posterior
    pymc3.traceplot(trace)

    Note como agora temos uma distribuição de probabilidades sobre cada um dos parâmetros. Agora podemos ter uma ideia de como a nossa incerteza sobre esses parâmetros está distribuída.

     

    Conclusão

    Em geral, programação probabilística é ideal para fazer inferência sobre modelos que você consegue escrever a distribuição de cada uma das variáveis. Embora fazer inferência desta maneira possa ser computacionalmente cara, em geral estes modelos podem ajudar ao colocar conhecimento de domínio específico (conhecimento sobre distribuições ou sobre relações das variáveis). O PyMC3 é uma ótima ferramenta para se familiarizar com o assunto (assim como o livro Probabilistic Programming for Hackers), e bibliotecas mais poderosas como Edward ou Pyro podem ser usadas para fazer inferência em modelos bem mais complicados. No DataLab, usamos modelos probabilísticos em ferramentas como o Dice, onde modelamos a relação entre variáveis como um modelo gráfico probabilístico.

    1 de dezembro de 2017

    Filtre por mês

    Filtre por autor