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

廿TT

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

ゼロ切断ポアソン分布によるセッション数のモデル化とデータ拡大による潜在利用者数の推定

ゼロ切断ポアソン分布によるセッション数のモデル - 廿TT ではウェブサイトのセッション数をゼロ切断ポアソン分布としてモデル化し、ウェブサイトにアクセスしなかった人も含めた潜在的な利用者数 N を推定しようとした。

しかし日によって N の推定値が大きく変動しており、解釈しにくい結果だった。

そこで別の方法を考えた。

参考にしたのは『BUGS で学ぶ階層モデリング入門』の6章のはじめに出てくる、データ拡大(data augmentation)の例題である。

BUGSで学ぶ階層モデリング入門: 個体群のベイズ解析

BUGSで学ぶ階層モデリング入門: 個体群のベイズ解析

データ拡大

ゼロが切断(truncate)されたポアソン分布を考える。

ポアソン分布に従う試行が N 回行われたが、観測されるのは N 回のうち、ゼロでない値になった n 回の試行のみである。

このようなデータから元のポアソン分布のパラメータを正しく推定するために、データ拡大という方法がある。

データ拡大ではまず、データセット (y_1,y_2,\ldots,y_n) に十分大きな数 m 個の 0 を付け足す。

この m 個の 0 のうち、いくつかは本来のポアソン分布にしたがって発生したであろう 0 になる。

のこりのいくつかは、余分に付け足された 0 になる。

「いくつか」というのを混合比率 ρ でパラメタライズして確率分布で表すと、

\displaystyle \Pr(Y_{j}=0)=\rho +(1-\rho )e^{{-\lambda }}
 \displaystyle \Pr(Y_{j}=y_j)=(1-\rho ){\frac {\lambda ^{y_{i}}e^{-\lambda }}{y_{i}!}},~y_{i}\geq 1

である。

こうなると未観測の N を推定する問題が、混合分布のパラメータを推定するのと同じことになる。

y_i に 0 を付け足したデータを y_i^\ast ~ (i=1,\ldots,n+m) と表すことにする。

さらに z_iy_i^\ast \neq 0 ならば 1、y_i^\ast =0 ならば 0 の値をとる変数とすると、上記の混合分布の尤度関数は、

\displaystyle L=\prod_{i=1}^{n+m} \left\{(1-\rho )\frac{\lambda ^{y_{i}}e^{-\lambda }}{y_{i}!}\right\}^{z_i}  \left\{ \rho +(1-\rho )e^{-\lambda}\right\}^{1-z_i}

とかける。

これはゼロ過剰ポアソン(zero inflated Poisson)分布の尤度関数と同じである。

対数をとって整理すると、

\displaystyle \log L=n \log(1-\rho)+\log \lambda \sum y_j - \lambda n \\
\displaystyle~~~+ m \log\{\rho + (1-\rho) e^{-\lambda}\}

となる。

ウェブサイトのセッション数のモデル

あるウェブサイトの利用者が全部で N 人いるとして、N 人のがある一定時間内に y_i~(i=1,\ldots,N) 回そのウェブサイトを利用するとする。

ここでは、この y_iポアソン分布に従うと仮定する。

一定時間内に一回以上サイトを訪れた利用者の数を n とすると、 n に相当するのがユーザー数(ユニークユーザー数)である。

一定時間内にサイトを訪れるユーザーののべ人数に相当するのがセッション数(訪問)なので、セッション数は  \sum y_i に相当する。

y_i それぞれの値がわかっていなくても、尤度の計算には \sum_{i=1}^{n} y_inm がわかっていればよい。

シミュレーション

上記の状況をシミュレーションで再現して、パラメータを推定してみる。

わかっているのはある日のセッション数 \sum_{i=1}^{n} y_i と、ユーザー数 n_i である。

ポアソン分布のパラメータ  \lambda_i は日によって異なるとし、一様分布で生成した。

m は日によらず一定の値を与える。

ベイズ推定にあたって、  \lambda_i区間 [0, ∞) の一様分布、混合比率  \rho_i区間 [0, 1]の一様分布を仮定した。

Stan コードは以下のようになった。

data {
  int<lower=0> TT;
  int<lower=0> n[TT];
  int<lower=0> m;
  int<lower=0> ysum[TT];
}

parameters {
  real<lower=0,upper=1> rho[TT];
  real<lower=0> lambda[TT];
}

model {
  for(i in 1:TT){
   target += n[i]*log(1-rho[TT])+ysum[i]*log(lambda[i])-n[i]*lambda[i]+
    m*log(rho[TT]+(1-rho[TT])*exp(-lambda[i])); 
  }
}

generated quantities {
  int N[TT];
  for(i in 1:TT){
    N[i] = neg_binomial_rng(n[i],1-exp(-lambda[i])); 
  }
}

generated quantities ブロックでは、ゼロ切断ポアソン分布によるセッション数のモデル - 廿TT でやったのと同様、負の二項分布を使って N を生成している。

R のコードはこうなった。

library(rstan)
library(cowplot)
TT <-100
lambda=runif(TT,1,10)
N =100
ysum <- n <-numeric(TT)
for(i in 1:TT){
  y <- rpois(N,lambda[i])
  y2 <- y[y != 0]#
  ysum[i] <- sum(y2)
  n[i] <- length(y2)
}
plot(ysum/n)
m <-10000 #number of zeros
rho=(m-N)/(n+m) #
dat <- list(TT=TT,n=n,m=m,ysum=ysum)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
ZIP_aug <-stan_model("~/dropbox/ZIP_aug.stan")
fit <-sampling(ZIP_aug,dat)
#traceplot(fit,pars=c("rho[1]","lambda[1]"))

lambdahat <-get_posterior_mean(fit,"lambda")[,"mean-all chains"]
rhohat <-get_posterior_mean(fit,"rho")[,"mean-all chains"]
ggplot()+
  geom_point(aes(lambda,lambdahat))+
  geom_abline()

ex<-extract(fit)
CI1 <-t(apply(ex$N,2,quantile,c(0.025,0.975)))
Nhat1 <-get_posterior_mean(fit,"N1")[,"mean-all chains"]
Ndf <-as.data.frame(cbind(1:TT,Nhat1,CI1))
colnames(Ndf)<-c("date","N1","lower1","upper1")
ggplot(Ndf,aes(x=date))+
  geom_line(aes(y=N1),size=1)+
  geom_ribbon(aes(ymin=lower1,ymax=upper1),alpha=0.2)+
  geom_hline(yintercept=100,linetype=2)

まず下の図は横軸に  \lambda_i の真値、縦軸に事後平均推定値をとったプロットである。

f:id:abrahamcow:20160806214409p:plain

おおむね近い値が求まっていることがわかる。

次の図は日ごとに異なる N の推定値の折れ線グラフで、グレーの帯はその 95%ベイズ信頼区間である。

f:id:abrahamcow:20160806214615p:plain

真値の 100 付近にちらばり、ベイズ信頼区間は真値をカバーしていることがわかる。

ケーススタディ

シミュレーション結果から、この方針で安心して推定に使えそうだぞ、と思ったので、Google アナリティクスデータを用いて N の推定を行った。

library(RGA)
library(scales)
authorize()
prof <-list_profiles()

datGA <-get_ga(profileId = prof$id[1],
              start.date = "2016-07-01",
              end.date = "2016-07-31",
              dimensions = "ga:date",
              metrics = "ga:users,ga:sessions")
datGAlist <- list(TT=31,n=datGA$users,ysum=datGA$sessions,m=m)
fitGA <-sampling(ZIP_aug,datGAlist)
traceplot(fitGA,pars=c("rho[1]","lambda[1]"))

lambdahatGA <-get_posterior_mean(fitGA,"lambda")[,"mean-all chains"]
exGA <- extract(fitGA)
CIlambdaGA <-t(apply(exGA$lambda,2,quantile,c(0.025,0.975)))
lambdadfGA <-data.frame(datGA$date,lambdahatGA,CIlambdaGA)
lambdadfGA <-data.frame(datGA$date,NhatGA,CIGA)
colnames(lambdadfGA)<-c("date","lambda","lower","upper")
ggplot(lambdadfGA,aes(x=as.Date(date)))+
  geom_line(aes(y=lambda),size=1)+
  geom_ribbon(aes(ymin=lower,ymax=upper),alpha=0.2)+
  scale_x_date(date_labels ="%m/%d")+
  ylim(c(0,NA))+xlab("")+ylab("")


CIGA <-t(apply(exGA$N,2,quantile,c(0.025,0.975)))
NhatGA <-get_posterior_mean(fitGA,"N")[,"mean-all chains"]
NdfGA <-data.frame(datGA$date,NhatGA,CIGA)
colnames(NdfGA)<-c("date","N1","lower1","upper1")
ggplot(NdfGA,aes(x=as.Date(date)))+
  geom_line(aes(y=N1),size=1)+
  geom_ribbon(aes(ymin=lower1,ymax=upper1),alpha=0.2)+
  geom_line(aes(y=datGA$users),size=1,colour="royalblue")+
  scale_x_date(date_labels ="%m/%d")+
  ylim(c(0,NA))+xlab("")+ylab("")

日ごとのポアソン分布のレートパラメータの推定値と、その95%ベイズ信頼区間がこちら。

f:id:abrahamcow:20160806215011p:plain

日ごとの N の推定値と、その95%ベイズ信頼区間がこちら。青い折れ線グラフは観測されたユーザー数、オレンジがセッション数である。

f:id:abrahamcow:20160806215302p:plain

フューチャーワーク

  • これがなんのやくにたつか考える。
  • Nλ の変動はゆるやかだろう、という事前知識を織り込んだモデルの検討。
  • 推定結果が m にどの程度依存するか調べる。