廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

非負値行列因子分解をRで(ベイズ推論による機械学習入門)2

abrahamcow.hatenablog.com

はい。

機械学習スタートアップシリーズ ベイズ推論による機械学習入門 (KS情報科学専門書)

機械学習スタートアップシリーズ ベイズ推論による機械学習入門 (KS情報科学専門書)

p.177 に「Sの期待値は更新アルゴリズムから除去できます。詳しくは Bayesian Inference for Nonnegative Matrix Factorisation Models」みたいなことが書いてあり、それを(おそらく)論文通りに実装したのがこちらの関数です。

格段に速くなりました。

NMFVB <-function(Y,L=2,alpha=1,beta=1,tol=1,maxit=200,seed=1234){
  set.seed(seed)
  N <- nrow(Y)
  K <- ncol(Y)
  M <- (!is.na(Y))*1
  Y[is.na(Y)] <- 0
  EW <- EelW <- matrix(rgamma(N*L,shape=alpha,scale=beta),N,L)
  EH <- EelH <- matrix(rgamma(L*K,shape=alpha,scale=beta),L,K)
for(iter in 2:maxit){
  Sw <- EW * (((Y*M)/(EelW %*% EelH)) %*% t(EelH))
  Sh <- EH * (t(EelW) %*% ((Y*M)/(EelW %*% EelH)))
  beta_W <- 1/(alpha/beta+M%*%t(EH))
  alpha_W <- alpha + Sw
  EW <- alpha_W*beta_W
  beta_H <- 1/(alpha/beta+t(EW)%*%M)
  alpha_H <- alpha + Sh
  EH <- alpha_H*beta_H
  EelW <-exp(digamma(alpha_W))*beta_W
  EelH <-exp(digamma(alpha_H))*beta_H
}
return(list(W=EW,H=EH,
            alpha_W=alpha_W,beta_W=beta_W,
            alpha_H=alpha_H,beta_H=beta_H))
}

たぶん、コードを眺めてもこのアルゴリズムがなにやってるかはさっぱりわからないと思います。論文を読んでください。

チェ・ゲバラ - Wikipedia の画像をランダムに欠測させて補完します。

library(jpeg)
Che<-readJPEG("CheHigh.jpg")
dim(Che)
mat4img <- function(mat1){
  t(apply(mat1,2,rev))
}

Y2 <- Y <-round(Che*1000)
dimy <-dim(Y)
set.seed(1234)
Y2[cbind(sample.int(dimy[1],1000000,replace = TRUE),
         sample.int(dimy[2],1000000,replace = TRUE))] <- NA
image(mat4img(Y2),axes=FALSE)

欠損させた画像。 1025×800。

f:id:abrahamcow:20180713232131p:plain

基底数3のとき。

out3 <-NMFVB(Y2,L=3)
image(mat4img(out3$W %*% out3$H),axes=FALSE)

f:id:abrahamcow:20180713232229p:plain

10のとき。

out10 <-NMFVB(Y2,L=10)
image(mat4img(out10$W %*% out10$H),axes=FALSE)

f:id:abrahamcow:20180713232243p:plain

30のとき。

out30 <-NMFVB(Y2,L=30)
image(mat4img(out30$W %*% out30$H),axes=FALSE)

f:id:abrahamcow:20180713232257p:plain