不倫報道の間隔を推定する

2016年は不倫報道が多いですね！全くどうでもいいんですが，これほど多いとデータとしてみなし，なんかの分析をしたくなる。

そこでデータをつくって分析してみます。

dataset <- as.data.frame(
matrix(c("ベッキー","2016-01-07",
         "狩野英孝","2016-02-03",
         "宮崎謙介議員","2016-02-10",
         "桂文枝","2016-03-16",
         "石井竜也","2016-03-17",
         "乙武洋匡","2016-03-23",
         "とにかく明るい安村","2016-04-01",
         "NMB48木下春奈","2016-05-08",
         "ファンキー加藤","2016-06-07",
         "三遊亭円楽","2016-06-10",
         "荻上チキ","2016-07-06",
         "小倉優子の夫","2016-08-02",
         "中村橋之助","2016-09-15",
         "ジュビロ磐田の藤田義明","2016-10-04",
         "浦沢直樹","2016-10-10",
         "矢島悠子アナウンサー","2016-10-13",
         "ラモス瑠偉","2016-10-18",
         "小渕健太郎","2016-10-20"
),ncol=2,byrow=T))
names(dataset) <- c("who","date")

# POSIXct型に変換
furin <- as.POSIXct(dataset$date)
# 不倫報道が生じるまでの日数を計算
span <- c()
for(i in 1:(length(furin)-1)){
  span[i] <- difftime(furin[i+1],furin[i])
}

「報道が生じる」のがポアソン分布，ポアソン分布で次の観測が生じるまでの感覚は指数分布に従うと考えられるので，これをモデル化。

# Stanで予測する
library(rstan)
options(mc.cores = parallel::detectCores())

# 

stancode <- '
data{
  int<lower=0> N;
  real<lower=0> X[N];
}
parameters{
  real<lower=0> lambda;
}
model{
  target += exponential_lpdf(X|lambda);
}
generated quantities{
  real<lower=0> mu;
  mu = inv(lambda);
}
'

# 推定
datastan <- list(N=length(span),X=span)
model <- stan_model(model_code=stancode)
fit <- sampling(model,data=datastan)

結果のmuの値が示すのが，次の不倫報道まで待つべき日数になります。

print(fit)

## Inference for Stan model: da8318251705d5281abfd8435af920f0.
## 4 chains, each with iter=2000; warmup=1000; thin=1; 
## post-warmup draws per chain=1000, total post-warmup draws=4000.
## 
##          mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
## lambda   0.06    0.00 0.01   0.04   0.05   0.06   0.07   0.09  2859    1
## mu      16.96    0.08 4.24  10.57  13.94  16.33  19.19  26.74  2806    1
## lp__   -68.34    0.02 0.72 -70.42 -68.51 -68.07 -67.89 -67.84  2252    1
## 
## Samples were drawn using NUTS(diag_e) at Fri Oct 21 16:04:18 2016.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

もちろんこれは，不倫「報道」までの日数予測であって，かなり人為的な要素があることです。予測の精度を検証しようってな話ではありませんが，占いによる予想よりはよっぽどマシでしょう，という笑い話としてお納めください。