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.

Copula 입문

TensorFlow.org에서보기

Google Colab에서 실행

GitHub에서 소그 보기

노트북 다운로드

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

[Copula](https://en.wikipedia.org/wiki/Copula_(probability_theory))는 확률 변수 간의 종속성을 캡처하기 위한 고전적인 접근 방식입니다. 더 공식적으로 코퓰러는 다변량 분포 $C(U_1, U_2, ...., U_n)$이므로 주변화(marginalizing)를 통해 $U_i \sim \text{Uniform}(0, 1)$을 얻게 됩니다.

코퓰러는 임의의 한계가 있는 다변량 분포를 만드는 데 사용할 수 있으므로 흥미롭습니다. 다음은 공식입니다.

• 확률 적분 변환을 사용하면 임의의 연속 R.V. $X$를 균일한 $F_X(X)$로 변환합니다. 여기서 $F_X$는 $X$의 CDF입니다.

• 코퓰러(이변량) $C(U, V)$가 주어지면 $U$와 $V$가 균일한 주변 분포를 갖게 됩니다.

• 이제 관심 있는 R.V의 $X, Y$가 주어지면 새 분포 $C'(X, Y) = C(F_X(X), F_Y(Y))$를 만듭니다. $X$와 $Y$의 한계는 원하는 한계입니다.

한계는 일변량이므로 측정 및/또는 모델링이 더 쉬울 수 있습니다. 코퓰러는 한계에서 시작하면서 차원 간의 임의의 상관관계를 달성할 수 있습니다.

가우시안 코퓰러

코퓰러가 어떻게 구성되는지 설명하려면 다변량 가우시안 상관관계에 따라 의존성을 캡처하는 경우를 고려하세요. 가우시안 코퓰러로는 $C(u_1, u_2, ...u_n) = \Phi_\Sigma(\Phi^{-1}(u_1), \Phi^{-1}(u_2), ... \Phi^{-1}(u_n))$가 주어졌습니다. 여기서 $\Phi_\Sigma$는 공분산 $\Sigma$와 평균 0을 가진 MultivariateNormal의 CDF를 나타내고, $\Phi^{-1}$는 표준 정규 분포를 위한 역 CDF입니다.

정규 분포의 역 CDF를 적용하면 균일한 차원이 정규 분포로 왜곡됩니다. 그런 다음 다변량 정규 분포의 CDF를 적용하면 분포가 어느 정도 균일하고 가우시안 상관관계를 갖도록 압축됩니다.

따라서 얻게 된 점은 가우시안 코퓰러가 균일한 한계를 갖는 단위 하이퍼큐브 $\LaTeX$[0, 1]^n$에 대한 분포라는 것입니다.

이와 같이 정의된 가우시안 코퓰러는 tfd.TransformedDistribution 및 적절한 Bijector로 구현할 수 있습니다. 즉, tfb.NormalCDF bijector로 구현된, 정규 분포의 역 CDF를 사용하여 MultivariateNormal을 변환합니다.

아래에서는 공분산이 콜레스키(Cholesky) 인자(따라서 MultivariateNormalTriL 대한 공분산)에 의해 매개변수화된다는 단순화 가정을 사용하여 가우시안 코퓰러를 구현합니다(다른 tf.linalg.LinearOperators를 사용하여 행렬이 없는 다양한 가정을 인코딩할 수 있음).

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

하지만 이러한 모델의 힘은 확률 적분 변환을 사용하여 임의의 R.V.에서 코퓰러를 사용하는 것입니다. 이런 방식으로 임의의 한계를 지정하고 코퓰러를 사용하여 함께 연결할 수 있습니다.

다음 모델로 시작합니다.

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

코퓰러를 사용하여 한계 Kumaraswamy 및 Gumbel을 갖는 이변량 R.V. $Z$를 얻습니다.

이들 두 R.V.로 생성된 제품 분포를 플롯하여 시작하겠습니다. 이는 코퓰러를 적용할 때 비교 지점 역할을 하기 위한 것입니다.

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

다양한 한계를 가진 결합 분포

이제 가우시안 코퓰러를 사용하여 분포를 함께 연결하고 플롯합니다. 선택한 도구는 TransformedDistribution이며, 이를 통해 적절한 Bijector를 적용하여 선택한 한계를 얻습니다.

특히 벡터의 여러 부분에 다양한 bijector를 적용하는 Blockwise bijector를 사용합니다(여전히 전단사 변환임).

이제 원하는 코퓰러를 정의할 수 있습니다. 목표 한계 목록(bijector로 인코딩됨)이 주어지면 코퓰러를 사용하고 지정된 한계를 갖는 새 분포를 쉽게 구성할 수 있습니다.

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

마지막으로 이 가우시안 코퓰러를 실제로 사용해 봅시다. $\begin{bmatrix}1 & 0\\rho & \sqrt{(1-\rho^2)}\end{bmatrix}$의 콜레스키를 사용합니다. 이는 분산 1에 해당하며, 다변량 정규 분포인 경우 상관관계 $\rho$에 해당합니다.

몇 가지 경우를 살펴보겠습니다.

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

마지막으로 원하는 한계를 실제로 얻었는지 확인해봅니다.

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

결론

이제 됐습니다. Bijector API로 가우시안 코퓰러를 생성할 수 있음을 입증했습니다.

더 일반적으로는 Bijector API로 bijector를 작성하고 이를 분포로 구성하면 유연한 모델링을 위한 풍부한 분포군을 만들 수 있습니다.