1 固有値による分類.

2次正方行列で表される一次変換について、固有値と固有ベクトルに注目して調べてきました。
概略は次のとおりです。

異なる2つの実固有値を持つ場合：
異なる2つの方向に拡大縮小する変換です。それぞれの方向に固有値の倍率で伸縮します。
異なる2つの複素固有値を持つ場合：
固有値が実数でない場合、変換は原点を中心とする回転を含みます。
重解の実固有値を持つ場合：
一次独立な固有ベクトルが2つのとき、相似変換 \(kE\) となり、
一次独立な固有ベクトルが1つのとき、せん断変換を含みます。

このように固有値と一次変換の関係について全体を概観しておくことは、予想を立てる上で有効です。
しかし、行列がどのような固有値と固有ベクトルを持つか、行列の見た目では分かりません。

代表的な行列について固有値を整理しておけば、見通し良く考えることができるでしょう。

1.1 実固有値を持つタイプ.

\(\begin{pmatrix} a & 0 \\ 0 & a \end{pmatrix}\)：原点を中心とした a 倍の相似変換.
固有値は \(a\) を重解に持ちます。
\(\begin{pmatrix} a & b \\ b & a \end{pmatrix}\)：対称行列タイプ１.
２つの固有値は \(a\pm b\)、固有ベクトルは　\(y=\pm x\) 方向です。
\(\begin{pmatrix} a & b \\ b & d \end{pmatrix}\)：対称行列タイプ２.
２つの異なる実固有値を持ち、固有ベクトルは直交します。

これらは、計算が簡明で分かりやすいので、重要な例として説明してきました。

対称行列は直交する固有ベクトルを持ちますが、一般の行列ではどうでしょうか。

固有ベクトルが直交しない場合でも、「6.3 固有ベクトル基底」で説明したように、

　　\(\overrightarrow{x}=u~\overrightarrow{p}+v~\overrightarrow{q}\) に対し、

　　　　　\(A\overrightarrow{x}=ku~\overrightarrow{p}+lv~\overrightarrow{q}\)

固有ベクトルを基底とした場合、一次変換によって座標はそれぞれの固有値倍になります。

1.1.1 三角行列による像.

対角成分が固有値、対角成分の積が行列式となっており、計算が容易です。
三角行列は計算の簡便性から、行列の分解に用いられます。

＜上三角行列＞

\(\begin{pmatrix} a & b \\ 0 & d \end{pmatrix}\)を上三角行列（Upper triangular matrix）と呼びます。

２つの固有値は \(a\) と \(d\)、固有ベクトルは x軸方向と \((b,d-a)\) です。

不動直線は、x軸と \(y=\dfrac{d-a}{b}x\)　です。

\(A=\begin{pmatrix} 1 & b \\ 0 & 1.5 \end{pmatrix}\) について、\(b=0.5\) の場合を考えます。

b <- 0.5
A <- matrix(c(1, 0, b, 1.5), 2, 2)
A

##      [,1] [,2]
## [1,]    1  0.5
## [2,]    0  1.5

eigen <- eigen(A)
v <- round(eigen$vectors,2)
v

##      [,1] [,2]
## [1,] 0.71    1
## [2,] 0.71    0

e <- round(eigen$values,2)
e

## [1] 1.5 1.0

固有値は行列の対角成分であり、x軸を不動直線に持つことが確かめられます。

＜固有ベクトルで張られる平行四辺形の像＞
平行四辺形内の点は、

　　　\(x_1+x_2=s\cdot v[,1]+t\cdot v[,2]\)　で表されます。

s <- runif(1000, min = 0, max = 1)
x_1 <- sapply(s, function(scalar) scalar * v[,1])
t <- runif(1000, min = 0, max = 1)
x_2 <- sapply(t, function(scalar) scalar * v[,2])
x <- x_1+x_2
data <- data.frame(x = x[1,], y = x[2,])
par(mfrow=c(1,2))
d_new <- A %*% t(data)
d_new <- data.frame(X=d_new[1,],Y=d_new[2,])
plot(data,asp=1,col="lightblue",cex=0.2,xlim=c(-2,2),ylim=c(-2,2))
abline(a = 0, b=v[2,1]/v[1,1],col = "lightblue", lty = 1)
abline(a = 0, b=v[2,2]/v[1,2],col = "lightblue", lty = 1)
plot(d_new,col="pink",cex=0.2,xlim=c(-2,2),ylim=c(-2,2),asp=1)
abline(a = 0, b=v[2,1]/v[1,1],col = "lightblue", lty = 1)
abline(a = 0, b=v[2,2]/v[1,2],col = "lightblue", lty = 1)
arrows(x0=0, y0=0, x1=v[1,1],  y1=v[2,1], code=2, col="blue")
arrows(x0=0, y0=0, x1=v[1,2],  y1=v[2,2], code=2, col="red")

不動直線は x軸と直線 \(y=x\) となっています。

＜下三角行列＞

\(\begin{pmatrix} a & 0 \\ b & d \end{pmatrix}\)を下三角行列 (Lower Triangular Matrix) と呼びます。

２つの固有値は \(a\)と \(d\)、固有ベクトルは \((d-a,-b)\) と y 軸方向です。
不動直線は、y軸と \(y=\dfrac{-b}{d-a}x\)　です。

\(B=\begin{pmatrix} 1 & 0 \\ b & 1.5 \end{pmatrix}\) について、\(b=0.5\) の場合を考えます。

b <- 0.5
B <- matrix(c(1, b, 0, 1.5), 2, 2)
B

##      [,1] [,2]
## [1,]  1.0  0.0
## [2,]  0.5  1.5

eigen <- eigen(B)
v <- round(eigen$vectors,1)
v

##      [,1] [,2]
## [1,]    0  0.7
## [2,]    1 -0.7

e <- round(eigen$values,1)
e

## [1] 1.5 1.0

固有値は行列の対角成分であり、y軸を不動直線に持つことが確かめられます。

＜固有ベクトルで張られる平行四辺形の像＞
平行四辺形内の点は、

　　　\(x_1+x_2=s\cdot v[,1]+t\cdot v[,2]\)　で表されます。

s <- runif(1000, min = 0, max = 1)
x_1 <- sapply(s, function(scalar) scalar * v[,1])
t <- runif(1000, min = 0, max = 1)
x_2 <- sapply(t, function(scalar) scalar * v[,2])
x <- x_1+x_2
data <- data.frame(x = x[1,], y = x[2,])
par(mfrow=c(1,2))
d_new <- B %*% t(data)
d_new <- data.frame(X=d_new[1,],Y=d_new[2,])
plot(data,asp=1,col="lightblue",cex=0.2,xlim=c(-2,2),ylim=c(-2,2))
abline(v = 0,col = "lightblue", lty = 1)
abline(a = 0, b=v[2,2]/v[1,2],col = "lightblue", lty = 1)
plot(d_new,col="pink",cex=0.2,xlim=c(-2,2),ylim=c(-2,2),asp=1)
abline(v = 0,col = "lightblue", lty = 1)
abline(a = 0, b=v[2,2]/v[1,2],col = "lightblue", lty = 1)
arrows(x0=0, y0=0, x1=v[1,1],  y1=v[2,1], code=2, col="blue")
arrows(x0=0, y0=0, x1=v[1,2],  y1=v[2,2], code=2, col="red")

不動直線は x軸と直線 \(y=-x\) となっています。

［注意］
　上三角行列の固有ベクトルと下三角行列の固有ベクトルは互いに直交しています。

また、上三角行列 A と下三角行列 B が互いに転置の関係にあるとき、
その積 \(C=BA\) はタイプ1の対称行列になります。

C <- B %*% A
C

##      [,1] [,2]
## [1,]  1.0  0.5
## [2,]  0.5  2.5

このことは、対称行列が下三角行列と上三角行列の積で表される可能性があることを示しています。

＜LU分解＞
　対称行列 \(\Sigma\) に対し、

下三角行列 \(A\) が存在して、\(\Sigma=AA'\)
上三角行列 \(A\) が存在して、\(\Sigma=A'A\)

このどちらの表現をとる場合もあります。
これを、対称行列のLU分解といいます。
（LはLowerの頭文字、UはUpperの頭文字です）
計算には後述のコレスキー分解などが用いられます。

これは、対称行列が三角行列から構成できることを意味します。

1.1.2 異なる固有値を持つ行列.

異なる固有値に対する固有ベクトルは一次独立です。
しかし、対称行列でない場合、固有ベクトルは直交するとは限りません。
この場合でも、一次変換は固有ベクトル方向の固有値倍の拡大縮小となっています。

実際にシミュレーションで確かめましょう。

\(\begin{pmatrix} 1 & b \\ b+1 & 1 \end{pmatrix}\)：固有ベクトルは直交せず斜交します。

斜交する固有ベクトルを基底とする座標を斜交座標といいます。

b <- 0.2
A <- matrix(c(1, b, b+1, 1), 2, 2)
A

##      [,1] [,2]
## [1,]  1.0  1.2
## [2,]  0.2  1.0

eigen <- eigen(A)
v <- eigen$vectors
v

##           [,1]       [,2]
## [1,] 0.9258201 -0.9258201
## [2,] 0.3779645  0.3779645

e <- eigen$values
e

## [1] 1.4898979 0.5101021

＜固有ベクトルで張られる平行四辺形の像＞.

s <- runif(1000, min = 0, max = 1)
x_1 <- sapply(s, function(scalar) scalar * v[,1])
t <- runif(1000, min = 0, max = 1)
x_2 <- sapply(t, function(scalar) scalar * v[,2])
x <- x_1+x_2
data <- data.frame(x = x[1,], y = x[2,])
par(mfrow=c(1,2))
d_new <- A %*% t(data)
d_new <- data.frame(X=d_new[1,],Y=d_new[2,])
plot(data,asp=1,col="lightblue",cex=0.2,xlim=c(-2,2),ylim=c(-2,2))
abline(a = 0, b=v[2,1]/v[1,1],col = "lightblue", lty = 1)
abline(a = 0, b=v[2,2]/v[1,2],col = "lightblue", lty = 1)
plot(d_new,col="pink",cex=0.2,xlim=c(-2,2),ylim=c(-2,2),asp=1)
abline(a = 0, b=v[2,1]/v[1,1],col = "lightblue", lty = 1)
abline(a = 0, b=v[2,2]/v[1,2],col = "lightblue", lty = 1)
arrows(x0=0, y0=0, x1=v[1,1],  y1=v[2,1], code=2, col="blue")
arrows(x0=0, y0=0, x1=v[1,2],  y1=v[2,2], code=2, col="red")

第１固有ベクトルの方向へ 1.5倍、第２個湯ベクトルの方向へ 0.5倍の変換になっていることが分かります。

このように、

２つの異なる実固有値を持つ場合、固有ベクトルは一次独立であり、
固有ベクトルを基底とすれば、一次変換は固有ベクトル方向の固有値倍の拡大縮小となります。

しかし、固有値が１つで、一次独立な固有ベクトルが存在しない場合はどうなるでしょうか。

1.1.3 固有値が重解の場合.

\(\begin{pmatrix} a & b \\ 0 & a \end{pmatrix}\)：上三角行列.
固有値は \(a\) を重解に持ちます。固有ベクトルは x軸方向のみです。

\(\begin{pmatrix} a & b \\ 0 & a \end{pmatrix} =a\begin{pmatrix} 1 & b/a \\ 0 & 1 \end{pmatrix}\)　なので、

まず　\(\begin{pmatrix} 1 & k \\ 0 & 1 \end{pmatrix}\)　タイプの変換を行い、（\(k=b/a\)）
次に a 倍の拡大変換をすることになります。

そこでまず、

　　　\(\begin{pmatrix} X \\ Y \end{pmatrix} =\begin{pmatrix} 1 & k \\ 0 & 1 \end{pmatrix} \begin{pmatrix} x \\ y \end{pmatrix} =\begin{pmatrix} x+ ky\\ y \end{pmatrix}\)　を考えます。

この変換では、

y座標は変化させず.
x座標はyに比例して \(ky\)だけ大きくなります。

点 \((x,y)\) は横に移動して点 \((X,Y)\) になります。

\((X,Y)\) は　\(Y=\dfrac{1}{k}(X-x)\)　を満たすので、
\((x,0)\) を通り、傾き \(\dfrac{1}{k}=\dfrac{a}{b}\) の直線上にあります。

　
＜例：上三角行列＞

　　\(\begin{pmatrix} 2 & 1 \\ 0 & 2 \end{pmatrix} =2\begin{pmatrix} 1 & 1/2 \\ 0 & 1 \end{pmatrix}\)　の場合を考えます。

A <- matrix(c(1, 0, 1/2, 1), 2, 2)
A

##      [,1] [,2]
## [1,]    1  0.5
## [2,]    0  1.0

eigen <- eigen(A)
v <- eigen$vectors
round(v,1)

##      [,1] [,2]
## [1,]    1   -1
## [2,]    0    0

e <- eigen$values
e

## [1] 1 1

\((0,0)\) \((0,1)\)を両端とする線分（青）を移動します。

x <- 0
y <- runif(1000, min = 0, max = 1)
data <- data.frame(x = x, y = y)
par(mfrow=c(1,2))
d_new <- A %*% t(data)
d_new <- data.frame(X=d_new[1,],Y=d_new[2,])
plot(data,asp=1,col="lightblue",cex=0.2,xlim=c(-1,1),ylim=c(-1,1))
abline(a = 0, b=v[2,1]/v[1,1],col = "lightblue", lty = 1)
abline(a = 0, b=v[2,2]/v[1,2],col = "lightblue", lty = 1)
plot(d_new,col="pink",cex=0.2,xlim=c(-1,1),ylim=c(-1,1),asp=1)
abline(a = 0, b=v[2,1]/v[1,1],col = "lightblue", lty = 1)

移動先は直線 \(Y=2X\)　の　\(0<Y<1\) に対応する線分（赤）です。

右図の直線の方程式は、
　　　\(Y=\dfrac{1}{k}(X-x)\)、即ち　\(Y=\dfrac{a}{b}(X-x)\) なので、

直線 \(x=0\) 上の点は、\(y=2x\) に移ります。

\(x\) を0 から 1 まで動かせば、赤の斜め線は横に平行に移動していきます。

単位正方形を変換してみましょう。

x <- runif(1000, min = 0, max = 1)
y <- runif(1000, min = 0, max = 1)
data <- data.frame(x = x, y = y)
par(mfrow=c(1,2))
d_new <- A %*% t(data)
d_new <- data.frame(X=d_new[1,],Y=d_new[2,])
plot(data,asp=1,col="lightblue",cex=0.2,xlim=c(-1.5,1.5),ylim=c(-1.5,1.5))
abline(a = 0, b=v[2,1]/v[1,1],col = "lightblue", lty = 1)
abline(a = 0, b=v[2,2]/v[1,2],col = "lightblue", lty = 1)
plot(d_new,col="pink",cex=0.2,xlim=c(-1.5,1.5),ylim=c(-1.5,1.5),asp=1)
abline(a = 0, b=v[2,1]/v[1,1],col = "lightblue", lty = 1)

このような変換をせん断といいます。 .
あとは、2倍の拡大を行えばよいことになります。

つまり、

重解の実固有値を持ち、一次独立な固有ベクトルが１つである場合、
一次変換は、せん断と拡大縮小の合成です。

\(\begin{pmatrix} 2 & 1 \\ 0 & 2 \end{pmatrix}\) で直接変換したものを以下に示します。

A <- matrix(c(2, 0, 1, 2), 2, 2)
A

##      [,1] [,2]
## [1,]    2    1
## [2,]    0    2

eigen <- eigen(A)
v <- eigen$vectors
e <- eigen$values
x <- runif(1000, min = 0, max = 1)
y <- runif(1000, min = 0, max = 1)
data <- data.frame(x = x, y = y)
par(mfrow=c(1,2))
d_new <- A %*% t(data)
d_new <- data.frame(X=d_new[1,],Y=d_new[2,])
plot(data,asp=1,col="lightblue",cex=0.2,xlim=c(-2,2),ylim=c(-2,2))
abline(a = 0, b=v[2,1]/v[1,1],col = "lightblue", lty = 1)
abline(a = 0, b=v[2,2]/v[1,2],col = "lightblue", lty = 1)
plot(d_new,col="pink",cex=0.2,xlim=c(-1,3),ylim=c(-2,2),asp=1)
abline(a = 0, b=v[2,1]/v[1,1],col = "lightblue", lty = 1)
abline(a = 0, b=v[2,2]/v[1,2],col = "lightblue", lty = 1)
arrows(x0=0, y0=0, x1=v[1,1],  y1=v[2,1], code=2, col="blue")
arrows(x0=0, y0=0, x1=v[1,2],  y1=v[2,2], code=2, col="red")

せん断による変換を更に2倍に拡大した図になっています。

＜例：下三角行列＞

\(\begin{pmatrix} 2 & 0 \\ 1 & 2 \end{pmatrix} =2\begin{pmatrix} 1 & 0 \\ 1/2 & 1 \end{pmatrix}\).

下三角行列の場合、せん断は y 軸方向になりますが、あとは同様です。

A <- matrix(c(2, 1, 0, 2), 2, 2)
A

##      [,1] [,2]
## [1,]    2    0
## [2,]    1    2

eigen <- eigen(A)
v <- eigen$vectors
round(v,1)

##      [,1] [,2]
## [1,]    0    0
## [2,]    1   -1

e <- eigen$values
e

## [1] 2 2

x <- runif(1000, min = 0, max = 1)
y <- runif(1000, min = 0, max = 1)
data <- data.frame(x = x, y = y)
par(mfrow=c(1,2))
d_new <- A %*% t(data)
d_new <- data.frame(X=d_new[1,],Y=d_new[2,])
plot(data,asp=1,col="lightblue",cex=0.2,xlim=c(-2,2),ylim=c(-2,2))
abline(v = 0,col = "lightblue", lty = 1)
abline(a = 0, b=v[2,2]/v[1,2],col = "lightblue", lty = 1)
plot(d_new,col="pink",cex=0.2,xlim=c(-2,2),ylim=c(-2,2),asp=1)
abline(v = 0,col = "lightblue", lty = 1)
abline(a = 0, b=v[2,2]/v[1,2],col = "lightblue", lty = 1)
arrows(x0=0, y0=0, x1=v[1,1],  y1=v[2,1], code=2, col="blue")
arrows(x0=0, y0=0, x1=v[1,2],  y1=v[2,2], code=2, col="red")

y軸方向へのせん断と2倍の拡大になっています。

1.2 複素固有値を持つタイプ.

\(\begin{pmatrix} a & -b \\ b & a \end{pmatrix}\)：回転と拡大.
異なる2つの複素固有値 \(a\pm bi\) を持ち、固有ベクトルを持ちません。

\(r=\sqrt{a^2+b^2}\) 、\(\theta\) は \((a,b)\) の偏角とすると、

　　　\(\begin{pmatrix} a & -b \\ b & a \end{pmatrix}=r\begin{pmatrix} a/r & -b/r \\ b/r & a/r \end{pmatrix} =r\begin{pmatrix} \cos\theta & -\sin\theta \\ \sin\theta & \cos\theta \end{pmatrix}\)　なので、.

原点の周りの \(\theta\)回転と\(r\)倍の拡大縮小を行う変換になります。

つまり、複素数の偏角の回転と複素数の大きさを拡大倍率とする相似変換の合成となります。

＜例：45度回転＞

b <- 1
A <- matrix(c(1, b, -b, 1), 2, 2)
A

##      [,1] [,2]
## [1,]    1   -1
## [2,]    1    1

eigen <- eigen(A)
v <- eigen$vectors
round(v,1)

##          [,1]     [,2]
## [1,] 0.7+0.0i 0.7+0.0i
## [2,] 0.0-0.7i 0.0+0.7i

e <- eigen$values
round(e,1)

## [1] 1+1i 1-1i

固有値、固有ベクトルとも複素数として算出されています。

単位正方形を変換してみます。

x <- runif(1000, min = 0, max = 1)
y <- runif(1000, min = 0, max = 1)
data <- data.frame(x = x, y = y)
d_new <- A %*% t(data)
d_new <- data.frame(X=d_new[1,],Y=d_new[2,])
par(mfrow=c(1,2))
plot(data,asp=1,col="lightblue",cex=0.1,xlim=c(-2,2),ylim=c(-2,2))
abline(h = 0, v=0,col = "lightblue", lty = 1)
plot(d_new,col="pink",cex=0.1,xlim=c(-2,2),ylim=c(-2,2),asp=1)
abline(h=0,v = 0, col = "lightblue")

原点を中心に 45°回転し、\(\sqrt{2}\fallingdotseq 1.4\)倍に拡大する変換になっています。

［注意］
次の行列は似ていますが回転ではありません。

\(\begin{pmatrix} a & b \\ b & -a \end{pmatrix}\)：直線に関する対称移動.
２つの固有値 \(\pm \sqrt{a^2+b^2}\) を持ち、固有ベクトルは直交します。

回転・拡大行列と似ていますが、これは対称行列です。

このように、簡明なタイプをシミュレーションで確かめておくことは、一次変換を理解する助けになります。

2 二変量正規分布.

\(X\sim N(\mu,\sigma^2)\)　に対し、\(Z=\dfrac{X-\mu}{\sigma}\sim N(0,1)\).

この操作を標準化といいます。この操作を逆に行うと、

\(Z\sim N(0,1)\)　に対し、\(X=\sigma Z +\mu\sim N(\mu,\sigma^2)\).

\(Z\) に一次関数を適用することで、一般の正規分布を構成できます。

同様の方法で 2変量の正規分布を構成します。

2.1 標準正規乱数.

１変量確率変数 x,y が独立に標準正規分布に従うとき、２変量確率変数 (x,y) は２変量標準正規分布に従います。

\(x,y \sim N(0,1):i.i.d.\)　のとき　\((x,y) \sim N_2(0,I_2)\).
x,y の分散共分散行列は　\(I_2=\begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix}\)　になります。

\(X_1,X_2\sim N(0,1)\:\:\:i.i.d.\)のとき、同時分布の確率密度関数は、\((X_1,X_2)\)の独立性から各々の確率密度関数の積となるので、

　　　\(\displaystyle f(x_1,x_2)= \dfrac{1}{\sqrt{2\pi}}\exp\left(-\dfrac{1}{2}x_1^2\right)\cdot \dfrac{1}{\sqrt{2\pi}}\exp\left(-\dfrac{1}{2}x_2^2\right)\).

この乗法定理は、\(X_1,X_2\)が独立のとき成立します。
式をまとめると、
　　　\(\displaystyle f(x_1,x_2)= \dfrac{1}{\sqrt{2\pi}}\exp\left[-\dfrac{1}{2}(x_1^2+x_2^2)\right]\)

2変量確率密度関数は、同心円を等確率曲線とするベル状の局面になります。

library(rgl)
options(rgl.printRglwidget = TRUE) 
close3d()
x1 <- rnorm(1000, 0,1)
x2 <- rnorm(1000, 0,1)
z <- dnorm(x1,0,1)*dnorm(x2,0,1)
plot3d(x1, x2, z, type = "s", col = "pink", size = 0.3,xlim=c(-3,3),ylim=c(-3,3),zlim=c(0,0.2))

ひょっこり統計「高校生の統計学」の「高校生の線形代数」では、ブラウザ上で回転やズームが可能です。

回転してz方向から俯瞰すると、2変量標準正規分布の分布を見ることができます。

点(x,y)を平面上に表示すると、次のような円上の分布になります。

n <- 10000
x <- rnorm(n,0,1)
y <- rnorm(n,0,1)
plot(x,y,asp=1,col="lightblue",cex=0.1)

data <- data.frame(X=x,Y=y)
round(cov(data),1)

##   X Y
## X 1 0
## Y 0 1

2変量確率変数では、分散共分散行列 cov(data) が 2変量データの散らばり具合を示す指標であり、2変量の関連の度合いを表す共分散を含みます。

2変量標準正規分布の分散共分散行列は単位行列になります。

2.2 分散共分散行列

2変量確率変数 \(X=\begin{pmatrix}X_1\\ X_2\end{pmatrix}\)　の期待値を \(m=\begin{pmatrix}m_1\\ m_2\end{pmatrix}\)　としたとき、

　　\(Cov\left(X\right)=E\left((X-m)(X-m)'\right)\) を分散共分散行列といます。

成分で計算すると、

　\(\begin{split}(X-m)(X-m)'&=\begin{pmatrix} X_1-m_1 \\ X_2-m_2 \end{pmatrix} \begin{pmatrix} X_1-m_1 & X_2-m_2 \end{pmatrix}\\&= \begin{pmatrix} (X_1-m_1)^2 & (X_1-m_1)(X_2-m_2) \\ (X_1-m_1)(X_2-m_2) & (X_2-m_2)^2 \end{pmatrix} \end{split}\)

\(\therefore E\left((X-m)(X-m)'\right)=\begin{pmatrix} V(X_1) & Cov(X_1,X_2) \\ Cov(X_1,X_2) & V(X_2) \end{pmatrix}\)

対角成分は分散、対角でない成分は共分散になります。
２変量の確率変数では分散のほかに変量相互の関連を示す共分散が重要です。

2変量での分散共分散行列は、１変量での分散に相当します。

１変量では、
　　\(Z\sim N(0,1)\) に対し　\(X=\sigma Z +\mu\sim N(\mu,\sigma^2)\)
　　標準正規分布に標準偏差を乗じることで確率変数を構成できました。

２変量では、
　　\(Z\sim N_2(0,I_2)\) に対し　\(X=AZ+m\) の分散共分散行列は、

　　\(\begin{split} 　　Cov(X)&=E\left((X-m)(X-m)'\right)=E(AZZ'A')\\ 　　&=A~E(ZZ')~A'=A~I_2~A'=AA' 　　\end{split}\)　　

　　\(X\sim N_2(m,\Sigma)\) となるには、
　　　　\(\Sigma=AA'\) を満たす\(A\) を用いて　\(X=AZ+m\) とする必要があります。.

\(AA'=\Sigma\) を満たす\(A\)は、行列 \(\Sigma\) の平方根に相当します。

2.3 行列の対角化

Xの分散共分散行列 \(\Sigma\) に対し、\(AA'=\Sigma\) を満たす\(A\)を求めます。

\(\Sigma\) を固有値と固有ベクトルを用いて対角化することで、 \(A\)を求めることができます。

\(\Sigma\) が２つの異なる固有値 \(\lambda_1,\lambda_2\) を持つとし、
標準化された固有ベクトルを \(\begin{pmatrix} x_{1} \\ y_{1} \end{pmatrix},\begin{pmatrix} x_{2} \\ y_{2} \end{pmatrix}\) とすると、

　　\(\Sigma\begin{pmatrix} x_{1} & x_{2} \\ y_{1} & y_{2} \end{pmatrix} =\lambda _{1}\begin{pmatrix} x_{1} & 0 \\ y_{1} & 0 \end{pmatrix}+\lambda _{2}\begin{pmatrix} 0 & x_{2} \\ 0 & y_{2} \end{pmatrix}\)

\(\begin{pmatrix} x_{1} & x_{2} \\ y_{1} & y_{2} \end{pmatrix}^{-1} \left[ \lambda _{1}\begin{pmatrix} x_{1} & 0 \\ y_{1} & 0 \end{pmatrix}+\lambda _{2}\begin{pmatrix} 0 & x_{2} \\ 0 & y_{2} \end{pmatrix}\right] =\lambda _{1}\begin{pmatrix} 1&0 \\ 0&0 \end{pmatrix}+\lambda _{2}\begin{pmatrix} 0&0 \\ 0&1 \end{pmatrix}\)

よって、　\(\begin{pmatrix} x_{1} & x_{2} \\ y_{1} & y_{2} \end{pmatrix}^{-1} \Sigma\begin{pmatrix} x_{1} & x_{2} \\ y_{1} & y_{2} \end{pmatrix} =\begin{pmatrix} \lambda_{1} & 0 \\ 0 & \lambda_{2} \end{pmatrix}\)

標準化された固有ベクトルを列とする行列を \(P=\begin{pmatrix} x_{1} &x_{2} \\ y_{1} &y_{2} \end{pmatrix}\)　とすると、
　　　\(P^{-1}\Sigma P =\begin{pmatrix} \lambda_{1} & 0 \\ 0 & \lambda_{2} \end{pmatrix}\) 　（対角成分は固有値）.

これを行列の対角化といいます。

［注意］.
固有ベクトルが標準化されていなくても、対角化の式は成立します。
分散共分散行列は対称行列なので、異なる固有値に対する固有ベクトルは直交し、標準化しておくと P は直交行列になります。

直交行列は、各列が直交する単位ベクトルで構成され、転置行列が逆行列になります。
このため計算が簡潔になります。

以下、具体例で見ていきましょう。

2.4 行列の平方根.

簡単のため、対称行列タイプ１の場合を考えます。

\(\Sigma= \begin{pmatrix} a & b \\ b & a \end{pmatrix}\) は、\(b \neq 0\) のとき、２つの異なる固有値を持ちます。
実際、
固有方程式　\(\lambda ^{2}-\left( 2a\right) \lambda +\left( a^2-b^2\right) =0\) の解は　\(\lambda=a \pm b\)

\(\Sigma= \begin{pmatrix} a & b \\ b & a \end{pmatrix} =\begin{pmatrix} V(X_1) & Cov(X_1,X_2) \\ Cov(X_1,X_2) & V(X_2) \end{pmatrix}\)　であり、

相関係数　\(r=\dfrac{Cov(X_1,X_2)}{\sqrt{V(X_1)}\sqrt{V(X_2)}}=\dfrac{b}{a}<1\) より　\(a>b\).

よって、\(\lambda_1=a+b>0,~~\lambda_2=a-b>0\)

タイプ１の分散共分散行列は２つの異なる正の固有値を持ちます。

固有ベクトルは、

　\(\begin{pmatrix} a & b \\ b & a \end{pmatrix}\begin{pmatrix} x \\ y \end{pmatrix}=\lambda \begin{pmatrix} x \\ y \end{pmatrix}\) から　\(\begin{pmatrix} a-\lambda & b \\ b & a-\lambda \end{pmatrix}\begin{pmatrix} x \\ y \end{pmatrix}=\begin{pmatrix} 0 \\ 0 \end{pmatrix}\)

\(\lambda_1=a+b\) のとき
\(\begin{pmatrix} -b & b \\ b & -b \end{pmatrix}\begin{pmatrix} x \\ y \end{pmatrix}=\begin{pmatrix} 0 \\ 0 \end{pmatrix}\)　\(\therefore~ \begin{pmatrix} x_{1} \\ y_{1} \end{pmatrix}=s \begin{pmatrix} 1 \\ 1 \end{pmatrix}\)
\(\lambda_2=a-b\)のとき.
\(\begin{pmatrix} b & b \\ b & b \end{pmatrix}\begin{pmatrix} x \\ y \end{pmatrix}=\begin{pmatrix} 0 \\ 0 \end{pmatrix}\)　　 \(\therefore~ \begin{pmatrix} x_{2} \\ y_{2} \end{pmatrix}=t \begin{pmatrix} -1 \\ 1 \end{pmatrix}\)

\(P=\begin{pmatrix} x_{1} &x_{2} \\ y_{1} &y_{2} \end{pmatrix}\)　は直交行列なので、\(s=t=\dfrac{1}{\sqrt{2}}\)

\(P=\begin{pmatrix} 1/\sqrt{2} & -1/\sqrt{2} \\ 1/\sqrt{2} & 1/\sqrt{2} \end{pmatrix}\)　の重要な性質としては、

行列式 \(\Delta=ad-bc=1\).
転置行列 \(P'=\begin{pmatrix} 1/\sqrt{2} & 1/\sqrt{2} \\ -1/\sqrt{2} & 1/\sqrt{2} \end{pmatrix}\) が Pの逆行列.
即ち、\(P^{-1}=P'\)

　（変換としては、Pは原点周りの45度回転を表します）

　　　\(P^{-1}\Sigma P =\begin{pmatrix} \lambda_{1} & 0 \\ 0 & \lambda_{2} \end{pmatrix}\) から　\(\Sigma =P\begin{pmatrix} \lambda_{1} & 0 \\ 0 & \lambda_{2} \end{pmatrix}P^{-1}\) .

目的は、

　　「分散共分散行列\(\Sigma\) に対し、\(AA'=\Sigma\)　を満たす\(A\)を求めること」

であるので、

　\(\Sigma=P\begin{pmatrix} \lambda_{1} & 0 \\ 0 & \lambda_{2} \end{pmatrix}P^{-1}= P \begin{pmatrix} \sqrt{\lambda_{1}} & 0 \\ 0 & \sqrt{\lambda_{2}} \end{pmatrix} \begin{pmatrix} \sqrt{\lambda_{1}} & 0 \\ 0 & \sqrt{\lambda_{2}} \end{pmatrix} P^{-1}\)

ここで \(A=P \begin{pmatrix} \sqrt{\lambda_{1}} & 0 \\ 0 & \sqrt{\lambda_{2}} \end{pmatrix}\)　とおくと、\(\Sigma=AA'\)　となります。

\(A\)は、x軸方向・y方向の拡大縮小と直交変換（原点の周りの回転）の合成です。
　

＜まとめ＞

　\(Z\sim N_2(0,I_2)\)　に対し、

　　　\(X=AZ+m\) とすると　\(X\sim N_2(m,AA')\)、

分散共分散行列 \(\Sigma\) が与えられたとき、

標準正規分布\(Z\)に　\(\Sigma=AA'\)　を満たす \(\Sigma\) の平方根\(A\)を乗ずることで、
\(N_2(m,\Sigma)\) に従う\(X\)を構成することができる。

イメージとしては、

2変量標準正規分布の円状分布を x軸方向・y方向の拡大縮小することで楕円分布となり、
回転することで楕円を回転した形状の分布となります。

つまり、

１変量では拡大縮小だけでしたが、
２変量では回転などの直交変換が加わります。

［注意］

対称行列タイプ２の場合、計算は難しくなります。

\(\Sigma\) が異なる２つの正の固有値を持つ場合、\(\Sigma\) は正定値対称行列と呼ばれます。

\(A=P \begin{pmatrix} \sqrt{\lambda_{1}} & 0 \\ 0 & \sqrt{\lambda_{2}} \end{pmatrix}\)　は、x軸方向・y方向の拡大縮小と直交変換（原点の周りの回転）の合成となりますが、回転角は直交行列 P によって変わります。

2.5 一意性の問題.

１変量の標準化の類比により、2変量正規分布を構成することができました。

しかし、\(\Sigma=AA'\)を満たす A は一意ではありません。

\(\Sigma\) が正定値行列のとき、

　　　\(A=P \begin{pmatrix} \sqrt{\lambda_{1}} & 0 \\ 0 & \sqrt{\lambda_{2}} \end{pmatrix}P^{-1}\)　もまた、\(\Sigma=AA'\)　を満たします。

\(A\) は対称行列であり、

　　\(A\begin{pmatrix} x_{1} \\ y_{1} \end{pmatrix} =P\begin{pmatrix} \sqrt{\lambda_{1}} & 0 \\ 0 & \sqrt{\lambda_{2}} \end{pmatrix}P^{-1} \begin{pmatrix} x_{1} \\ y_{1} \end{pmatrix} =\sqrt{\lambda_{1}} \begin{pmatrix} x_{1} \\ y_{1} \end{pmatrix}\)

　　\(A\begin{pmatrix} x_{2} \\ y_{2} \end{pmatrix} =P\begin{pmatrix} \sqrt{\lambda_{1}} & 0 \\ 0 & \sqrt{\lambda_{2}} \end{pmatrix}P^{-1} \begin{pmatrix} x_{2} \\ y_{2} \end{pmatrix} =\sqrt{\lambda_{2}} \begin{pmatrix} x_{2} \\ y_{2} \end{pmatrix}\)

\(A\) の固有ベクトルは \(\Sigma\) と同じですが、固有値は \(\sqrt{\lambda_1},~\sqrt{\lambda_2}\) になります。

\(A\) による変換では、

変換は互いに垂直な固有ベクトル方向への\(\sqrt{\lambda_1},~\sqrt{\lambda_2}\) 倍の拡大縮小になり、
変換により固有ベクトル方向の分散は\({\lambda_1},~{\lambda_2}\) 倍になります。

拡大縮小は２つの方向になりますが、それぞれの方向の分散は拡大縮小倍率の２乗倍になります。
この仕組みは１変量の一次式の分散が係数の２乗倍になることと同じです。

A による一次変換では、固有ベクトル方向の分散は固有値倍になる。

２変量標準正規分布では各変数の分散は１なので、A によって変換すれば固有ベクトル方向の分散は \({\lambda_1},~{\lambda_2}\) になります。

＜対称行列タイプ１＞

分散共分散行列が \(\Sigma= \begin{pmatrix} a & b \\ b & a \end{pmatrix}\) となる２変量正規分布を、2変量標準正規分布からの一次変換によって構成します。

一次変換として、
　　　\(A=P \begin{pmatrix} \sqrt{\lambda_{1}} & 0 \\ 0 & \sqrt{\lambda_{2}} \end{pmatrix}P^{-1}\)　を用いることにしましょう。　

拡大縮小の方向を与える直交行列 P と倍率を与える固有値が決まれば、\(A\) は決まります。

対称行列タイプ１の場合、\(\lambda_1=a+b\)、\(\lambda_2=a-b\) であり、

　　　\(P=\begin{pmatrix} 1/\sqrt{2} & -1/\sqrt{2} \\ 1/\sqrt{2} & 1/\sqrt{2} \end{pmatrix}\)　　\(P^{-1}=\begin{pmatrix} 1/\sqrt{2} & 1/\sqrt{2} \\ -1/\sqrt{2} & 1/\sqrt{2} \end{pmatrix}\).

y=x 方向の標準偏差 \(\sigma_1\)、y=-x 方向の標準偏差 \(\sigma_2\) とする2変量正規分布は、
2変量標準正規分布 \(Z \sim N(0,I_2)\) に対し.

　\(A=P \begin{pmatrix} \sigma_1 & 0 \\ 0 & \sigma_2 \end{pmatrix}P^{-1} =\dfrac{1}{2} \begin{pmatrix} 1 & -1 \\ 1 & 1 \end{pmatrix} \begin{pmatrix} \sigma_1 & 0 \\ 0 & \sigma_2 \end{pmatrix} \begin{pmatrix} 1 & 1 \\ -1 & 1 \end{pmatrix}\).

を用い　\(X=AZ\)　によって構成できます。

（但し、\(\sigma_1=\sqrt{a+b}\)、　\(\sigma_2=\sqrt{a-b}\) ）

このとき、 \(A=\dfrac{1}{2} \begin{pmatrix} \sigma_1+\sigma_2 & \sigma_1-\sigma_2 \\ \sigma_1-\sigma_2 & \sigma_1+\sigma_2 \end{pmatrix}\)　となります。

シミュレーションで確かめましょう。

タイプ１は計算が容易です。

\(\Sigma=\begin{pmatrix} 1 & 0.5 \\ 0.5 & 1 \end{pmatrix}\)
を分散共分散行列とする2変量正規分布を構成します。.

a <- 1
b <- 0.5
Sigma <- matrix(c(a, b, b, a), 2, 2)
Sigma

##      [,1] [,2]
## [1,]  1.0  0.5
## [2,]  0.5  1.0

固有値は　\(k1=a+b=1.5\)、　\(k2=a-b=0.5\) です。

行列 A を計算します。

k1 <- a+b
k2 <- a-b
K <- matrix(c(sqrt(k1), 0, 0, sqrt(k2)), 2, 2)
P <- matrix(c(1, 1, -1, 1), 2, 2)
A <- P %*% K %*% t(P)*(1/2)
A

##           [,1]      [,2]
## [1,] 0.9659258 0.2588190
## [2,] 0.2588190 0.9659258

一方、

　\(A=\dfrac{1}{2} \begin{pmatrix} \sigma_1+\sigma_2 & \sigma_1-\sigma_2 \\ \sigma_1-\sigma_2 & \sigma_1+\sigma_2 \end{pmatrix} =\dfrac{1}{2} \begin{pmatrix} \sqrt{1.5}+\sqrt{0.5} & \sqrt{1.5}-\sqrt{0.5} \\ \sqrt{1.5}-\sqrt{0.5} & \sqrt{1.5}+\sqrt{0.5} \end{pmatrix}\).

\(\sqrt{1.5}+\sqrt{0.5}=0.9659258\)で合致しています。

標準正規乱数 data を A によって一次変換し描画します。

d_new <- A %*% t(data)
d_new <- data.frame(X=d_new[1,],Y=d_new[2,])
par(mfrow=c(1,2))
plot(data,asp=1,col="lightblue",cex=0.1,xlim=c(-3,3),ylim=c(-3,3))
abline(h = 0, v=0,col = "lightblue", lty = 1)
plot(d_new,col="pink",cex=0.1,xlim=c(-3,3),ylim=c(-3,3),asp=1)
abline(a = 0, b=1,col = "pink", lty = 1)
abline(a = 0, b=-1,col = "pink", lty = 1)

左は標準正規乱数、右は変換後のデータです。
変換後のデータの分散共分散行列を cov( ) 関数で計算します。

round(cov(d_new),1)

##     X   Y
## X 1.0 0.5
## Y 0.5 1.0

確かに、目的とした分散共分散行列に一致しています。

＜対称行列タイプ２＞

タイプ２では、固有ベクトルは y=x,y=-x 方向とはなりませんが、直交します。

対称行列は、直交する固有ベクトル方向への固有値倍の拡大縮小であり.
その平方根は、直交する固有ベクトル方向への固有値の平方根倍の拡大縮小.

であることは同じです。.

3 コレスキー分解.

前節では、

　「固有値と固有ベクトルを用いて \(\Sigma=AA'\) を満たす A を求める」

方法を学習しました。

本節では、固有値と固有ベクトルに頼らない方法として、コレスキー分解を紹介します。

コレスキー分解は、\(\Sigma=AA'\)を満たす \(A\) を下三角行列として求める方法です。

まず、 \(\Sigma\) が相関行列のときを考えます。

\(\Sigma= \begin{pmatrix} 1 & \rho \\ \rho & 1 \end{pmatrix}\:,\:\: A = \begin{pmatrix} \, l_{11} & 0 \\ \, l_{21} & l_{22}\\ \end{pmatrix}\)　とおくと、 .

\(A=\begin{pmatrix} 1 & 0 \\ \rho & \sqrt{1-\rho^{2}} \end{pmatrix}\) のとき、\(\Sigma=AA'\) となります。
実際、
　　\(\begin{pmatrix} 1 & \rho \\ \rho & 1 \end{pmatrix}= \begin{pmatrix} l_{11} & 0\\ l_{21} & l_{22}\\ \end{pmatrix} \begin{pmatrix} l_{11} & l_{21}\\ 0 & l_{22}\\ \end{pmatrix} = \begin{pmatrix} l_{11}^2 & l_{11}l_{21}\\ l_{11}l_{21} & l_{21}^2+l_{22}^2\\ \end{pmatrix}\)

　　\(1= l_{11}^2\:\:,\:\:\rho = l_{11}l_{21}\)　より、　　\(l_{11}=1\:\:,\:\:l_{21}=\rho\).

　　\(1= l_{21}^2+l_{22}^2\)　　より　　 \(l_{22}=\sqrt{1-\rho^2}\)

解は一意ではありませんが、計算は簡便で、次元が上がっても容易に計算できます。
このように、

相関行列のコレスキー行列は、相関係数で表されます。

3.1 相関行列の分解.

\(\Sigma=\begin{pmatrix} 1 & 0.5 \\ 0.5 & 1 \end{pmatrix}\)
を分散共分散行列とする2変量正規分布をコレスキー分解で構成します。

b <- 0.5
A <- matrix(c(1, b, 0, sqrt(1-b^2)), 2, 2)
A

##      [,1]      [,2]
## [1,]  1.0 0.0000000
## [2,]  0.5 0.8660254

相関係数から計算したコレスキー行列です。
\(AA'\) が \(\Sigma\) になることを確かめます。

Sigma <- A %*% t(A)
Sigma

##      [,1] [,2]
## [1,]  1.0  0.5
## [2,]  0.5  1.0

コレスキー行列は y 軸を不動直線に持ち、対角成分が固有値になっていることも確かめられます。

eigen <- eigen(A)
v <- eigen$vectors
v

##           [,1] [,2]
## [1,] 0.2588190    0
## [2,] 0.9659258    1

e <- eigen$values
e

## [1] 1.0000000 0.8660254

\(Z \sim N(0,I_2)\) に対し、コレスキー行列 \(A\) によって \(X=AZ\) と変換すると、

　　　\(X\sim N_2(0,\Sigma)\) となることを確かめておきましょう。

n <- 5000
library(MASS)
mu <- c(0, 0)
Sigma <- matrix(c(1, 0, 0, 1), 2, 2)
data <- mvrnorm(n, mu, Sigma)
d <- data.frame(X=data[,1],Y=data[,2])
d_new <- A %*% t(data)
d_new <- data.frame(X=d_new[1,],Y=d_new[2,])
par(mfrow=c(1,2))
plot(data,asp=1,col="lightblue",cex=0.1,xlim=c(-3,3),ylim=c(-3,3))
abline(h = 0, v=0,col = "lightblue", lty = 1)
plot(d_new,col="pink",cex=0.1,xlim=c(-3,3),ylim=c(-3,3),asp=1)
abline(a = 0, b=1,col = "pink", lty = 1)
abline(a = 0, b=-1,col = "pink", lty = 1)

変換されたデータ d_new の分散共分散行列は \(\Sigma\) となります。

round(cov(d_new),1)

##     X   Y
## X 1.0 0.5
## Y 0.5 1.0

目的とした分散共分散行列に一致しています。

3.2 分散共分散行列の分解

一般の分散共分散行列 \(\Sigma\) については、

　　　\(\Sigma=\begin{pmatrix} \sigma_x^2 & \sigma_{xy}\\ \sigma_{xy} & \sigma_y^2 \end{pmatrix} =\begin{pmatrix} \sigma_x^2 & \rho\sigma_x \sigma_y \\ \rho\sigma_x \sigma_y & \sigma_y^2 \end{pmatrix}\)　が成立します。

\(\Sigma= \begin{pmatrix} \sigma_x & 0 \\ 0 & \sigma_y \end{pmatrix} \begin{pmatrix} 1 & \rho \\ \rho & 1 \end{pmatrix} \begin{pmatrix} \sigma_x & 0 \\ 0 & \sigma_y \end{pmatrix}\)　であるので、相関行列のコレスキー分解を用いて分解ができます。
即ち、\(A=\begin{pmatrix} \sigma_x & 0 \\ 0 & \sigma_y \end{pmatrix} \begin{pmatrix} 1 & 0 \\ \rho & \sqrt{1-\rho^{2}} \end{pmatrix}\)　が \(\Sigma\) のコレスキー分解です。

\(\Sigma=\begin{pmatrix} 1 & 0.5 \\ 0.5 & 4 \end{pmatrix}\)
を分散共分散行列とする2変量正規分布をコレスキー分解で構成します。

Sigma <- matrix(c(1, 0.5, 0.5, 4), 2, 2)
Sigma

##      [,1] [,2]
## [1,]  1.0  0.5
## [2,]  0.5  4.0

\(S=\begin{pmatrix} \sigma_x & 0 \\ 0 & \sigma_y \end{pmatrix} =\begin{pmatrix} 1 & 0 \\ 0 & 2 \end{pmatrix}\)、 \(R=\begin{pmatrix} 1 & \rho \\ \rho & 1 \end{pmatrix}\) として \(SRS\) を計算し、\(\Sigma\) となることを確かめます。
（ \(\sigma_x=1,~\sigma_y=2,~\rho\sigma_x \sigma_y=0.5\)　より、　\(\rho=0.25\) ）

S <- matrix(c(1, 0, 0, 2), 2, 2)
R <- matrix(c(1, 0.25, 0.25, 1), 2, 2)
S %*% R %*% S

##      [,1] [,2]
## [1,]  1.0  0.5
## [2,]  0.5  4.0

コレスキー行列を求めます。

b <- 0.25
A <- matrix(c(1, b, 0, sqrt(1-b^2)), 2, 2)
S %*% A

##      [,1]     [,2]
## [1,]  1.0 0.000000
## [2,]  0.5 1.936492

これが \(\Sigma\) のコレスキー分解の行列です。

R の chol( ) 関数で計算することもできます。
chol( )関数は上三角行列を出力するので、下三角行列にするには転置する必要があります。

chol(Sigma)

##      [,1]     [,2]
## [1,]    1 0.500000
## [2,]    0 1.936492

A <- t(chol(Sigma))
A

##      [,1]     [,2]
## [1,]  1.0 0.000000
## [2,]  0.5 1.936492

A %*% t(A)

##      [,1] [,2]
## [1,]  1.0  0.5
## [2,]  0.5  4.0

\(A\) の固有値、固有ベクトルを求めます。

eigen <- eigen(A)
v <- round(eigen$vectors,1)
v

##      [,1] [,2]
## [1,]    0  0.9
## [2,]    1 -0.5

e <- eigen$values
e

## [1] 1.936492 1.000000

d_new <- A %*% t(data)
d_new <- data.frame(X=d_new[1,],Y=d_new[2,])
par(mfrow=c(1,2))
plot(data,asp=1,col="lightblue",cex=0.1,xlim=c(-3,3),ylim=c(-3,3))
abline(h = 0, v=0,col = "lightblue", lty = 1)
plot(d_new,col="pink",cex=0.1,asp=1)
abline(a = 0, b=v[2,2]/v[1,2],col = "lightblue", lty = 1)
abline(v = 0,col = "lightblue", lty = 1)
arrows(x0=0, y0=0, x1=v[1,1],  y1=v[2,1], code=2, col="blue")
arrows(x0=0, y0=0, x1=v[1,2],  y1=v[2,2], code=2, col="red")

3.3 白色化.

1変量正規分布では、

\(X\sim N(\mu,\sigma^2)\)　に対し、\(Z=\dfrac{X-\mu}{\sigma}\sim N(0,1)\)
によって、\(X\) を標準正規分布に変換することができます。

この操作を標準化といいます。これは、偏差を分散の平方根で割る計算です。

２変量量正規分布では、変量相互の相関があるので注意が必要です。
各変数の偏差を分散の平方根で割っても変量間の相関が残るからです。

二変量正規分布には、次の２通りの定義の方法があります。

確率密度関数による定義.

　　　平均 \(\mu\)、分散共分散行列 \(\Sigma\) に対し、
　　　　\(f(X)=\dfrac{1}{\left(\sqrt{2\pi}\right)^2\sqrt{|\Sigma|}}\exp\left[{-\dfrac{1}{2}(X-\mu)'\ \Sigma^{-1}(X-\mu)}\right]\).
　　　　このとき、\(X\) は２変量正規分布 \(N(\mu,\Sigma)\) に従います。

一次変換による定義.

　　　２次元標準正規分布 Z 、正則行列 A に対し、
　　　　\(X=AZ+\mu \sim N(\mu,AA')\).
　　　　\(Σ=AA′\) は分散共分散行列です。

一次変換による定義を用いると、

\(X\sim N(\mu,AA')\)　に対し、\(Z=A^{-1}(X-\mu)\)
とすれば、Z は２変量標準正規分布に従います。
この操作を白色化と呼びます。

白色化は、すべての周波数を均等に含む光は白く見えることに由来します。
二変量標準正規分布では、各成分は独立で均等です。

\(Σ=AA′\) を満たす A を求める方法として、

固有値と固有ベクトルを用いる方法.
コレスキー分解を用いる方法.

について学びました。
実際の応用については、多変量解析や機械学習の分野で学習してください。

高校生の線形代数_3