はじめに
こちらの本を読みました。
こちらには非負値行列因子分解(NMF)の話も出てきます。
前から、各要素がガウス分布やポアソン分布に従うとする場合が多いですが、ゼロ過剰ポアソン分布に従うとする方が適切だろうと思う場合が多かったです。
そのため、各要素がゼロ過剰ポアソン分布に従う非負値行列因子分解を考えて、Rstanで実装してみました。
ZIP NMF
NMF
非負値行列因子分解(NMF : Nonnegative Matrix Factorization)とは、各要素が非負値から構成されている行列を、
と
の二つの行列に分解する解析手法です。
データを低次元空間に写像する方法であるとも呼ばれています。
ZIP
ゼロ過剰ポアソン分布(Zero-Inflated Poisson model)とは、0が多いポアソン分布のことです。
具体的には、ベルヌーイ分布により0のみの分布?とポアソン分布が切り替わる混合分布になっています。
ここで特筆すべきなのは、0の中にはベルヌーイ分布により発生したされた0と、ポアソン分布から発生した0が混ざっていることです。
ZIP NMF
の各要素がZIPに従うと仮定したNMFを考えます。
確率で0、
でポアソン分布のNMFにあるようなNMFであり、次のような確率分布の式で表現することができます。
[tex: W{d,m} ]と[tex: H{m,n} ]は2つの行列に分解された際の各要素であり、ガンマ分布に従っているとします。
Rstanでやってみる
使うデータはこんな感じのものです。
library(tidyverse) library(rstan) library(Rcpp) library(bayesplot) # 並列化のおまじない rstan_options(auto_write=TRUE) options(mc.cores=parallel::detectCores()) mat = c(5,0,2,2,2,3,0,0,2,0, 3,2,0,1,0,3,0,1,0,3, 4,3,0,0,0,5,0,0,2,0, 1,0,2,3,1,4,1,0,2,0, 2,0,4,2,3,0,0,2,0,0, 1,0,1,3,2,4,0,0,0,1) dat = matrix(mat, 10,6) %>% t
分布を見ると、いい感じに0過剰ポアソン分布に従うようになっています。
data.frame(X = mat) %>% ggplot(aes(x=X))+ geom_histogram(binwidth = 1, alpha=0.5) + ggtitle("Data histogram for NMF")
これに対して、ポアソン分布を仮定したNMFをこちらの記事のStanコードを用いて実行してみました。
潜在因子数は3としました。
復元されたデータの分布は次にようになりました。
0過剰の部分が表現されていません。
復元誤差(RMSE)は次にようになりました。 (stanオブジェクトからパラメータをtidyに扱う方法が見つからなかったので、良い方法があれば教えていただきたいです。)
> (summary(stan.fit)$summary[K:(K+length(mat)-1),"50%"] -mat)^2 %>% mean() %>% sqrt [1] 0.9831921
ここからが、ZIP NMFです。
上記のZIP NMFをstanでモデル化したのは次のような感じです。
data{
int<lower=1> D;
int<lower=1> N;
int<lower=1> M;
int X[D,N];
}
parameters{
matrix<lower=0>[D, M] W;
matrix<lower=0>[M, N] H;
matrix<lower=-10, upper=10>[D, N] theta;
real<lower=0> alpha;
real<lower=0> beta;
}
model{
alpha ~ gamma(1, 2);
beta ~ gamma(1,2);
for(d in 1:D){
for(n in 1:N){
theta[d,n] ~ uniform(-10,10);
}
}
for(d in 1:D){
for(m in 1:M){
W[d,m] ~ gamma(alpha, beta);
}
}
for(m in 1:M){
for(n in 1:N){
H[m,n] ~ gamma(alpha, beta);
}
}
for(d in 1:D){
for(n in 1:N){
if(X[d,n] == 0){
target += log_sum_exp(
bernoulli_lpmf(0|inv_logit(theta[d, n])),
bernoulli_lpmf(1|inv_logit(theta[d, n])) + poisson_lpmf(0|W[d,:]*H[:,n])
);
}else {
target += bernoulli_lpmf(1|inv_logit(theta[d, n])) + poisson_lpmf(X[d,n]|W[d,:]*H[:,n]);
}
}
}
}
generated quantities {
int Xpred[D,N];
real Xlogit[D,N];
for(d in 1:D){
for(n in 1:N){
Xlogit[d,n] = inv_logit(theta[d, n]);
if(Xlogit[d,n] < 0.5){
Xpred[d,n] = 0;
}else{
Xpred[d,n] = poisson_rng(W[d,:]*H[:,n]);
}
}
}
}
キックするコードは次にような感じです。 潜在因子数は3としました。
run_data <- list(D=nrow(dat),N=ncol(dat),M=3,X=dat) stan.fit <- stan("ZIP_NMF.stan", data = run_data, seed = 42, iter = 4000, warmup = 2000, chains = 4)
収束したか見てみます。
mcmc_rhat(rhat(stan.fit))
ちゃんと収束しているようにですね
サンプリングに自己相関がないことも確認しました。
復元された値の分布を確認してみました。
0過剰の部分が表現できているようです。
復元誤差も確認してみました。
> (summary(stan.fit)$summary[K:(K+length(mat)-1),"50%"] -mat)^2 %>% mean() %>% sqrt [1] 0.8164966
単にポアソン分布を仮定した場合よりも当てはまりが良い状態にすることができました。
0の発生の潜在構造を捉えていないため、ベルヌーイ分布のパラメータの表列を行列分解により表現するするように拡張する必要があるかなと感じました。
