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.

Introducción a las cópulas

Ver en TensorFlow.org

Ejecutar en Google Colab

Ver fuente en GitHub

Descargar el bloc de notas

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

Una [cópula](https://en.wikipedia.org/wiki/Copula_(probability_theory)) es un enfoque clásico para capturar la dependencia entre variables aleatorias. Más formalmente, una cópula es una distribución multivariada $C(U_1, U_2 , ...., U_n)$ tal que al marginar se obtiene $U_i \sim \text{Uniform}(0, 1)$.

Las cópulas son interesantes porque podemos usarlas para crear distribuciones multivariadas con marginales arbitrarios. Estos son los pasos a seguir:

• El uso de la transformación integral de probabilidad convierte una RV $X$ continua arbitraria en una uniforme $F_X(X)$, donde $F_X$ es el CDF de $X$.

• Dada una cópula (digamos bivariada) $C(U, V)$, tenemos que $U$ y $V$ tienen distribuciones marginales uniformes.

• Ahora, dadas nuestras RV de interés $X, Y$, cree una nueva distribución $C'(X, Y) = C(F_X(X), F_Y(Y))$. Los marginales para $X$ y $Y$ son los que deseábamos.

Los marginales son univariados y por lo tanto pueden ser más fáciles de medir y/o modelar. Una cópula permite partir de marginales, pero también lograr una correlación arbitraria entre dimensiones.

Cópula gaussiana

Para ilustrar cómo se construyen las cópulas, consideremos el caso de capturar la dependencia según correlaciones gaussianas multivariadas. Una cópula gaussiana es aquella dada por $C(u_1, u_2, ...u_n) = \Phi_\Sigma(\Phi^{-1}(u_1), \Phi^{-1}(u_2), ... \Phi^{-1}(u_n))$ donde $\Phi_\Sigma$ representa la CDF de un MultivariateNormal, con covarianza $\Sigma$ y media 0, y $\Phi^{-1}$ es la CDF inversa de la normal estándar.

La aplicación de la CDF inversa de la normal deforma las dimensiones uniformes para que se distribuyan normalmente. La aplicación de la CDF de la normal multivariada luego aplasta la distribución para que sea marginalmente uniforme y con correlaciones gaussianas.

De este modo, lo que obtenemos es que la cópula gaussiana es una distribución sobre el hipercubo unitario $[0, 1]^n$ con marginales uniformes.

Definida como tal, la cópula gaussiana se puede implementar con tfd.TransformedDistribution y Bijector apropiados. Es decir, estamos transformando una MultivariateNormal, mediante el uso de la CDF inversa de la distribución normal, implementada por el biyector tfb.NormalCDF.

A continuación, implementamos una cópula gaussiana con una suposición simplificadora: que la covarianza e parametriza mediante un factor de Cholesky (de ahí una covarianza para MultivariateNormalTriL). (Se podrían usar otros tf.linalg.LinearOperators para codificar diferentes suposiciones sin matrices).

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);

Sin embargo, el poder de tal modelo es usar la Transformada Integral de Probabilidad, para usar la cópula en RV arbitrarias. De esta manera, podemos especificar marginales arbitrarios y usar la cópula para unirlos.

Empezamos con un modelo:

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

y usamos la cópula para obtener una RV $Z$ bivariada, que tiene marginales Kumaraswamy y Gumbel.

Para empezar, trazaremos la distribución de productos generada por esas dos variables aleatorias. Esto es solo para que sirva como punto de comparación cuando aplicamos la 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);

Distribución conjunta con diferentes marginales

Ahora usamos una cópula gaussiana para unir las distribuciones y trazarla. Nuevamente, nuestra herramienta preferida es TransformedDistribution y aplicamos el Bijector apropiado para obtener los marginales elegidos.

Específicamente, utilizamos un biyector Blockwise que aplica diferentes biyectores en diferentes partes del vector (que sigue siendo una transformación biyectiva).

Ahora podemos definir la cópula que queremos. Dada una lista de marginales objetivo (codificados como biyectores), podemos construir fácilmente una nueva distribución que use la cópula y tenga los marginales especificados.

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")

Finalmente, usemos esta cópula gaussiana. Usaremos un factor de Cholesky de $\begin{bmatrix}1 & 0\\\rho & \sqrt{(1-\rho^2)}\end{bmatrix}$, que corresponderá a las varianzas 1 y la correlación $\rho$ para la normal multivariada.

Veremos algunos 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))

Finalmente, verifiquemos que realmente obtengamos los marginales que queremos.

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))
  

Conclusión

¡Y ahí está! Hemos demostrado que podemos construir cópulas gaussianas mediante el uso de la API Bijector.

De manera más general, escribir biyectores con la API Bijector y componerlos con una distribución puede crear complejas familias de distribuciones para un modelado flexible.