PCA（主成分分析）のド発展版に相当する、ガウス過程を用いたGPLVMをRからサクッと使うまでの備忘録です。
GPLVMの説明で分かりやすいのは、以下の統計数理研究所のH26年度公開講座「ガウス過程の基礎と応用」の持橋先生と大羽先生の発表資料です。
- [1] 統計数理研究所 H26年度公開講座「ガウス過程の基礎と応用」 (web) 元論文は以下です。
- [2] M. K. Titsias and N. D. Lawrence (2010) Bayesian Gaussian Process Latent Variable Model. Thirteenth International Conference on Artificial Intelligence and Statistics (AISTATS), JMLR: W&CP 9, pp. 844-851. (pdf)
そして、提唱者であるLawrenceさんがグラントをとって、GPLVMを実装してライブラリ（vargplvm）をGithub上に公開しています。メインはMatlabコードです。GithubのURLとR版とREADMEのURLは以下です。
- [3] https://github.com/SheffieldML/vargplvm
- [4] https://github.com/SheffieldML/vargplvm/tree/master/vargplvmR
- [5] https://github.com/SheffieldML/vargplvm/blob/master/vargplvm/matlab/vargplvm_README.txt
RコードはMatlabコードをJie Haoさん（Nicolas Durrandeさんが協力）がC++（内部）とR（インターフェース）に移植したようです。dynamical modelの部分と発展的なデモは移植されていないので、何かあったらMatlabコードの方を見る必要があるでしょう。Stanもそうですが、こういうグラントdrivenでいいライブラリが出てくるの嬉しいですね。何よりグラントの結果は、そうそう簡単に作者がキレて消えたりしないのがいいところです。
このライブラリは推論部分に変分ベイズを使っており、一般的にデータポイントが1000を超えてくると段々つらくなるガウス過程にもかかわらず、動作はとても高速です。例えば、Rコードを使った限りでは、3クラスの1000ポイントx12特徴量のデータぐらいならショボいcore i5とかでも計算は数分です。
それではインストール方法から説明します。まずC++コードをコンパイルするので、Windows環境ならRoolsなどをインストールしておく必要があります。次に{devtools}をインストールしていない人は以下でインストールしてください。
install.packages("devtools")
次にGithubから{vargplvm}をインストールします。
devtools::install_github("SheffieldML/vargplvm/vargplvmR")
次に使い方です。vargplvmはRのライブラリとしてはかなり不親切で、??vargplvmとかして関数の説明見てもサッパリ分かりません。R版のデモが2つありますのでそれらを参考にしました。
https://github.com/SheffieldML/vargplvm/blob/master/vargplvmR/R/demOilVargplvm1.R
https://github.com/SheffieldML/vargplvm/blob/master/vargplvmR/R/demmESCVargplvm.R
以下では説明のためirisでやります（ほんとはもっとデータがあって特徴量も多いものでないとGPLVMの意味がないです）。
library(vargplvm) set.seed(1) Y <- as.matrix(scale(iris[,1:4])) options <- vargplvmOptions("dtcvar") options$kern <- c("rbfard2", "bias", "white") options$numActive <- 10 options$optimiser <- "scg" latentDim <- 3 model <- vargplvmCreate(latentDim, ncol(Y), Y, options) model <- vargplvmParamInit(model, model$m, model$X) covD <- dim(model$vardist$covars) model$vardist$covars <- 0.5*matrix(1, covD[1], covD[2]) + 0.001*matrix(rnorm(covD[1]*covD[2]), covD[1], covD[2]) model$learnBeta <- 1 model <- vargplvmOptimise(model, display=1, iters=2) mm <- vargplvmReduceModel(model,2) library(ggplot2) df <- data.frame(X1=mm$X[,1], X2=mm$X[,2], Species=iris$Species) p <- ggplot(data=df, aes(x=X1, y=X2, color=Species)) + geom_point(alpha=0.5) ggsave(p, file='result.png', dpi=300, w=4, h=3)
- 4行目:データは
matrixでデータポイント×特徴量の形でないといけません。念のためscalingしていますが、デモを見る限りscaleしなくてもよさそうです。 - 6行目:カーネルを指定しています。カーネルを変えるときは第一引数の
rbfard2のところを変えます。 - 7行目:持橋先生の資料で言うところの43ページ目の「active set」の数を指定します。データポイントの数以下でないといけません。
- 9行目:潜在変数の次元を指定します。特徴量の数以下でないといけません。ちなみに、デモでは1000x12のデータに対して、
numActive=50、latentDim=10となっていました。 - 19行目:実際の処理です。itersを増やすとどんどん細かい刻み幅にして最適化をすすめていきます。そんなに増やさなくてもよさそうです。
- 20行目:計算した潜在変数の空間を射影しているようです。詳しいことはもう少し読み込まないと分かりません。すみません。
- 24行目:よく論文に出てきている図は、この射影後の1成分目と2成分目をプロットしているものになります。
結果は以下です。
濃淡を描くときは、mm$vardist$meansとmm$vardist$covarsを元に描いていると思います。
