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

廿TT

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

状態空間非定常ポアソン(NHPP using Stan)

R 確率過程 Stan

ポアソン過程は再発事象のモデルとしてよく使われる。

ポアソン過程ではイベントが観測された時刻を T_t (t=1,\ldots,N) とすると、イベントの生起間隔  \Delta T_t = T_t - T_{t-1} は独立にパラメータ λ の指数分布に従う。

ポアソン過程の拡張としてパラメータλ が時間に依存して変化する非定常ポアソン過程というのが考えられる。

今回は以下のように一回差分のトレンドで時間変化を表現する。

 \mu_t \sim {\rm Normal}(\mu_{t-1},\sigma^2),~t=2,\ldots,N
 \Delta T_t \sim {\rm Exponential}(\exp(\mu_t)) ~ t=2,\ldots,N

\mu_1\sigma には幅の広い一様分布を仮定する。

Stan コードはこうなりました。

//NHPP.stan
data {
  int<lower=1> n;
  real<lower=0> D[n];
}
parameters {
  vector[n] loglambda;
  real<lower=0>sigma;
}
model {
  loglambda[2:n] ~ normal(loglambda[1:(n-1)],sigma);
  D ~ exponential(exp(loglambda));
}

R の boot パッケージには炭鉱災害の発生を記録した coal データというのが入っている。

data("coal",package = "boot")
p1 <-ggplot(coal)+
  geom_step(aes(x=date,y=1:length(date)))+
  ylab("Cumulative recurrences")
print(p1)

f:id:abrahamcow:20161126085429p:plain

これに上記の NHPP(Non-Homogeneous Poisson Process; 非定常ポアソン過程)を当てはめる。

rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
NHPPmodel <-stan_model("~/Documents/NHPP.stan")
dat <-list(D=diff(ti),n=n-1)
NHPPfit <-sampling(NHPPmodel,dat,iter=10000)
lambdahat <- exp(get_posterior_mean(NHPPfit,pars="loglambda")[,"mean-all chains"])
exNHPP <-rstan::extract(NHPPfit)
lower <-exp(apply(exNHPP$loglambda, 2, quantile,0.2))
upper <-exp(apply(exNHPP$loglambda, 2, quantile,0.8))
p1+geom_line(aes(x=cumsum(c(date[1],1/lambdahat)),y=1:length(date)),colour="royalblue",size=1)+
  geom_line(aes(x=cumsum(c(date[1],1/lower)),y=1:length(date)),colour="royalblue",size=1,linetype=2)+
  geom_line(aes(x=cumsum(c(date[1],1/upper)),y=1:length(date)),colour="royalblue",size=1,linetype=2)

f:id:abrahamcow:20161126085730p:plain

点線は60%のベイズ信頼区間です。

これがなんの役に立つのかは後で考える。

2回差分のトレンドでもやってみました。

data {
  int<lower=1> n;
  real<lower=0> D[n];
}
parameters {
  vector[n] loglambda;
  real<lower=0>sigma;
}
model {
  loglambda[3:n] ~ normal(2*loglambda[2:(n-1)]-loglambda[1:(n-2)],sigma);
  D ~ exponential(exp(loglambda));
  sigma ~ exponential(0.1);
}

sigma が収束しにくかったので指数分布を入れて縛ってます。

f:id:abrahamcow:20161126090235p:plain

そのせいかベイズ信頼区間の幅はちょっと狭くなった。

グラフを見た限りでは、インテンシティ(loglambda)の経時的な変化はどちらでも大差ない感じ。

StanとRでベイズ統計モデリング (Wonderful R)

StanとRでベイズ統計モデリング (Wonderful R)

abrahamcow.hatenablog.com

(R+Google アナリティクス)エラーバーで信頼下限をプロット

Google アナリティクス R グラフ 乱数

場面設定

当サイトは女性の訪問者が少ないので、女性の訪問者を増やしたいと思っている。

サイト作りの参考にしようと女性の新規訪問の割合が多いランディングページをリストアップしたい。

そこで新規セッション率で降順にソートをかけると、セッション数10〜20くらいのページばかりが上の方にきた。

分母の数が小さいときに割り算値を額面通りに受け取ってよいものか。

そういうときには区間推定をして、信頼下限で判定してやろう。

区間推定では小さく見積もればこれくらい、大きく見積もればこれくらい、という値を分母の数に応じて出すことができる。

実践

信頼度95%の信頼区間を棒で、新規率を点で表したプロットが以下である。

f:id:abrahamcow:20161126060421p:plain

点の大きさでセッション数も合わせて表現した。

セッション数が多くなるとき、信頼区間の幅は狭くなることがわかる。

#パッケージをまとめて読み込み
library(RGA)
library(dplyr)
library(tidyr)
library(cowplot)
#パッケージ読み込み終わり

#データ取得
authorize()
prof <-list_profiles()
dat_ga <-get_ga(profileId = prof$id[1],
                start.date = "2016-10-01",
                end.date = "2016-10-31",
                dimensions = "ga:landingPagePath,ga:userGender",
                metrics = "ga:sessions,ga:newUsers")
#データ取得終わり

#信頼区間を計算
CI <-t(mapply(function(x,n)binom.test(x=x,n=n)$conf.int, dat_ga$newUsers,dat_ga$sessions))
LPdf <- data.frame(dat_ga,CI) %>% 
  rename(lower=X1,upper=X2) %>% 
  mutate(percent_new=newUsers/sessions)

#女性の新規率の信頼下限が0.8以上のページを抜き出す
LPdf085 <-filter(LPdf,landingPagePath %in% landingPagePath[userGender=="female"& lower>=0.8])

#以下プロット
theme_set(theme_cowplot(font_family = "HiraKakuProN-W3"))
ggplot(LPdf085,aes(x=landingPagePath,y=percent_new,colour=userGender))+
  geom_pointrange(aes(ymin=lower,ymax=upper),position=position_dodge(width = 0.3))+
  geom_point(aes(size=sessions),position=position_dodge(width = 0.3))+
  scale_y_continuous(label=scales::percent)+
  coord_flip()+ylab("新規率")+xlab("")+
  geom_hline(yintercept=0.9,linetype=2)

新規率以外の比率、直帰率やコンバージョン率でも同じプロットが使える。

補足:シミュレーション

信頼度95%の信頼区間は、100回推定したら内95回くらいはそのなかに真のパラメータを含むだろう。

シミュレーションしてみる。

ここではセッション数がポアソン分布に従ってランダムに変動するとして、二項分布の成功確率パラメータの信頼区間を繰り返し構成した。

N <-rpois(10000,100)
x <- rbinom(10000,N,0.8)
CIsim <-t(mapply(function(x,n)binom.test(x=x,n=n)$conf.int, x,N))
cat("coverage probability is" ,round(mean(CIsim[,1]<=0.8 & CIsim[,2]>=0.8),2))

Coverage probability(信頼区間がその範囲内に真値を含む割合)は0.96と95%に近い値になった。

品質管理のための統計手法 (日経文庫)

品質管理のための統計手法 (日経文庫)