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

廿TT

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

[rstan]横浜駅SFの移入と離脱のモデル

R Stan

横浜駅SF (カドカワBOOKS)

横浜駅SF (カドカワBOOKS)

経緯


はじめに:とりあえず棒グラフ

まずは rvest で 横浜駅SF(柞刈湯葉)のアクセス数 - カクヨム からPV数のデータを抜き出してきます。

library(rvest)
library(pipeR)
read1 <- read_html("https://kakuyomu.jp/works/4852201425154905871/accesses")
tab1 <-html_table(read1)
tab2 <-strsplit(tab1[[1]]$X2,"\n")
tab3 <-do.call("rbind",tab2)
PVchar <-gsub(",","",tab3[,1])
dat <-data.frame(number=1:length(PVchar),title=tab1[[1]]$X1,PV= as.integer(PVchar))

横浜駅SF、各話のPV数は下図の通りです。

f:id:abrahamcow:20160512033807p:plain

装飾的な要素を廃した棒グラフは次の図です。

f:id:abrahamcow:20160512034001p:plain

各エピソードに上から順に通し番号をふっています。

モデルと推定

閉じた形の関数で説明するのは諦め、確率モデルを考えました。

PV数は、

 y_i = r_i + p_i y_{i-1}

すなわち、

エピソード i のPV 数(y_i)= 移入数(r_i)+ 前エピソードからの遷移(p_i y_{i-1}

という形に分解できると仮定しました。

移入数はそのエピソードから閲覧を開始するPVの数です。前エピソードからの遷移は前エピソードのPV数に遷移する割合(0 \le p_i \le 1)を乗じた形で表すことにします。

PV数の分布はポアソン分布を仮定しました。

となりあったエピソードは条件が近いだろう、ということから、\log r_i の事前分布は平均 \log r_{i-1}正規分布{\rm logit}( p_i) の事前分布は平均 {\rm logit}( p_{i-1})正規分布としました。

無情報事前分布として、正規分布の分散は幅の広い一様分布にしました。

モデルを記述した Stan のコードはこれです。

data {
  int<lower=1> N_site;
  int<lower=0> Y[N_site];
}

parameters {
  real r[N_site];
  real p0[N_site];
  real<lower=0,upper=1e+7> s_r;
  real<lower=0,upper=1e+7> s_p;
}

model {
  Y[1] ~ poisson_log(r[1]);
  for (j in 2:N_site){
   r[j] ~ normal(r[j-1], s_r);
   p0[j] ~ normal(p0[j-1], s_p);
   Y[j] ~ poisson(exp(r[j])+inv_logit(p0[j])*Y[j-1]);
  }
}

generated quantities {
  real<lower=0, upper=1> p[N_site];
  real<lower=0> expr[N_site];
  real<lower=0> Y_mean[N_site];
  for (j in 1:N_site){
    expr[j] <- exp(r[j]);
    p[j] <- inv_logit(p0[j]);
  }
  Y_mean[1] <- exp(r[1]);
  for (j in 2:N_site){
    Y_mean[j] <- exp(r[j])+inv_logit(p0[j])*Y[j-1];
  }
}

R から Stan を操作します。

dat4stan <-list(Y=dat$PV,N_site=nrow(dat))
library(rstan)
model1 <-stan_model("yokohamamodel.stan")
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
fit1 <-sampling(model1,dat4stan,iter=6000,warmup=5000,
                control=list(adapt_delta=0.9,max_treedepth=15))

#traceplot(fit1,c("s_r","s_p"))
#traceplot(fit1,"p0")
#traceplot(fit1,"r")
Y_mean <-get_posterior_mean(fit1,"Y_mean")[,5]
expr <-get_posterior_mean(fit1,"expr")[,5]
p <-get_posterior_mean(fit1,"p")[,5]

expr <-get_posterior_mean(fit1,"expr")[,5]
ggplot(dat)+
  geom_bar(aes(x=number,y=PV),stat = "identity",fill="grey60")+
  geom_line(aes(x=number,Y_mean),colour="blue",size=1)+
  geom_line(aes(x=number,expr),colour="red",size=1)+
  scale_y_continuous(labels=comma)

当てはめた結果は下図のようになりました。

f:id:abrahamcow:20160512035752p:plain

赤い線は移入数、青い線は移入数と前エピソードからの遷移をあわせたPV数の推定値です。

考察

係数 p_i は全エピソードから遷移してくる割合を示しているため、1-p_i は離脱率(そのページで読み進めるのをやめる割合)と解釈できます。

各エピソードの離脱率とその95%信用区間は下図のようになりました。

f:id:abrahamcow:20160512040137p:plain

ex <-extract(fit1)
upper_p <-apply(1-ex$p,2,function(x)quantile(x,0.975))
lower_p <-apply(1-ex$p,2,function(x)quantile(x,0.025))

ggplot(dat)+
  geom_bar(aes(x=reorder(title,-number),y=1-p),stat = "identity",alpha=0.6)+
  geom_errorbar(aes(x=rev(number),ymin=lower_p,ymax=upper_p),colour="red")+
  coord_flip()+
  xlab("")+ylab("離脱率")+
  scale_y_continuous(labels=percent)

序盤は離脱率が高く徐々に離脱率は下がっていくのですが、12話、16話、21話など、ちょこっと離脱率が高くなっているエピソードも見受けられます。

そのあたりがユーザーが疲れてくるタイミングなのかもしれません。

各エピソードの移入数とその95%信用区間は下図のようになりました。

f:id:abrahamcow:20160512040857p:plain

upper_r <-apply(ex$expr,2,function(x)quantile(x,0.975))
lower_r <-apply(ex$expr,2,function(x)quantile(x,0.025))

ggplot(dat)+
  geom_bar(aes(x=reorder(title,-number),y=expr),stat = "identity",alpha=0.6)+
  geom_errorbar(aes(x=rev(number),ymin=lower_r,ymax=upper_r),colour="red")+
  coord_flip()+
  xlab("")+ylab("移入数")+
  scale_y_continuous(labels=comma)

やはり最終話から読みはじめるユーザーがけっこういるようです。

アマゾンアフィリエイトのコーナー

地球礁 (河出文庫)

地球礁 (河出文庫)

ぼくの好きなSF『地球礁』です。

子どもたちが韻を踏んだ詩で人を殺す話です。

訳文もリズミカルで読みやすいです。

買ってください。よろしくお願いします。