【R】２次元正規分布を観察 - ケィオスの時系列解析メモランダム

2次元正規分布は、2つの連続変数が同時に正規分布に従う場合の確率分布を表すものです。これは、1変数の正規分布を2次元空間に拡張したものであり、2変数間の相関関係も含めて記述しています。

たとえば、身長と体重、温度と湿度、脳活動の2つの信号など、同時に測定される連続データの関係性を捉える際に広く用いられます。

基礎

2次元正規分布では、以下の情報をもとに、データの広がり方や偏り方を記述します：

平均ベクトル（それぞれの変数の平均値）

$\displaystyle{ \boldsymbol{\mu} = \begin{pmatrix} \mu_x \\ \mu_y \end{pmatrix} }$

共分散行列（それぞれの分散と、変数間の相関）

$\displaystyle{ \Sigma = \begin{pmatrix} \sigma_x^2 & \rho \sigma_x \sigma_y \\ \rho \sigma_x \sigma_y & \sigma_y^2 \end{pmatrix} }$

確率密度関数（PDF） 2次元空間上の点 \$(x, y)\$ に対する確率密度は以下で表されます：

$\displaystyle{ f(x, y) = \frac{1}{2\pi \sigma_x \sigma_y \sqrt{1 - \rho^2}} \exp\left[ -\frac{1}{2(1 - \rho^2)} \left( \frac{(x - \mu_x)^2}{\sigma_x^2} + \frac{(y - \mu_y)^2}{\sigma_y^2} - \frac{2\rho(x - \mu_x)(y - \mu_y)}{\sigma_x \sigma_y} \right) \right] }$

Rスクリプトで分布を描画

２つの変数の相関係数rhoが与えられたときの､２次元正規分布を描くRスクリプトの例です。

# ----------------------------
# 2次元正規分布の立体プロット（4方向）
# ----------------------------

# 平均と共分散
mu_x <- 0
mu_y <- 0
sigma_x <- 1
sigma_y <- 1
rho <- 0  # 相関係数

# グリッド生成
x <- seq(-3, 3, length = 100)
y <- seq(-3, 3, length = 100)
z <- matrix(0, nrow = length(x), ncol = length(y))

# 定数項
norm_const <- 1 / (2 * pi * sigma_x * sigma_y * sqrt(1 - rho^2))

# 密度関数の計算
for (i in 1:length(x)) {
  for (j in 1:length(y)) {
    dx <- x[i] - mu_x
    dy <- y[j] - mu_y
    exponent <- -1 / (2 * (1 - rho^2)) * (
      (dx^2) / sigma_x^2 +
      (dy^2) / sigma_y^2 -
      2 * rho * dx * dy / (sigma_x * sigma_y)
    )
    z[i, j] <- norm_const * exp(exponent)
  }
}

# 画面を2x2に分割
par(mfrow = c(2, 2))

# 4方向から描画
persp(x, y, z,
      theta = 30, phi = 30,
      expand = 0.5,
      col = "lightblue",
      xlab = "X", ylab = "Y", zlab = "Density",
      main = "View 1: theta=30, phi=30")

persp(x, y, z,
      theta = 60, phi = 20,
      expand = 0.5,
      col = "lightgreen",
      xlab = "X", ylab = "Y", zlab = "Density",
      main = "View 2: theta=60, phi=20")

persp(x, y, z,
      theta = 120, phi = 40,
      expand = 0.5,
      col = "lightcoral",
      xlab = "X", ylab = "Y", zlab = "Density",
      main = "View 3: theta=120, phi=40")

persp(x, y, z,
      theta = 300, phi = 25,
      expand = 0.5,
      col = "lightyellow",
      xlab = "X", ylab = "Y", zlab = "Density",
      main = "View 4: theta=300, phi=25")