ไพรเมอร์ Copulas

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(U1,U2,....,Un) ดังกล่าวว่า marginalizing ให้ UiUniform(0,1)

Copulas นั้นน่าสนใจเพราะเราสามารถใช้พวกมันเพื่อสร้างการแจกแจงแบบหลายตัวแปรด้วยส่วนเพิ่มตามอำเภอใจ นี่คือสูตร:

  • ใช้ เป็นส่วนประกอบสำคัญน่าจะเป็นเปลี่ยน ผลัดพล RV อย่างต่อเนื่อง X เข้าไปในเครื่องแบบหนึ่ง FX(X)ที่ FX เป็น CDF ของ X
  • ได้รับการเชื่อม (พูด bivariate) C(U,V)เรามีที่ U และ V มีการกระจายร่อแร่เครื่องแบบ
  • ตอนนี้ได้รับ RV ของเราที่น่าสนใจ X,Yสร้างใหม่กระจาย C(X,Y)=C(FX(X),FY(Y))มาร์จินสำหรับ X และ Y จะเป็นคนที่เราต้องการ

ระยะขอบไม่แปรผัน ดังนั้นจึงอาจวัดและ/หรือแบบจำลองได้ง่ายกว่า คอปูลาช่วยให้เริ่มต้นจากส่วนขอบ แต่ยังบรรลุความสัมพันธ์ตามอำเภอใจระหว่างมิติ

Gaussian Copula

เพื่อแสดงให้เห็นวิธีการสร้าง copulas ให้พิจารณากรณีของการจับการพึ่งพาตามสหสัมพันธ์แบบเกาส์เซียนพหุตัวแปร Gaussian เชื่อมเป็นหนึ่งที่ได้รับจาก C(u1,u2,...un)=ΦΣ(Φ1(u1),Φ1(u2),...Φ1(un)) ที่ ΦΣ หมายถึง CDF ของ MultivariateNormal ที่มีความแปรปรวน Σ และค่าเฉลี่ย 0 และ Φ1 เป็น CDF ผกผันสำหรับมาตรฐานปกติ

การใช้ CDF ผกผันปกติจะบิดเบือนมิติที่สม่ำเสมอเพื่อกระจายตามปกติ การใช้ CDF ของค่าปกติแบบพหุตัวแปรแล้วบีบการแจกแจงให้เป็นแบบเดียวกันเล็กน้อยและมีความสัมพันธ์แบบเกาส์เซียน

ดังนั้นสิ่งที่เราได้รับคือว่า Gaussian เชื่อมคือการกระจายมากกว่าหน่วย hypercube [0,1]n กับมาร์จินเครื่องแบบ

กำหนดให้เป็นเช่นนั้นเสียนเชื่อมสามารถดำเนินการกับ tfd.TransformedDistribution และเหมาะสม Bijector นั่นก็คือเรามีการปรับเปลี่ยน MultivariateNormal, ผ่านการใช้งานของการกระจายปกติของผกผัน CDF ที่ดำเนินการโดย tfb.NormalCDF bijector

ด้านล่างเราใช้ Gaussian เชื่อมกับสมมติฐานการลดความซับซ้อน: ที่แปรปรวนจะแปรโดยปัจจัย 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);

png

อย่างไรก็ตาม พลังจากแบบจำลองดังกล่าวกำลังใช้ Probability Integral Transform เพื่อใช้ copula กับ RVs ที่กำหนดเอง ด้วยวิธีนี้ เราสามารถระบุ Marginals ตามอำเภอใจ และใช้ copula เพื่อต่อเข้าด้วยกัน

เราเริ่มต้นด้วยแบบจำลอง:

XKumaraswamy(a,b)YGumbel(μ,β)

และใช้เชื่อมที่จะได้รับ bivariate RV Zซึ่งมีมาร์จิน Kumaraswamy และ กัมเบล

เราจะเริ่มต้นด้วยการวางแผนการกระจายผลิตภัณฑ์ที่สร้างโดย RV สองตัวนี้ นี่เป็นเพียงเพื่อใช้เป็นจุดเปรียบเทียบเมื่อเราใช้ Copula

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

การกระจายร่วมกับขอบต่างกัน

ตอนนี้เราใช้คอปูลาเกาส์เซียนเพื่อรวมการแจกแจงเข้าด้วยกัน และพลอตมัน อีกครั้งเครื่องมือของทางเลือกของเราคือการ TransformedDistribution ใช้ที่เหมาะสม Bijector ที่จะได้รับมาร์จินที่เลือก

โดยเฉพาะเราใช้ Blockwise bijector ซึ่งใช้ bijectors แตกต่างกันในส่วนต่าง ๆ ของเวกเตอร์ (ซึ่งยังคงเป็นการเปลี่ยนแปลง bijective)

ตอนนี้เราสามารถกำหนด Copula ที่เราต้องการได้ จากรายชื่อเป้าหมายเป้าหมาย (เข้ารหัสเป็น bijectors) เราสามารถสร้างการกระจายใหม่ที่ใช้ copula และมีระยะขอบที่ระบุได้อย่างง่ายดาย

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

สุดท้ายนี้ เรามาลองใช้ Gaussian Copula กันดีกว่า เราจะใช้ Cholesky ของ [10ρ(1ρ2)]ซึ่งจะสอดคล้องกับความแปรปรวนที่ 1 และความสัมพันธ์ ρ สำหรับหลายตัวแปรปกติ

เราจะพิจารณาบางกรณี:

# 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

บทสรุป

และเราไปกันเลย! เราได้แสดงให้เห็นว่าเราสามารถสร้างแบบเกาส์ Copulas ใช้ Bijector API

โดยทั่วไปการเขียน bijectors ใช้ Bijector API และเขียนพวกเขาด้วยการจัดจำหน่ายที่สามารถสร้างครอบครัวที่อุดมไปด้วยการกระจายสำหรับการสร้างแบบจำลองที่มีความยืดหยุ่น