StatsBeginner: 初学者の統計学習ノート

統計学およびR、Pythonでのプログラミングの勉強の過程をメモっていくノート。たまにMacの話題。

Rで3次元グラフを描く練習

主成分の分散最大化のグラフを描いてみようと思う

 Rの操作方法の練習として、前回の主成分分析の学習ノート(前回のエントリへ)のところに出てきた、主成分の分散を表す式、


 { \displaystyle

V_{z_1} = {a_1}^{2} {a_2}^{2} + 2r_{{x_1}{x_2}}{a_1}{a_2}

}  ・・・(9.7)


 のグラフを描いてみることにする。{a_1}とか{a_2}が動くと{V_{z_1}}がどういうふうに動くの化を、可視化するわけである。rは相関係数なのだが、適当に0.5とでもしておこう。
 その上で、


 { \displaystyle {a_1}^{2} + {a_2}^{2} = 1 }  ・・・(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軸の平面に格子を描いて、この格子を2変数関数にしたがって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がどんな行列かみてみる(ただしデータが多いので5行・5列目まで)
   [,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にしておいた。
)


f:id:midnightseminar:20140727005940p:plain:w400


 これは、thetaとphiの値をいじることで、角度をかえて眺めることができる。


 f:id:midnightseminar:20140727010353p:plain:w400
 f:id:midnightseminar:20140727010408p:plain:w400
 
 

制約条件と最大値の可視化

 さて次に、制約条件を表す式、


 { \displaystyle {a_1}^{2} + {a_2}^{2} = 1 }


 を可視化してみる。
 この式は、よく考えたら円を表す公式だ。ということは、さっき描画した3次元のグラフでいうと、真上からみて{ ({a_1}, {a_2}) }の平面に円を描くような感じになるんだろう。


 ではまず、円周を表してくれるであろう点の集合のデータをつくってみる。

> 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の関係が円であることを確認


 f:id:midnightseminar:20140727012757p:plain


 たしかに円周になってる。
 で、この円の関係にある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)


 あとは超簡単で、3次元の座標に対応するデータ(ベクトル)をscatterplot3d()関数に投げればいいだけだ。

> scatterplot3d(a1, a2, V, type="l", angle=60) # 線の形と眺める角度も引数として与えてある


 f:id:midnightseminar:20140727011554p:plain:w400


 こういう図になったが、たしかに真上からみれば円になってそうな感じで、最初につくった平面のグラフと円柱が交わるところを結んだような感じのグラフになっている。
 分散を最大化させる{a_1}{a_2}の組み合わせが2個あることも、一目瞭然である。(上から見て円になる線のグラフしか見なかった場合は、あくまでこの線は制約範囲の外周を表しているだけなので、ひょっとしたら中心あたりが上向きに突出している面になってるかもと思うことはできるが。)

*1:perspectiveの略かな?