廿TT

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

Hawkes 過程を Web サイトへの訪問時刻に当てはめてみる

今日の川柳

点過程の時系列解析 (統計学One Point)

点過程の時系列解析 (統計学One Point)

『点過程の時系列解析』の例題です。

『点過程の時系列解析』はこれまでなかった本です。

統計的推定とモデリングの視点から書かれた点過程の入門書は日本初じゃないでしょうか。

p.148 の Hawkes 過程の尤度関数(8.54式)はたぶん誤植で(-bが抜けてる)指数関数カーネルを用いた Hawkes 過程の対数尤度関数は

 \sum_{i} \log \{\mu +\sum_{j < i} ab \exp [ - b (t_i-t_j) ]\}- \{ \mu T \sum_{i} a (1-\exp[-b(t_i-t_j)]) \}

です。

この尤度を使って Google アナリティクスのデータに Hawkes 過程を当てはめてみます。

f:id:abrahamcow:20190618204336p:plain

f:id:abrahamcow:20190618204442p:plain

intensityは単位時間(ここでは1分)あたりのセッション数、countは累積のセッション数です。

ただ当てはめてみただけで別にだからなんだというのはないんですけど、Rのコードを貼ります。

library("googleAnalyticsR")
library("tidyverse")
ga_auth()
account_list <- ga_account_list()
ga_id <- account_list[3,'viewId']

dat_ga <- google_analytics(ga_id,
                             date_range = c("2019-06-13","2019-06-13"),
                             metrics = c("sessions"),
                             dimensions = c("date","hour","minute"))

dat_ga2 <-dat_ga %>% 
  mutate(time=as.POSIXct(paste(date,hour,minute),format="%Y-%m-%d %H %M")) %>%
  arrange(time) %>% 
  mutate(count=row_number())
theme_set(theme_bw(16))
ggplot(dat_ga2,aes(x=time))+
  geom_step(aes(y=count))+
  geom_rug()
ti <- as.numeric(dat_ga2$time-as.POSIXct("2019-06-13 00:00:00"))
Tau <- 24*60
ll <- function(par){
  mu <- exp(par[1])
  b <- exp(par[2])
  a <- exp(par[3])
  sum(log(mu+sapply(ti, function(t)sum(a*b*exp(-b*(t-ti[ti<t])))))) - 
    (mu*Tau+sum(a*(1-exp(-b*(Tau-ti)))))
}
opt <- optim(rep(0,3),ll,control = list(fnscale=-1))

mu <- exp(opt$par[1])
b <- exp(opt$par[2])
a <- exp(opt$par[3])
xv <- seq(0,Tau,by=1)
lv <- mu+sapply(xv, function(t)sum(a*b*exp(-b*(t-ti[ti<t]))))
Lv <- sapply(xv, function(t){mu*t+sum(a*(1-exp(-b*(t-ti[ti<t]))))})
df <- data.frame(x=xv,intensity=lv,L=Lv)
ggplot(df)+
  geom_line(aes(x=x*60+as.POSIXct("2019-06-13 00:00:00"),y=intensity))+
  geom_rug(data=dat_ga2,aes(x=time))+
  xlab("time")
ggplot(dat_ga2,aes(x=time,y=count))+
  geom_step(size=1)+
  geom_line(data=df,aes(x=x*60+as.POSIXct("2019-06-13 00:00:00"),y=L),colour="royalblue",size=1)