Bayesian_Switchpoint_Analysis.ipynb

Kernel: Python 3

Copyright 2019 The TensorFlow Probability Authors.

Licensed under the Apache License, Version 2.0 (the "License");

#@title Licensed under the Apache License, Version 2.0 (the "License"); { display-mode: "form" }
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# https://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

Análise bayesiana de pontos de alteração

Veja em TensorFlow.org

Executar no Google Colab

Ver fonte em GitHub

Baixar notebook

Este notebook reimplementa e estende o exemplo “Análise bayesiana de pontos de mudança” da documentação de pymc3.

Pré-requisitos

import tensorflow.compat.v2 as tf
tf.enable_v2_behavior()
import tensorflow_probability as tfp
tfd = tfp.distributions
tfb = tfp.bijectors
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = (15,8)
%config InlineBackend.figure_format = 'retina'
import numpy as np
import pandas as pd

Dataset

O dataset é este aqui. Observação: há outra versão deste exemplo disponível em outro lugar, mas há dados “ausentes” e, nesse caso, você precisaria incluir os valores ausentes (caso contrário, seu modelo jamais sairá dos parâmetros iniciais porque a função de verossimilhança não estará definida).

disaster_data = np.array([ 4, 5, 4, 0, 1, 4, 3, 4, 0, 6, 3, 3, 4, 0, 2, 6,
                           3, 3, 5, 4, 5, 3, 1, 4, 4, 1, 5, 5, 3, 4, 2, 5,
                           2, 2, 3, 4, 2, 1, 3, 2, 2, 1, 1, 1, 1, 3, 0, 0,
                           1, 0, 1, 1, 0, 0, 3, 1, 0, 3, 2, 2, 0, 1, 1, 1,
                           0, 1, 0, 1, 0, 0, 0, 2, 1, 0, 0, 0, 1, 1, 0, 2,
                           3, 3, 1, 1, 2, 1, 1, 1, 1, 2, 4, 2, 0, 0, 1, 4,
                           0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1])
years = np.arange(1851, 1962)
plt.plot(years, disaster_data, 'o', markersize=8);
plt.ylabel('Disaster count')
plt.xlabel('Year')
plt.title('Mining disaster data set')
plt.show()

Modelo probabilístico

O modelo pressupõe um “ponto de alteração” (por exemplo: um ano durante o qual as regulações de segurança mudaram) e a taxa de desastre com distribuição de Poisson com taxas constantes (mas potencialmente diferentes) antes e após esse ponto de alteração.

A contagem de desastres real é fixa (observada); toda amostra desse modelo precisará especificar tanto o ponto de alteração quanto a taxa de desastres “inicial” e “posterior”.

Modelo original no exemplo da documentação de pymc3:

$$\begin{align*} (D_t|s,e,l)&\sim \text{Poisson}(r_t), \\ & \,\quad\text{with}\; r_t = \begin{cases}e & \text{if}\; t < s\\l &\text{if}\; t \ge s\end{cases} \\ s&\sim\text{Discrete Uniform}(t_l,\,t_h) \\ e&\sim\text{Exponential}(r_e)\\ l&\sim\text{Exponential}(r_l) \end{align*}$$

Entretanto, a taxa de desastres média $r_t$ tem uma descontinuidade no ponto de alteração $s$, o que a torna diferenciável. Portanto, ela não fornece nenhum sinal de gradiente ao algoritmo Monte Carlo Hamiltoniano (HMC) – porém, como o prior $s$ é contínuo, o fallback do HMC para um random-walk é suficiente para encontrar as áreas de massa de alta probabilidade neste exemplo.

Em um segundo modelo, modificamos o original usando a “alteração” de sigmoide entre e e l para tornar a transição diferenciável e usamos uma distribuição uniforme contínua para o ponto de alteração $s$ (é possível argumentar que este modelo é mais fiel à realidade, pois uma “alteração” na taxa média provavelmente seria distribuída ao longo de vários anos). Portanto, o novo modelo é:

$$\begin{align*} (D_t|s,e,l)&\sim\text{Poisson}(r_t), \\ & \,\quad \text{with}\; r_t = e + \frac{1}{1+\exp(s-t)}(l-e) \\ s&\sim\text{Uniform}(t_l,\,t_h) \\ e&\sim\text{Exponential}(r_e)\\ l&\sim\text{Exponential}(r_l) \end{align*}$$

Na ausência de mais informações, pressupomos que $r_e = r_l = 1$ como parâmetros para os priors. Vamos executar os dois modelos e comparar os resultados da inferência.

def disaster_count_model(disaster_rate_fn):
  disaster_count = tfd.JointDistributionNamed(dict(
    e=tfd.Exponential(rate=1.),
    l=tfd.Exponential(rate=1.),
    s=tfd.Uniform(0., high=len(years)),
    d_t=lambda s, l, e: tfd.Independent(
        tfd.Poisson(rate=disaster_rate_fn(np.arange(len(years)), s, l, e)),
        reinterpreted_batch_ndims=1)
  ))
  return disaster_count

def disaster_rate_switch(ys, s, l, e):
  return tf.where(ys < s, e, l)

def disaster_rate_sigmoid(ys, s, l, e):
  return e + tf.sigmoid(ys - s) * (l - e)

model_switch = disaster_count_model(disaster_rate_switch)
model_sigmoid = disaster_count_model(disaster_rate_sigmoid)

O código acima define o modelo por meio de distribuições JointDistributionSequential. As funções disaster_rate são chamadas com um array de [0, ..., len(years)-1] para produzir um vetor de len(years) variáveis aleatórias – os anos antes do switchpoint (ponto de alteração) são early_disaster_rate, e os anos após, são late_disaster_rate (módulo da transição sigmoide).

Aqui está a prova real de que a função de probabilidade logarítmica é verdadeira:

def target_log_prob_fn(model, s, e, l):
  return model.log_prob(s=s, e=e, l=l, d_t=disaster_data)

models = [model_switch, model_sigmoid]
print([target_log_prob_fn(m, 40., 3., .9).numpy() for m in models])  # Somewhat likely result
print([target_log_prob_fn(m, 60., 1., 5.).numpy() for m in models])  # Rather unlikely result
print([target_log_prob_fn(m, -10., 1., 1.).numpy() for m in models]) # Impossible result

[-176.94559, -176.28717]
[-371.3125, -366.8816]
[-inf, -inf]

HMC para fazer inferência bayesiana

Definimos o número de resultados e passos de burn-in necessários; o código é modelado principalmente de acordo com a documentação de tfp.mcmc.HamiltonianMonteCarlo e usa um tamanho adaptativo de passo (caso contrário, o resultado é muito sensível ao valor do tamanho do passo escolhido). Usamos o valor 1 como estado inicial da cadeia.

Mas essa não é a história completa. Se você rever a definição do modelo acima, verá que algumas das distribuições de probabilidade não são bem definidas em todo o eixo de números reais. Portanto, restringimos o espaço para que o HMC examine por meio do encapsulamento do kernel do HMC com um TransformedTransitionKernel, que especifica os bijetores forward para transformar os números reais para o domínio no qual a distribuição de probabilidade está definida (conforme os comentários no código abaixo).

num_results = 10000
num_burnin_steps = 3000

@tf.function(autograph=False, jit_compile=True)
def make_chain(target_log_prob_fn):
   kernel = tfp.mcmc.TransformedTransitionKernel(
       inner_kernel=tfp.mcmc.HamiltonianMonteCarlo(
          target_log_prob_fn=target_log_prob_fn,
          step_size=0.05,
          num_leapfrog_steps=3),
       bijector=[
          # The switchpoint is constrained between zero and len(years).
          # Hence we supply a bijector that maps the real numbers (in a
          # differentiable way) to the interval (0;len(yers))
          tfb.Sigmoid(low=0., high=tf.cast(len(years), dtype=tf.float32)),
          # Early and late disaster rate: The exponential distribution is
          # defined on the positive real numbers
          tfb.Softplus(),
          tfb.Softplus(),
      ])
   kernel = tfp.mcmc.SimpleStepSizeAdaptation(
        inner_kernel=kernel,
        num_adaptation_steps=int(0.8*num_burnin_steps))

   states = tfp.mcmc.sample_chain(
      num_results=num_results,
      num_burnin_steps=num_burnin_steps,
      current_state=[
          # The three latent variables
          tf.ones([], name='init_switchpoint'),
          tf.ones([], name='init_early_disaster_rate'),
          tf.ones([], name='init_late_disaster_rate'),
      ],
      trace_fn=None,
      kernel=kernel)
   return states

switch_samples = [s.numpy() for s in make_chain(
    lambda *args: target_log_prob_fn(model_switch, *args))]
sigmoid_samples = [s.numpy() for s in make_chain(
    lambda *args: target_log_prob_fn(model_sigmoid, *args))]

switchpoint, early_disaster_rate, late_disaster_rate = zip(
    switch_samples, sigmoid_samples)

Execute os dois modelos em paralelo:

Visualize o resultado

Visualizamos o resultado como histogramas das amostras da distribuição posterior para a taxa de desastre inicial e posterior, bem como o ponto de alteração. É adicionada aos histogramas uma linha sólida representando a média das amostras, bem como o 95º percentil dos limites do intervalo de credibilidade como linhas tracejadas.

def _desc(v):
  return '(median: {}; 95%ile CI: $[{}, {}]$)'.format(
      *np.round(np.percentile(v, [50, 2.5, 97.5]), 2))

for t, v in [
    ('Early disaster rate ($e$) posterior samples', early_disaster_rate),
    ('Late disaster rate ($l$) posterior samples', late_disaster_rate),
    ('Switch point ($s$) posterior samples', years[0] + switchpoint),
]:
  fig, ax = plt.subplots(nrows=1, ncols=2, sharex=True)
  for (m, i) in (('Switch', 0), ('Sigmoid', 1)):
    a = ax[i]
    a.hist(v[i], bins=50)
    a.axvline(x=np.percentile(v[i], 50), color='k')
    a.axvline(x=np.percentile(v[i], 2.5), color='k', ls='dashed', alpha=.5)
    a.axvline(x=np.percentile(v[i], 97.5), color='k', ls='dashed', alpha=.5)
    a.set_title(m + ' model ' + _desc(v[i]))
  fig.suptitle(t)
  plt.show()