Gaussian_Copula.ipynb

Kernel: Python 3

Copyright 2018 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.

Guia sobre cópulas

Veja em TensorFlow.org

Executar no Google Colab

Ver fonte em GitHub

Baixar notebook

import numpy as np
import matplotlib.pyplot as plt
import tensorflow.compat.v2 as tf
tf.enable_v2_behavior()

import tensorflow_probability as tfp
tfd = tfp.distributions
tfb = tfp.bijectors

Uma [cópula](https://pt.wikipedia.org/wiki/Cópula_(estatística) é uma estratégia básica para obter a dependência entre variáveis aleatórias. Mais formalmente, uma cópula é uma distribuição multivariada $C(U_1, U_2, ...., U_n)$ em que a marginalização fornece $U_i \sim \text{Uniform}(0, 1)$.

As cópulas são interessantes porque podemos usá-las para criar distribuições multivariadas com distribuições marginais arbitrárias. Esta é a receita:

• Usar a Transformação Integral da Probabilidade transforma uma variável aleatória contínua arbitrária $X$ em outra uniforme $F_X(X)$, em que $F_X$ é a função de distribuição acumulada (CDF, na sigla em inglês) de $X$.

• Dada uma cópula (bivariada, digamos) $C(U, V)$, $U$ e $V$ têm distribuições marginais uniformes.

• Agora, dadas nossas variáveis aleatórias de interesse $X, Y$, crie uma nova distribuição $C'(X, Y) = C(F_X(X), F_Y(Y))$. As distribuições marginais para $X$ e $Y$ são as que desejamos.

As distribuições marginais são univariadas e, portanto, podem ser mais fáceis de mensurar e/ou modelar. Uma cópula permite começar pelas distribuições marginais, mas ainda assim obter a correlação arbitrária dentre as dimensões.

Cópula gaussiana

Para ilustrar como as cópulas são construídas, considere o caso de obter a dependência segundo as correlações gaussianas multivariadas. Uma cópula gaussiana é aquela fornecida por $C(u_1, u_2, ...u_n) = \Phi_\Sigma(\Phi^{-1}(u_1), \Phi^{-1}(u_2), ... \Phi^{-1}(u_n))$, em que $\Phi_\Sigma$ representa a função de distribuição acumulada de uma distribuição normal multivariada (MultivariateNormal), com covariância $\Sigma$ e média 0, e $\Phi^{-1}$ é a função de distribuição acumulada inversa para a distribuição normal padrão.

Aplicar a CDF inversa da distribuição normal deforma as dimensões uniformes de modo a terem distribuição normal. Em seguida, aplicar a CDF da distribuição normal multivariada espreme a distribuição de forma a se tornar uniforme marginalmente e com correlações gaussianas.

Portanto, a cópula gaussiana é uma distribuição ao longo do hipercubo unitário $[0, 1]^n$ com distribuições marginais uniformes.

Definida desta forma, a cópula gaussiana pode ser implementada com tfd.TransformedDistribution e um Bijector apropriado. Ou seja, estamos transformando uma distribuição normal multivariada pelo uso da CDF inversa da distribuição normal, implementada pelo bijetor tfb.NormalCDF.

Abaixo, implementamos uma cópula gaussiana com uma suposição simplificadora: que a covariância é parametrizada por um fator de Cholesky (portanto, uma covariância de MultivariateNormalTriL). Também é possível usar outros tf.linalg.LinearOperators para codificar diferentes suposições sem matrizes.

class GaussianCopulaTriL(tfd.TransformedDistribution):
  """Takes a location, and lower triangular matrix for the Cholesky factor."""
  def __init__(self, loc, scale_tril):
    super(GaussianCopulaTriL, self).__init__(
        distribution=tfd.MultivariateNormalTriL(
            loc=loc,
            scale_tril=scale_tril),
        bijector=tfb.NormalCDF(),
        validate_args=False,
        name="GaussianCopulaTriLUniform")


# Plot an example of this.
unit_interval = np.linspace(0.01, 0.99, num=200, dtype=np.float32)
x_grid, y_grid = np.meshgrid(unit_interval, unit_interval)
coordinates = np.concatenate(
    [x_grid[..., np.newaxis],
     y_grid[..., np.newaxis]], axis=-1)

pdf = GaussianCopulaTriL(
    loc=[0., 0.],
    scale_tril=[[1., 0.8], [0., 0.6]],
).prob(coordinates)

# Plot its density.

plt.contour(x_grid, y_grid, pdf, 100, cmap=plt.cm.jet);

Portanto, o poder de um modelo como esse é usar a Transformação Integral da Probabilidade para utilizar a cópula em variáveis aleatórias arbitrárias. Dessa forma, podemos especificar distribuições marginais arbitrárias e usar a cópula para uni-las.

Começamos com o seguinte modelo:

$$\begin{align*} X &\sim \text{Kumaraswamy}(a, b) \\ Y &\sim \text{Gumbel}(\mu, \beta) \end{align*}$$

E depois usamos a cópula para obter uma variável aleatória bivariada $Z$, que tem distribuições marginais Kumaraswamy e Gumbel.

Vamos começar plotando a distribuição do produto gerada por essas duas variáveis aleatórias apenas para servir como ponto de comparação quando aplicarmos a cópula.

a = 2.0
b = 2.0
gloc = 0.
gscale = 1.

x = tfd.Kumaraswamy(a, b)
y = tfd.Gumbel(loc=gloc, scale=gscale)

# Plot the distributions, assuming independence
x_axis_interval = np.linspace(0.01, 0.99, num=200, dtype=np.float32)
y_axis_interval = np.linspace(-2., 3., num=200, dtype=np.float32)
x_grid, y_grid = np.meshgrid(x_axis_interval, y_axis_interval)

pdf = x.prob(x_grid) * y.prob(y_grid)

# Plot its density

plt.contour(x_grid, y_grid, pdf, 100, cmap=plt.cm.jet);

Distribuição conjunta com distribuições marginais diferentes

Agora, usamos uma cópula gaussiana para unir as distribuições e plotamos o resultado. Novamente, escolhemos a ferramenta TransformedDistribution, aplicando o Bijector adequado para obter as distribuições marginais escolhidas.

Especificamente, usamos o bijetor Blockwise, que aplica diferentes bijetores em diferentes partes do vetor (isso ainda é uma transformação bijetora).

Agora, podemos definir a cópula que desejamos. Dada uma lista de distribuições marginais alvo (codificadas como bijetores), podemos construir facilmente uma nova distribuição que use a cópula e tenha as distribuições marginais especificadas.

class WarpedGaussianCopula(tfd.TransformedDistribution):
  """Application of a Gaussian Copula on a list of target marginals.

  This implements an application of a Gaussian Copula. Given [x_0, ... x_n]
  which are distributed marginally (with CDF) [F_0, ... F_n],
  `GaussianCopula` represents an application of the Copula, such that the
  resulting multivariate distribution has the above specified marginals.

  The marginals are specified by `marginal_bijectors`: These are
  bijectors whose `inverse` encodes the CDF and `forward` the inverse CDF.

  block_sizes is a 1-D Tensor to determine splits for `marginal_bijectors`
  length should be same as length of `marginal_bijectors`.
  See tfb.Blockwise for details
  """
  def __init__(self, loc, scale_tril, marginal_bijectors, block_sizes=None):
    super(WarpedGaussianCopula, self).__init__(
        distribution=GaussianCopulaTriL(loc=loc, scale_tril=scale_tril),
        bijector=tfb.Blockwise(bijectors=marginal_bijectors,
                               block_sizes=block_sizes),
        validate_args=False,
        name="GaussianCopula")

Por fim, vamos usar a cópula gaussiana. Usaremos um Cholesky de $\begin{bmatrix}1 & 0\\\rho & \sqrt{(1-\rho^2)}\end{bmatrix}$, que corresponderá a variâncias 1 e correlação $\rho$ para a distribuição normal multivariada.

Vamos conferir alguns casos:

# Create our coordinates:
coordinates = np.concatenate(
    [x_grid[..., np.newaxis], y_grid[..., np.newaxis]], -1)


def create_gaussian_copula(correlation):
  # Use Gaussian Copula to add dependence.
  return WarpedGaussianCopula(
      loc=[0.,  0.],
      scale_tril=[[1., 0.], [correlation, tf.sqrt(1. - correlation ** 2)]],
      # These encode the marginals we want. In this case we want X_0 has
      # Kumaraswamy marginal, and X_1 has Gumbel marginal.

      marginal_bijectors=[
          tfb.Invert(tfb.KumaraswamyCDF(a, b)),
          tfb.Invert(tfb.GumbelCDF(loc=0., scale=1.))])


# Note that the zero case will correspond to independent marginals!
correlations = [0., -0.8, 0.8]
copulas = []
probs = []
for correlation in correlations:
  copula = create_gaussian_copula(correlation)
  copulas.append(copula)
  probs.append(copula.prob(coordinates))


# Plot it's density

for correlation, copula_prob in zip(correlations, probs):
  plt.figure()
  plt.contour(x_grid, y_grid, copula_prob, 100, cmap=plt.cm.jet)
  plt.title('Correlation {}'.format(correlation))

Por fim, vamos verificar se obtivemos as distribuições marginais que desejamos.

def kumaraswamy_pdf(x):
    return tfd.Kumaraswamy(a, b).prob(np.float32(x))

def gumbel_pdf(x):
    return tfd.Gumbel(gloc, gscale).prob(np.float32(x))


copula_samples = []
for copula in copulas:
  copula_samples.append(copula.sample(10000))

plot_rows = len(correlations)
plot_cols = 2  # for 2  densities [kumarswamy, gumbel]
fig, axes = plt.subplots(plot_rows, plot_cols, sharex='col', figsize=(18,12))

# Let's marginalize out on each, and plot the samples.

for i, (correlation, copula_sample) in enumerate(zip(correlations, copula_samples)):
  k = copula_sample[..., 0].numpy()
  g = copula_sample[..., 1].numpy()


  _, bins, _ = axes[i, 0].hist(k, bins=100, density=True)
  axes[i, 0].plot(bins, kumaraswamy_pdf(bins), 'r--')
  axes[i, 0].set_title('Kumaraswamy from Copula with correlation {}'.format(correlation))

  _, bins, _ = axes[i, 1].hist(g, bins=100, density=True)
  axes[i, 1].plot(bins, gumbel_pdf(bins), 'r--')
  axes[i, 1].set_title('Gumbel from Copula with correlation {}'.format(correlation))
  

Conclusão

Pronto! Demonstramos que podemos criar cópulas gaussianas usando a API Bijector.

De forma mais geral, ao escrever bijetores usando a API Bijector e uni-los com uma distribuição, é possível criar famílias avançadas de distribuições, o que proporciona uma modelagem flexível.