読者です 読者をやめる 読者になる 読者になる

StatModeling Memorandum

StanとRでベイズ統計モデリングします. たまに書評.

infer.netの例題シリーズ4 Bayesian PCA

BUGS R

今回は「Short Examples: Bayesian PCA and Factor Analysis」をやります。みんな大好きPRMLの「12章 連続潜在変数 | 12.2 確率的主成分分析 | 12.2.3 ベイズ的主成分分析」にバッチリ対応します。

infer.netの元ネタと同一の真のパラメータからデータを作って推定してみます。グラフィカルモデルは以下の通りです(infer.netの図を引用)。

データY(図中X)は1000個×10次元(N×D)です。真の潜在変数z.m(図中z)は1000個×3component(N×M)の行列で、線形変換の重みtrue.w.m(図中W)は3component×10次元(M×D)の行列です。

infer.netの解説ではいったんM=6componentとしてwを推定し、wの要素の絶対値の平均が小さい行を不要なcomponentとみなせる、という話の展開をしています。これはPRMLではp.301 図12.14(ヒントン図; 行列の要素の大小を視覚的に表したもの)の話に相当します。しかしながらcomponentの数の決め方は現時点で決定打があるわけではない様子で、今回infer.netではこれを試した、という位置づけのようです。

BUGSコードは以下の通りです。

model {
   for (n in 1:N) {
      for (d in 1:D) {
         Y[n, d] ~ dnorm(u[n, d], tau[d])
         u[n, d] <- t[n, d] + mu[d]
         t[n, d] <- inprod(z[n,], w[,d])
      }
   }
   for (d in 1:D) {
      tau[d] ~ dgamma(2.0, 2.0)
      mu[d] ~ dnorm(0, 0.01)
   }
   for (n in 1:N) {
      for (m in 1:M) {
         z[n, m] ~ dnorm(0, 1)
      }
   }
   for (m in 1:M) {
      alpha[m] ~ dgamma(2.0, 2.0)
      for (d in 1:D){
         w[m, d] ~ dnorm(0, alpha[m])
      }
   }
}

今回はJAGSを使いました。6行目、BUGSには行列の積の演算子がないためinprod()で対応しています。実行するRコードは以下になります。

library(MASS)
library(rjags)
library(R2WinBUGS)

set.seed(123)

get.example.data <- function(N){
   true.w.m <- matrix(c(-0.30, 0.40, 0.20, -0.15, 0.20, -0.25, -0.50, -0.10, -0.25, 0.10, -0.10, -0.20, 0.40, 0.50, 0.15, -0.35, 0.05, 0.20, 0.20, -0.15, 0.15, 0.05, 0.15, -0.10, -0.15, 0.25, -0.10, 0.15, -0.30, -0.55), 3, 10, byrow=T)
   true.mu.v <- c(-0.95, 0.75, -0.20, 0.20, 0.30, -0.35, 0.65, 0.20, 0.25, 0.40)
   true.tau.v <- c(8.0, 9.0, 10.0, 11.0, 10.0, 9.0, 8.0, 9.0, 10.0, 11.0)
   M <- nrow(true.w.m)
   D <- ncol(true.w.m)
   z.m <- matrix(rnorm(N*M, 0, 1), N, M, byrow=T)
   t.m <- z.m %*% true.w.m
   data <- t(sapply(1:N, function(n){
      u.v <- t.m[n, ] + true.mu.v
      rnorm(D, u.v, 1/sqrt(true.tau.v))
   }))
}

N <- 1000
Y <- get.example.data(N)
D <- ncol(Y)
M <- 6

data <- list(
   N = N,
   Y = Y,
   D = D,
   M = M
)

inits <- function(){
   list(
      z = matrix(rnorm(N*M,0,1),N,M),
      w = matrix(rnorm(M*D,0,1),M,D),
      mu = rep(0, D),
      tau = rep(1, D),
      alpha = rep(1, M)
   )
}

m <- jags.model("model/model.bugs", data, inits, n.chains=3, n.adapt=100)
update(m, 5000)
mcmc.list <- coda.samples(m, c("w", "mu", "tau"), n.iter=200000, thin=200)
  • 8行目:真のwのデータを与えています。3component×10次元(M×D)の行列です。

以下は結果になります。

mutauはかなりうまく推定できています。infer.netに書いてあったようにwの絶対値の大きさで評価すると、w[4,*], w[5,*]のcomponentは落とせそうかなと思いたくなります。しかし、推定されたMCMCサンプルのsdが大きいため、そうはうまくいっていないと思います。次元削減をする場合は次元数もパラメータに含まれるようなモデルがよさそうです。