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

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

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

f:id:StatModeling:20201114161918j:plain

データ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（↓に掲載。ヒントン図と呼ばれます。行列の要素の大小を視覚的に表したもの）の話に相当します。 f:id:StatModeling:20201114163431p:plain

しかしながら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）の行列です。

以下は結果になります。

f:id:StatModeling:20201114161922p:plain

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

StatModeling Memorandum

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

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