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

廿TT

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

ゼロ切断ポアソン分布によるセッション数のモデル

R Google アナリティクス

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

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

一定時間内にサイトを訪れるユーザーののべ人数に相当するのがセッション数ですから、セッション数は  \sum k_i に相当します。

k_iポアソン分布に従うとすると、ポアソン分布のパラメータは、
ゼロ切断ポアソン分布のパラメータの最尤推定 - 廿TT でみたように、  \sum k_in があれば求まるので、セッション数とユーザー数から推定できます。

この場合ポアソン分布のパラメータは単位時間のユーザー一人あたりのセッション回数と解釈できます。

また、潜在的なウェブサイトの利用者数 N も推定することができます。

ただしこれも ゼロ切断ポアソン分布のパラメータの最尤推定 - 廿TT でみたように、Nプラグイン推定してやるとバイアスが出るようなので、パラメトリックブートストラップで求めることにしました。

f:id:abrahamcow:20160805050411p:plain

一枚目のグラフはポアソン分布のパラメータの推定値と、その95%信頼区間です。

ここではポアソン分布のパラメータは日によって異なるとしています。

信頼区間パラメトリックブートストラップによるもの(青)と正規近似によるもの(赤)の二種類を重ねてプロットしていますが、両者はほとんど重なっています。

f:id:abrahamcow:20160805052120p:plain

二枚目のグラフは潜在的なウェブサイトの利用者数 Nパラメトリックブートストラップによる推定値と95%信頼区間です。

下の青い線はその日のユーザー数を示しています。

ぼくとしては「ポアソン分布のパラメータは日々変化するけれども、潜在的N は日によってあまり変動しない」というような結果が出て欲しかったのですが、 N も結構変動してます。

こうなると解釈が難しい、というかモデルに妥当性がないのかもしれません。

以上の計算をやった R のコードを貼っておきます。

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

#Google アナリティクスデータの取得
dat1 <-get_ga(profileId = prof$id[1],
              start.date = "2016-07-01",
              end.date = "2016-07-31",
              dimensions = "ga:date",
              metrics = "ga:users,ga:sessions")

#尤度関数
logL0 <- function(ysum,n){
  function(lambda) {log(lambda)*ysum - n*lambda - n*log(1-exp(-lambda))}
}

#パラメトリックブートストラップを行う関数
simulater <- function(n,lambda,iter=1000){
  logL0 <- function(ysum,n){
    function(lambda) {log(lambda)*ysum - n*lambda - n*log(1-exp(-lambda))}
  }
  lambdahat <- N <- numeric(iter)
  for(i in 1:iter){
    y <- rpois(n*30,lambda) # 多めに出しておく
    N[i] <-which.max(cumsum(y!=0)==n)
    y2 <- y[y != 0][1:n]
    ysum <- sum(y2)
    logL <- logL0(ysum,n)
    opt1 <- optimise(logL,c(0,100),maximum = TRUE)
    lambdahat[i] <-opt1$maximum
  }
  cbind(lambdahat,N)
}

TT <- nrow(dat1)
lambdas <- numeric(TT)
simlist <- vector("list",TT)

#日ごとの計算
for(i in 1:TT){
  logL <- logL0(dat1$sessions[i],dat1$users[i])
  opt1 <- optimise(logL,c(0,100),maximum = TRUE)
  lambdas[i] <-opt1$maximum
  simlist[[i]] <- simulater(dat1$users[i],lambdas[i])
}

#ブートストラップ推定値
Nhat <-sapply(simlist,function(x)mean(x[,2]))
Nlower <-sapply(simlist,function(x)quantile(x[,2],0.025))
Nupper <-sapply(simlist,function(x)quantile(x[,2],0.975))

lamhat <-sapply(simlist,function(x)mean(x[,1]))
lamlower <-sapply(simlist,function(x)quantile(x[,1],0.025))
lamupper <-sapply(simlist,function(x)quantile(x[,1],0.975))

#正規近似による信頼区間を求める関数
ZTPCI <- function(lambda,n,level=0.95){
  NFI <- (n/(1-exp(-lambda)))*(1/lambda - exp(-lambda)/(1-exp(-lambda)))
  alpha = 1-level
  lower=qnorm(alpha/2,lambda,sqrt(1/NFI))
  upper=qnorm(1-alpha/2,lambda,sqrt(1/NFI))
  cbind(lower,upper)
}

CI = ZTPCI(lambdas,dat1$users)

estlam <- data.frame(date=dat1$date,lambda=lambdas,
                   upper1=lamupper,lower1=lamlower,
                   upper2=CI[,2],lower2=CI[,1])
#グラフ(一枚目)
ggplot(estlam ,aes(x=as.Date(date),y=lambda))+
  geom_point()+geom_line()+
  geom_ribbon(aes(ymin=lower1,ymax=upper1),alpha=0.2,fill="cornflowerblue")+
  geom_ribbon(aes(ymin=lower2,ymax=upper2),alpha=0.2,fill="tomato")+
  scale_x_date(date_labels ="%m/%d")+
  xlab("")+ylim(c(0,NA))

estN <- data.frame(date=dat1$date,population=Nhat,
                  upper=Nupper,lower=Nlower,users=dat1$users)
#グラフ(二枚目)
ggplot(estN,aes(x=as.Date(date),y=population))+
  geom_point()+geom_line()+
  geom_ribbon(aes(ymin=lower,ymax=upper),alpha=0.2)+
  geom_line(aes(y=users),colour="blue")+
  scale_x_date(date_labels ="%m/%d")+
  xlab("")+ylim(c(0,NA))

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