コピュラ入門書

TensorFlow.orgで表示 GoogleColabで実行 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

[コピュラ](https://en.wikipedia.org/wiki/Copula_(probability_theory%29)は、確率変数間の依存関係をキャプチャするための古典的なアプローチです。より正式には、コピュラは多変量分布$ C(U_1、U_2 、....、U_n)$は、マージナライズによって$ U_i \ sim \ text {Uniform}(0、1)$が得られるようにします。

コピュラは、任意の周辺分布を持つ多変量分布を作成するために使用できるため、興味深いものです。これがレシピです:

  • 確率積分変換を使用すると、任意の連続RV $ X $が均一な$ F_X(X)$に変わります。ここで、$ F_X $は$ X $のCDFです。
  • コピュラ(たとえば2変量)$ C(U、V)$が与えられると、$ U $と$ V $の周辺分布は均一になります。
  • ここで、関心のあるRV $ 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を適用すると、分布が押しつぶされてわずかに均一になり、ガウス相関が得られます。

したがって、ガウスコピュラは、均一な周辺分布を持つ単位超立方体$ [0、1] ^ n $にわたる分布であることがわかります。

そのように定義されたガウスコピュラは、 tfd.TransformedDistributionと適切なBijectortfd.TransformedDistributionして実装できます。つまり、 tfb.NormalCDFバイジェクターによって実装された正規分布の逆CDFを使用して、MultivariateNormalを変換しています。

以下では、共分散がコレスキー因子によってパラメーター化されるという1つの単純化された仮定でガウスコピュラを実装します(したがって、 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);

png

ただし、このようなモデルの能力は、確率積分変換を使用して、任意のRVでコピュラを使用することです。このようにして、任意の周辺分布を指定し、コピュラを使用してそれらをつなぎ合わせることができます。

モデルから始めます。

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

コピュラを使用して、周辺分布のKumaraswamyGumbelを持つ2変量RV $ Z $を取得します。

これらの2つのRVによって生成された製品分布をプロットすることから始めます。これは、コピュラを適用するときの比較ポイントとして機能するためのものです。

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

png

異なる周辺分布との同時分布

ここで、ガウスコピュラを使用して分布を結合し、それをプロットします。ここでも、選択したツールは、適切なBijectorを適用して選択した周辺Bijectorを取得するTransformedDistributionです。

具体的には、ベクトルのさまざまな部分にさまざまな全単射を適用するBlockwiseワイズ全単射を使用します(これは全単射変換です)。

これで、必要なコピュラを定義できます。 (バイジェクターとしてエンコードされた)ターゲット周辺分布のリストが与えられると、コピュラを使用し、指定された周辺分布を持つ新しい分布を簡単に構築できます。

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

png

png

png

最後に、必要な限界が実際に得られることを確認しましょう。

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

png

結論

そして、そこに行きます!私たちは、私たちが使用してガウスコピュラを構築することができることを実証してきましたBijector APIを。

より一般的には、 Bijector APIを使用してBijectorイジェクターを作成し、それらをディストリビューションで構成すると、柔軟なモデリングのための豊富なディストリビューションファミリーを作成できます。