主成分の分散最大化のグラフを描いてみようと思う
Rの操作方法の練習として、前回の主成分分析の学習ノート（前回のエントリへ）のところに出てきた、主成分の分散を表す式、
・・・(9.7)
のグラフを描いてみることにする。とか
が動くと
がどういうふうに動くの化を、可視化するわけである。rは相関係数なのだが、適当に0.5とでもしておこう。
その上で、
・・・(9.8)
という制約条件を付けることの意味も、可視化してみようと思う。
それにより、分散の「最大値」をグラフからイメージしてみることにする。
分散をあらわす平面の描画
まず、上記の数式をRが読める形で記述すると、
# V = a1^2 + a2^2 + (2 * r * a1 * a2) ←(9.7)式
# a1^1 +a2^2 = 1 ←(9.8)式
となる点に注意する。
さて、f(x, y)という2変数からなる関数の描画には、Rではpersp()という関数*1を使うのだが、これはけっこう、使い方のセオリーを知ってないといけない代物だ。
カンタンにいうと、最初にx軸・y軸の平面に格子を描いて、この格子を２変数関数にしたがってz方向に変形していくわけである。
まず、昇順かつ等間隔になるようにxとyのデータを生成して、格子の座標を与えてみる。
> a1 <- seq(-1, 1, length=30) # -1から1までの範囲で、等間隔に30個のデータを生成する。
> a2 <- a1 # 格子をつくるだけだから、a1と同じでよい。
>
> print(a1) # a1がどんなデータか一応みておく
（ブログに貼ると汚いので省略）
次に、2変数関数f(x, y)を定義する。
> f <- function(a1, a2){a1^2 + a2^2 + (2 * 0.5 * a1 * a2)}
> # 0.5というのは相関係数を適当に置いてみたもの。
で、今度はVの軸に入れるデータを得なければならないわけだが、これにはouter()という関数を使う。この関数を使うと、あるxとyに対応するf(x, y)の値を計算して行列を与えてくれる。
> V <- outer(a1, a2, f) # 3つめの引数が関数を表している。
> print(V[1:5, 1:5]) # zがどんな行列かみてみる（ただしデータが多いので５行・５列目まで）
[,1] [,2] [,3] [,4] [,5]
[1,] 3.000 2.798 2.605 2.422 2.249
[2,] 2.798 2.600 2.413 2.234 2.065
[3,] 2.605 2.413 2.229 2.056 1.892
[4,] 2.422 2.234 2.056 1.887 1.728
[5,] 2.249 2.065 1.892 1.728 1.573
最後に、persp()関数をつかって、グラフを描画してみる。
persp(a1, a2, V,
xlab="a1", ylab="a2", zlab="V", # 軸のラベルを指定
theta = 30, # 横回転の角度。ここをいじって表示のアングルを変える。
phi = 30, # 縦回転の角度。ここをいじって表示のアングルを変える。
ticktype = "simple", # 線の種類
lwd=0.5, # 線が太いと見づらかったりするので0.5にしておいた。
)
これは、thetaとphiの値をいじることで、角度をかえて眺めることができる。
制約条件と最大値の可視化
さて次に、制約条件を表す式、
を可視化してみる。
この式は、よく考えたら円を表す公式だ。ということは、さっき描画した３次元のグラフでいうと、真上からみての平面に円を描くような感じになるんだろう。
ではまず、円周を表してくれるであろう点の集合のデータをつくってみる。
> a <- rep(seq(-1, 1, length=100), 100) # データ1000個とかにすると計算が遅いのでこれが限界
> a1 <- c(a, a)
> a2 <- c(sqrt(1-a^2), -sqrt(1-a^2))
1つのa1の値に対してa2の値が2つ存在し、原点を中心とする円だから単純に符号の正負を変えればよいという点に注意。
なおこの時点で、さっき格子状の面を描くときに作ったa1、a2というオブジェクトはなくなって入れ替わっている。一連の分析として実行する時は、ほんとはオブジェクト名を分ければいいんだけど。
で、これを一応プロットしてみよう。
> plot(a1, a2, xlab="a1", ylab="a2", type="p", cex=0.5, pch=19) # a1とa2の関係が円であることを確認
たしかに円周になってる。
で、この円の関係にあるa1とa2のデータを、さっきのVの式に入れてみる。
> V <- a1^2 + a2^2 + (2 * 0.5 * a1 * a2) # 相関係数はまた0.5にしといた
これでV、a1、a2のデータが揃ったのでこれを空間上にプロットすればいいのだが、3次元空間上に「点」をプロットするときは、scatterplot3dというパッケージを使う。
> install.packages("scatterplot3d")
> library(scatterplot3d)
あとは超簡単で、３次元の座標に対応するデータ（ベクトル）をscatterplot3d()関数に投げればいいだけだ。
> scatterplot3d(a1, a2, V, type="l", angle=60) # 線の形と眺める角度も引数として与えてある
こういう図になったが、たしかに真上からみれば円になってそうな感じで、最初につくった平面のグラフと円柱が交わるところを結んだような感じのグラフになっている。
分散を最大化させると
の組み合わせが２個あることも、一目瞭然である。（上から見て円になる線のグラフしか見なかった場合は、あくまでこの線は制約範囲の外周を表しているだけなので、ひょっとしたら中心あたりが上向きに突出している面になってるかもと思うことはできるが。）
*1:perspectiveの略かな？
