Hatena::ブログ(Diary)

驚異のアニヲタ社会復帰への道

Prima Project

2018-06-25

MikuHatsune2018-06-25

W杯の試合観戦中にトイレはいついくべきか

こんなツイートを観測した。

 

ハーフタイムに水道使用料が増えているのがわかる。

試合中に離席すると一番盛り上がる得点シーンを見逃してしまうため、試合中はなかなかトイレやお風呂にいけない。

というわけで試合中に離席するにはどの時間帯が一番よいかを調べる。

 

高校サッカーの点差を解析したときと同様に、過去のW杯の試合結果から得点が入った時刻を取得する。ここで、1930年から2014年大会までの20大会(1942年と1946年は中止)について、836試合あり、得点シーンは2373だった(wikiFIFA W杯のページをパースしたため、本当にそうなのかはわからない)。

1970年大会が6ゴール取得できてなかったようである。

1930 1934 1938 1950 1954 1958 1962 1966 1970 1974 1978 1982 1986 1990 1994 1998 2002 2006 2010 2014 
  70   70   84   88  140  126   89   89   89   97  102  146  132  115  141  171  161  147  145  171 

 

いつ得点が入るかをベタに分布をとると、高校サッカー選手権とほぼ同様で、どの時間帯でもほぼ一様な感じである。ただし、アディショナルタイムはすべて45分もしくは90分に換算しており、後半90分はほかの時間帯に比べて突出して得点が入っているので、試合終了間際は離席せずに見届けたほうがよい。

前半では全体の44%、後半では56% の得点が分布するので、どちらかというと前半のほうに離席するのがよい。

f:id:MikuHatsune:20180625210132p:image

 

ある時間帯に得点が入ってから、次の得点が入るまでの時間の分布は、ガンマ分布のようになる。このため、ガンマ分布で推定すると、パラメータ(1.266, 0.0554) をもつ以下のようなガンマ分布になる。

得点が入ってから17分間で、50% の確率で次の点がはいるようなので、得点が入った直後は油断しないほうがよい。

f:id:MikuHatsune:20180625210136p:image

 

ある時間帯にN点差がついているときに、そのままリードして試合終了する確率も選手権の解析と同様にしてみたが、ほとんど高校サッカーと同じ結果になった。試合終了間際まで1点リードしているとき、劇的ゴールで追いつくのは2.4% しかない(89分時点で1点リードしてそのまま90分=試合終了までリードを保つ確率が97.6%)。

勝ったな(確信)と思って風呂にはいったり寝たりするのは、90%の確信度でいくならば1点差なら後半75分あたり、2点差なら前半30分あたりでよい。

f:id:MikuHatsune:20180625210139p:image

 

今回のデータ取りはstringr を使ってR 内で完結させた。

# データとり
# wiki から
# W杯のトップページより、各グループステージ、決勝トーナメントの
# ページからデータをパースするほうが、フォーマットが整っていて効率がよい

url <- NULL
y1 <- seq(2014, 1998, -4)
lab <- c(sprintf("_Group_%s", LETTERS[1:8]), "_knockout_stage")
url <- c(url, mapply(function(z) sprintf("https://en.wikipedia.org/wiki/%d_FIFA_World_Cup%s", z, lab), y1))

y1 <- seq(1994, 1986, -4)
lab <- c(sprintf("_Group_%s", LETTERS[1:6]), "_knockout_stage")
url <- c(url, mapply(function(z) sprintf("https://en.wikipedia.org/wiki/%d_FIFA_World_Cup%s", z, lab), y1))

y1 <- 1982
lab <- c(sprintf("_Group_%s", c(1:6, LETTERS[1:4])), "_knockout_stage")
url <- c(url, mapply(function(z) sprintf("https://en.wikipedia.org/wiki/%d_FIFA_World_Cup%s", z, lab), y1))

y1 <- c(1978, 1974)
lab <- c(sprintf("_Group_%s", c(1:4, LETTERS[1:2])), "_knockout_stage")
url <- c(url, mapply(function(z) sprintf("https://en.wikipedia.org/wiki/%d_FIFA_World_Cup%s", z, lab), y1))

y1 <- c(seq(1970, 1954, -4), 1930)
lab <- c(sprintf("_Group_%s", 1:4), "_knockout_stage")
url <- c(url, mapply(function(z) sprintf("https://en.wikipedia.org/wiki/%d_FIFA_World_Cup%s", z, lab), y1))

y1 <- 1950
lab <- c(sprintf("_Group_%s", 1:4), "_final_round")
url <- c(url, mapply(function(z) sprintf("https://en.wikipedia.org/wiki/%d_FIFA_World_Cup%s", z, lab), y1))

y1 <- c(1938, 1934)
lab <- c(sprintf("_Group_%s", 1:4), "_final_tournament")
url <- c(url, mapply(function(z) sprintf("https://en.wikipedia.org/wiki/%d_FIFA_World_Cup%s", z, lab), y1))
write(url, "list.txt")

# パース
library(stringr)
fi <- list.files(pattern="FIFA")
Res <- NULL
g <- 0
for(f in fi){
  txt <- readLines(f)
  flag <- c("score"=0, "right"=1)
  res <- NULL
  for(tmp in txt){
    #if(str_detect(tmp, "<th style.*\\d+&#8211;\\d+.*</th>")){
    if(str_detect(tmp, "<th style.*&#8211;.*</th>")){ # 1986年のグループステージが半角スペースがあっておかしい
      flag["score"] <- 1
      goaltime <- vector("list", 2)
    }
    if(flag["score"] == 1){
      if(str_detect(tmp, "[\\d+]+'")){
        goaltime[[ flag["right"] ]] <- c(goaltime[[ flag["right"] ]], gsub("'", "", str_extract_all(tmp, "[\\d+]+'")[[1]]))
      }
      if(str_detect(tmp, "Report")){
        flag["right"] <- 2
      }
      if(str_detect(tmp, "</table>")){
        res <- c(res, list(goaltime))
        flag <- c("score"=0, "right"=1)
      }
    }
  }
  year <- as.numeric(str_extract(f, "\\d+"))
  for(i in seq(res)){
    g <- g + 1
    for(j in seq(res[[i]])){
      for(k in res[[i]][[j]]){
        l <- as.numeric(strsplit(k, "\\+")[[1]])
        hoge <- c(year, g, j, l, rep(0, 2-length(l)))
        Res <- rbind(Res, hoge)
      }
    }
  }
}
colnames(Res) <- c("year", "gameID", "HA", "time", "extra")
# 解析
dat <- read.table("score.txt", header=TRUE)
# 試合数の確認
Ngame <- c(18, 17, 18, 22, 26, 35, rep(32, 3), 38, 38, rep(52, 4), rep(64, 5))
mapply(function(z) length(unique(z$gameID)), split(dat, dat$year))

# 得点時間分布
sb <- subset(dat, time <= 90)
t1 <- 1:45
t2 <- 46:90
tab <- table(factor(sb$time, c(t1, t2)))/nrow(sb) * 100

cols <- c("red", "green")
par(mar=c(5, 5, 2, 2), cex.lab=1.6, cex.main=2)
b <- barplot(tab, col=c(mapply(rep, cols, each=sapply(list(t1, t2), length))), las=1, main="得点時間分布", ylab="頻度[%]", ylim=c(0, 3))
mtext("前半(分)", 1, line=3, at=mean(t1))
mtext("後半(分)", 1, line=3, at=mean(t2)+15)

# 得点間隔
sb <- dat
Tgoal <- lapply(split(sb$time, sb$gameID), sort)
difftime <- lapply(Tgoal, function(z) tail(c(0, z), -1) - head(c(0, z), -1))
difftime <- unlist(difftime)

library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

# ガンマ分布による得点時間間隔
code <- "
data{
  int N;
  vector<lower=0>[N] Time;
}
parameters{
  real<lower=0, upper=50> p[2];
}
model{
  Time ~ gamma(p[1], p[2]);
}
"

m <- stan_model(model_code=code)
standata <- list(N=length(difftime), Time=difftime)
fit <- sampling(m, data=standata, iter=6000, warmup=1500, seed=1234)
ex <- extract(fit)

tab <- table(factor(difftime, 0:max(difftime)))/length(difftime) * 100
b <- barplot(tab, las=1)
ps <- apply(ex$p, 2, median)
x0 <- seq(0, max(difftime), length=1000)
y0 <- dgamma(x0, ps[1], ps[2])

alpha <- c(0.25, 0.5, 0.75, 0.9)
d0 <- c(5, 0.2)
par(mar=c(5, 5, 2, 2), cex.lab=1.6, cex.main=2)
plot(tab, xlab="前の得点からの経過時間(分)", ylab="頻度[%]", las=1)
lines(x0, y0*100, col=2, lwd=3)
for(i in seq(alpha)){
  x1 <- qgamma(alpha[i], ps[1], ps[2])
  y1 <- dgamma(x1, ps[1], ps[2])*100
  arrows(x1+d0[1], y1+d0[2], x1, y1, length=0.1, lwd=3, col=4)
  points(x1, y1, pch=16, col=4)
  text(x1+d0[1], y1+d0[2], sprintf("%2d %s (%.1f min)", alpha[i]*100, "%", x1), pos=4)
  title("次の得点の時間の分布")
}

# 得点差勝利確率
sb <- subset(dat, time <= 90)
mat <- mat0 <- mat1 <- mat2 <- matrix(0, max(dat$gameID), 90)
idx <- c(1, -1)
for(i in unique(sb$gameID)){
  tmp <- subset(sb, gameID==i)
  tmp <- tmp[order(tmp$time),]
  for(j in 1:nrow(tmp)){
    mat[i, tmp$time[j] ] <- mat[i, tmp$time[j]] + idx[tmp$HA[j]]
  }
}

mat0[,1] <- mat[,1]
for(j in 2:ncol(mat)){
  mat0[,j] <- mat0[, j-1] + mat[, j]
}

# 1点差以上
mat1[abs(mat0) > 0] <- 1
prob1 <- prob2 <- rep(0, 89)
for(j in 1:89){
  j1 <- mat1[, j] > 0
  prob1[j] <- mean(sapply(apply(mat1[j1, j:90], 1, unique), length) <= 1)
}

mat2 <- abs(mat0)
for(j in 1:89){
  j1 <- mat2[, j] > 1
  prob2[j] <- mean(apply(mat2[j1, j:90] > 0, 1, all))
}

p <- cbind(prob1, prob2)
par(mar=c(5, 5, 2, 4), cex.lab=1.6, cex.main=2, las=1)
matplot(p, xlab="時間", ylab="勝利確率", col=cols, pch=16, xaxp=c(0, 90, 6))
abline(v=45, lty=3)
legend("bottomright", legend=sprintf("%d 点差", 2:1), col=rev(cols), pch=16, cex=2)
text(par()$usr[2], tail(p[,1], 1), sprintf("%.1f %s", tail(p[,1], 1)*100, "%"), pos=4, xpd=TRUE)

2018-06-14

MikuHatsune2018-06-14

決勝トーナメントに向けて初戦が大事というが初戦はどれくらい大事なのか

2018FIFAワールドカップが始まった。日本の決勝トーナメント進出は初戦に勝利できるかにかかっている! とかなんとかよく聞くが、実際に初戦に勝利するのはどれだけ大事かは定量的に言われていない気がする。

2002年と2010年に決勝トーナメントに進出したので8年周期のジンクスがある! みたいに言っていたTV があったが何を言っているのかさっぱりわからなかった。

 

というわけで、32チームを4チーム8グループに分けて、上位2チームが決勝トーナメントに進出するような仕様になった1998年から2014年まで5大会分の延べ160チームの勝敗表から、初戦に勝利すると決勝トーナメント進出にどれだけ有利になるかを調べる。

初戦に勝つのをW、決勝トーナメントに進出するのをT とすると、初戦に勝って決勝トーナメントに進出するのはP(T|W) である。

地味に数えればよい。行は初戦での勝点、つまり上から負け、引き分け、勝ち、であり、列は決勝トーナメントに進出したからのFalse/True である。

初戦に勝ったうえで、決勝トーナメントに進出しているかは、3の行を見ればいいから、51/(9+51) の0.85 の確率で決勝トーナメントに進出している。逆に、初戦に勝っても0.15 は決勝トーナメントに進出できていない。

     0  1
  0 53  7
  1 18 22
  3  9 51

逆に、初戦で負けた場合に決勝トーナメントに進出できるかどうかは、0 の行を見ればよい。初戦で負けると 7/(53+7) の0.12 程度しか決勝トーナメントに進出できない。

 

勝点をどれくらいとれば決勝トーナメントに進出できるかを見てみると、勝点3 でも決勝トーナメントに進出したことが1回だけある(1998年Bグループのチリ)。

基本的には勝点5、つまり1勝2分ならば他のチームの星の取り合いを考慮しても決勝トーナメントには進出できる。勝点4の1勝1分1敗パターンは微妙で、決勝トーナメントに進出できるかどうかは5分5分である。というのも、このときは他のチームの勝敗パターンが同様に勝点4のチームがいることが多くて、2位と3位になっていることが多い。

points  0  1
     0 12  0
     1 24  0
     2  8  0
     3 22  1
     4 14 15
     5  0 15
     6  0 14
     7  0 21
     9  0 14

ROC をしたけどあまり意味がない。

f:id:MikuHatsune:20180614142431p:image

 

dat <- read.table(text=txt, header=TRUE)
# 初戦の勝敗と決勝トーナメント進出
tab <- table(dat$g1, dat$final)

library(pROC)
points <- rowSums(dat[, grep("g[1-3]", colnames(dat))])
# 勝点と決勝トーナメントの進出
table(points, dat$final)

r <- roc(dat$final, points)
thres <- coords(r, c(0:7, 9), "threshold")

lot(r, lwd=3)
points(thres[2,], thres[3,], pch=15, cex=1.5)
text(thres[2,1:4], thres[3,1:4], colnames(thres)[1:4], adj=c(NA, 1.8), cex=1.5)
text(thres[2,5:9], thres[3,5:9], colnames(thres)[5:9], adj=c(-1, 1.8), cex=1.5)

txt <- "
year group teamID g1 g2 g3 final
1998 A Brazil 3 3 0 1
1998 A Norway 1 1 3 1
1998 A Morocco 1 0 3 0
1998 A Scotland 0 1 0 0
1998 B Italy 1 3 3 1
1998 B Chile 1 1 1 1
1998 B Austria 1 1 0 0
1998 B Cameroon 1 0 1 0
1998 C France 3 3 3 1
1998 C Denmark 3 1 0 1
1998 C South_Africa 0 1 1 0
1998 C Saudi_Arabia 0 0 1 0
1998 D Nigeria 3 3 0 1
1998 D Paraguay 1 1 3 1
1998 D Spain 0 1 3 0
1998 D Bulgaria 1 0 0 0
1998 E Netherlands 1 3 1 1
1998 E Mexico 3 1 1 1
1998 E Melgium 1 1 1 0
1998 E South_Korea 0 0 1 0
1998 F Germany 3 1 3 1
1998 F Yugoslavia 3 1 3 1
1998 F Iran 0 3 0 0
1998 F United_States 0 0 0 0
1998 G Romania 3 3 1 1
1998 G England 3 0 3 1
1998 G Colombia 0 3 0 0
1998 G Tunisia 0 0 1 0
1998 H Argentina 3 3 3 1
1998 H Croatia 3 3 0 1
1998 H Jamica 0 0 3 0
1998 H Japan 0 0 1 0
2002 A Denmark 3 1 3 1
2002 A Senetal 3 1 1 1
2002 A Uruguay 0 1 1 0
2002 A France 0 1 0 0
2002 B Spain 3 3 3 1
2002 B Paraguay 1 0 3 1
2002 B South_Africa 1 3 0 0
2002 B Slovenia 0 0 0 0
2002 C Brazil 3 3 3 1
2002 C Turkey 0 1 3 1
2002 C Costa_Rica 3 1 0 0
2002 C China_PR 0 0 0 0
2002 D South_Korea 3 1 3 1
2002 D United_States 3 1 0 1
2002 D Portugal 0 3 0 0
2002 D Poland 0 0 3 0
2002 E Germany 3 1 3 1
2002 E Ireland 1 1 3 1
2002 E Cameroon 1 3 0 0
2002 E Saudi_Arabia 0 0 0 0
2002 F Sweden 1 3 1 1 
2002 F England 1 3 1 1
2002 F Argentina 3 0 1 0
2002 F Nigeria 0 0 1 0
2002 G Mexico 3 3 1 1
2002 G Italy 3 0 1 1
2002 G Croatia 0 3 0 0
2002 G Ecuador 0 0 3 0
2002 H Japan 1 3 3 1
2002 H Belgium 1 1 3 1
2002 H Russia 3 0 0 0
2002 H Tunisia 0 1 0 0
2006 A Germany 3 3 3 1
2006 A Ecuador 3 3 0 1
2006 A Poland 0 0 3 0
2006 A Costa_Rica 0 0 0 0
2006 B England 3 3 1 1
2006 B Sweden 1 3 1 1
2006 B Paraguay 0 0 3 0
2006 B Trinidad_and_Tobago 1 0 0 0
2006 C Argentina 3 3 1 1
2006 C Netherlands 3 3 1 1
2006 C Ivory_Coast 0 0 3 0
2006 C Serbia_and_Montenegro 0 0 0 0
2006 D Portugal 3 3 3 1
2006 D Mexico 3 1 0 1
2006 D Angola 0 1 1 0
2006 D Iran 0 0 1 0
2006 E Iraly 3 1 3 1
2006 E Ghana 0 3 3 1
2006 E Czech_Republic 3 0 0 0
2006 E United_States 0 1 0 0
2006 F Brazil 3 3 3 1
2006 F Australia 3 0 1 1
2006 F Croatia 0 1 1 0
2006 F Japan 0 1 0 0
2006 G Swizerland 1 3 3 1
2006 G France 1 1 3 1
2006 G South_Korea 3 1 0 0
2006 G Togo 0 0 0 0
2006 H Spain 3 3 3 1
2006 H Ukraine 0 3 3 1
2006 H Tunisia 1 0 0 0
2006 H Saudi_Arabia 1 0 0 0
2010 A Uruguay 1 3 3 1
2010 A Mexico 1 3 0 1
2010 A South_Africa 1 0 3 0
2010 A France 1 0 0 0
2010 B Argentina 3 3 3 1
2010 B South_Korea 3 0 1 1
2010 B Greece 0 3 0 0
2010 B Nigeria 0 0 1 0
2010 C United_States 1 1 3 1
2010 C England 1 1 3 1
2010 C Slovenia 3 1 0 0
2010 C Algeria 0 1 0 0
2010 D Germany 3 0 3 1
2010 D Ghana 3 1 0 1
2010 D Australia 0 1 3 0
2010 D Serbia 0 3 0 0
2010 E Netherlands 3 3 3 1
2010 E Japan 3 0 3 1
2010 E Denmark 0 3 0 0
2010 E Cameroon 0 0 0 0
2010 F Paraguay 1 3 1 1
2010 F Slovakia 1 0 3 1
2010 F New_Zealand 1 1 1 0
2010 F Italy 1 1 0 0
2010 G Brazil 3 3 1 1
2010 G Portugal 1 3 1 1
2010 G Ivory_Coast 1 0 3 0
2010 G North_Korea 0 0 0 0
2010 H Spain 0 3 3 1
2010 H Chile 3 3 0 1
2010 H Swizerland 3 0 1 0
2010 H Honduras 0 0 1 0
2014 A Brazil 3 1 3 1
2014 A Mexico 3 1 3 1
2014 A Croatia 0 3 0 0
2014 A Cameroon 0 0 0 0
2014 B Netherlands 3 3 3 1
2014 B Chile 3 3 0 1
2014 B Spain 0 0 3 0
2014 B Australia 0 0 0 0
2014 C Colombia 3 3 3 1
2014 C Greece 0 1 3 1
2014 C Ivory_Coast 3 0 0 0
2014 C Japan 0 1 0 0
2014 D Costa_Rica 3 3 1 1
2014 D Uruguay 0 3 3 1
2014 D Italy 3 0 0 0
2014 D England 0 0 1 0
2014 E France 3 3 1 1
2014 E Swizerland 3 0 3 1
2014 E Ecuador 0 3 1 0
2014 E Honduras 0 0 0 0
2014 F Argentina 3 3 3 1
2014 F Nigeria 1 3 0 1
2014 F Bosnia_and_Herzegovina 0 0 3 0
2014 F Iran 1 0 0 0
2014 G Germany 3 1 3 1
2014 G United_States 3 1 0 1
2014 G Portugal 0 1 3 0
2014 G Ghana 0 1 0 0
2014 H Belgium 3 3 3 1
2014 H Algeria 0 3 1 1
2014 H Russia 1 0 1 0
2014 H South_Korea 1 0 0 0
"

2018-06-12

MikuHatsune2018-06-12

なぜ人は新刊を落とすのか

こんなツイートを見かけた。

C93 で買いそびれていた。

 

808 のサークルにアンケート調査をして、新刊を落としたサークルと落とさなかったサークルで、制作にとりかかった時期の分布を出している。

ここで、落とさなかったサークルは平均64日前から、落としたサークルは平均48日前から制作に取り組んでいる、と考察しているが、ヒストグラムの分布を見るにニ峰性で、要約統計量として平均を出すのも、中央値や最頻値を出すのもよろしくない。

落としたサークルの平均48日が、ヒストグラムのみぎから3番目の40-59日のビンに入るので、まあ、悪くないとして、落としたことのないサークルの平均64日が、最も頻度の少ない60-79日のビンに収まってしまうことに違和感がバリバリである。

f:id:MikuHatsune:20180612162623p:image

 

ここで、数値から離れて、コミケ同人誌を制作するというモチベーションを考えてみると、60-70日前というのは、夏冬ともにコミケ当選の時期と被っており、当選前から制作しようと考えるか、当選後に制作しようと考えるか、という仮定をおいてみる。

すると、このヒストグラムは、新刊を落としたことのある/ない、に加えて、コミケ当選(60-70日程度前)の前/後で制作を開始する、が混合した分布と考えられる。ビンの高さは左右非対称のように見えるので、正規分布ではなくガンマ分布を想定した。

 

ガンマ分布は形パラメータをふたつとり、G(s_1,s_2) と書くと、コミケ当選後に制作を始めるサークルの割合をpとすると、コミケ当選前から制作を始めるサークルの割合は1-p となり、混合分布はpG_1+(1-p)G_2 となる。これをヒストグラムの累積分布から、落としたことのあるサークル、落としたことのないサークル、のパラメータ各2つずつ、とp ゴリ押しで推定する。

 

推定した結果、pコミケ当選後から制作するサークルのガンマ分布パラメータふたつ、コミケ当選前から制作するサークルのガンマ分布パラメータふたつ、の順で並べると以下の通りである。

[[1]]
[1]  0.8531444  4.0315755  0.1158495 95.3276615  0.9989746

[[2]]
[1]  0.7308913  7.7311609  0.1909049 98.0969186  1.0017727

混合ガンマ分布での平均値は以下の通りで、コミケ当選前後で制作を開始するサークルを別にしなければ、12日程度の差がある。

[1] 43.69920 55.95106

さて、コミケ当選前後で制作を開始するという分類を入れたのでこの別で見てみる。

f:id:MikuHatsune:20180612163620p:image

以下の行列は制作開始前の平均日数である。行はコミケ当選後に制作を始めるかコミケ当選前から制作を始めるサークルかで、列は新刊を落としたことがあるか、ないか、である。コミケ当選後では6日程度の差があるが、コミケ当選前では2日である。総製作期間に対する割合を考慮するとあまり差がないと考えられる。

         [,1]     [,2]
[1,] 34.80011 40.49744
[2,] 95.42551 97.92333

声優統計を書いていたときはコミケ当選前からネタ探しはして、実際にとりかかるのはコミケ当選が発表されてからすぐだったのでだいたい60日前くらいだろうか。8回書かせてもらって落としたことはない(ドヤァア

txt <- "
19 71 23
39 169 105
59 117 119
79 8 11
99 57 63
100 22 43
"

d <- as.matrix(read.table(text=txt, row=1))

day <- c(0, as.numeric(rownames(d)))

cols <- c("skyblue", "hotpink")
led <- sprintf("新刊を落とした経験が%s", c("ある", "ない"))
yl <- c(0, 180)
ab <- seq(0, 180, by=20)
barplot(t(d), beside=TRUE, col=NA, border=NA, las=1, ylim=yl, xaxt="n", yaxt="n", xlab="日数")
abline(h=ab, col=grey(0.8))
barplot(t(d), beside=TRUE, col=cols, las=1, ylim=yl, add=TRUE, yaxt="n")
axis(2, at=ab, labels=ab, las=1)
legend("topright", legend=led, col=cols, pch=15, bg="white")

# cumulative density
dens <- sweep(apply(d, 2, cumsum), 2, colSums(d), "/")

# 最適化でゴリ押しする
y <- head(dens[,1], -1)
resid <- function(par){
  yhat <- par[1]*pgamma(day[2:(length(day)-1)], par[2], par[3]) + (1-par[1])*pgamma(day[2:(length(day)-1)], par[4], par[5])
  sum((y-yhat)^2)
}

niter <- 10000
res <- vector("list", niter)
fun <- function() optim(c(runif(1), rpois(1, 40), runif(1, 0, 10), rpois(1, 90), runif(1, 0, 10)), resid, method="BFGS")
for(i in seq(niter)){
  res[[i]] <- try(fun(), silent=TRUE)
}

res <- res[sapply(res, length) > 1] # 推定がうまくいかなかった場合を除く
res[[which.min(sapply(res, "[[", "value"))]]
ps <- res[[which.min(sapply(res, "[[", "value"))]]$par
# ps には新刊を落としたことがある場合と落としたことがない場合をlist でいれる

X <- 0:150
es <- mapply(function(z) z[1]*dgamma(X, z[2], z[3]) + (1-z[1])*dgamma(X, z[4], z[5]), ps)
matplot(es, type="l", col=cols, lty=1, lwd=4, xaxt="n", xlab="日数", ylab="Density")
abline(v=seq(0, 100, by=20), lty=3, col=1)
legend("topright", legend=led, col=cols, pch=15, bg="white")
axis(1, at=seq(0, max(X), by=10), seq(0, max(X), by=10))

# コミケ当選前後をわけない
colSums(es * X)
# コミケ当選前後をわける
mapply(function(z) c(z[2]/z[3], z[4]/z[5]), ps)

2018-05-30

ホーム&アウェー方式のアウェーゴールは妥当なのか

チャンピオンズリーグを見ていた。

CL というか本拠地があってホームでもアウェーでも試合を行うタイプのスポーツにあるのが、先に言ったホームとアウェーで戦い方というか試合に影響する部分が大きい、ということである。

サッカーは特にホーム&アウェーの影響が大きい(アウェーゴール)らしいのだが、ホームとアウェーで各々入れ替えて試合をするのでそれはいいのでは…とは思うが、それでも、先にホームで戦うのかそれともアウェーなのかも影響する、と言われているらしい。

というわけでデータを集計して解析した。

データはwikipediaのCL のページをクローリングするが、1992年から試合形式が変わっているようなので、1992年から2017年まで、ホーム&アウェー形式で試合をしている部分をパースしてくる。

ちなみに決勝戦は事前に決められた中立地で一発勝負で行われるため、決勝戦は除く。

# R
i <- 1992:2017
j <- c(93:99,0:18)
sprintf("https://en.wikipedia.org/wiki/%d%s%02d_UEFA_Champions_League", i, "&#8211;", j)
# ターミナル
wget -i url.txt
for i in $(ls [0-9]*)
do
  name=$(echo $i | grep -o '^[0-9]*')
  mv $i $name
done

その後、たぶんエンダッシュ – を目印にひたすら試合結果をパースする(後述)。

 

ホーム&アウェーで各2回ずつあるので分けて集計しておく。

何も考えずにホームとアウェーで平均ゴール数を出すと

[1] 1.623472 1.050530

となる。単純に1.5倍くらいホームのほうが点を取りやすい。

f:id:MikuHatsune:20180529204335p:image

 

1試合あたりのゴール数はポアソン分布に従うと仮定すると、平均ゴール数はホームのほうがよく取れる。

poisson.test(colSums(goal))
	Comparison of Poisson rates

data:  colSums(goal) time base: 1
count1 = 3984, expected count1 = 3281, p-value < 2.2e-16
alternative hypothesis: true rate ratio is not equal to 1
95 percent confidence interval:
 1.470317 1.624500
sample estimates:
rate ratio 
  1.545384 

 

最初にホーム/アウェーにくるか、2試合目に来るか(1st leg と2nd leg という)で得点数に差があるかを検討する。単純にホーム&アウェーの順番も格納しておいたので個別にやるだけ。

   H_1st    A_1st    A_2nd    H_2nd 
1.443358 1.064385 1.036675 1.803586

得点数の分布は、2nd leg でホームにくる場合が最も多く点をとっているようである。平均ゴール数も1.8点と、2試合合計で本当に勝負がかかっているときがもっとも点をとる、ようである。

f:id:MikuHatsune:20180529204332p:image

個別にポアソン検定をすると、

H_1st vs A_1st では、別のチームの比較だが、1st leg でホームのほうが点を取る。

H_1st vs A_2nd では、同じチームの比較だが、1st leg でホームのときのほうが点を取る。

H_1st vs H_2nd では、別のチームの比較だが、2nd leg でホームのほうが点を取る。

A_1st vs A_2nd では、別のチームの比較だが、アウェーでの点を取る数は変わらない。

A_1st vs H_2nd では、同じチームの比較だが、2nd leg でホームのときのほうが点を取る。

A_2nd vs H_2nd では、別のチームの比較だが、2nd leg でホームのほうが点を取る。

というわけで、2nd leg でホームのときが一番点を取れそうな気がするようである。

 

時系列でデータをとったので、Elo rating によるチームの強さ定量化をしたかったがやれていない。

[1] "H_1st vs A_1st"

	Comparison of Poisson rates

data:  colSums(score)[cmb[, i]] time base: 1
count1.H_1st = 1771, expected count1 = 1538.5, p-value < 2.2e-16
alternative hypothesis: true rate ratio is not equal to 1
95 percent confidence interval:
 1.261775 1.457655
sample estimates:
rate ratio.H_1st 
        1.356049 

[1] "H_1st vs A_2nd"

	Comparison of Poisson rates

data:  colSums(score)[cmb[, i]] time base: 1
count1.H_1st = 1771, expected count1 = 1521.5, p-value < 2.2e-16
alternative hypothesis: true rate ratio is not equal to 1
95 percent confidence interval:
 1.294792 1.497468
sample estimates:
rate ratio.H_1st 
        1.392296 

[1] "H_1st vs H_2nd"

	Comparison of Poisson rates

data:  colSums(score)[cmb[, i]] time base: 1
count1.H_1st = 1771, expected count1 = 1992, p-value = 2.675e-12
alternative hypothesis: true rate ratio is not equal to 1
95 percent confidence interval:
 0.7513686 0.8522624
sample estimates:
rate ratio.H_1st 
       0.8002711 

[1] "A_1st vs A_2nd"

	Comparison of Poisson rates

data:  colSums(score)[cmb[, i]] time base: 1
count1.A_1st = 1306, expected count1 = 1289, p-value = 0.5157
alternative hypothesis: true rate ratio is not equal to 1
95 percent confidence interval:
 0.9497088 1.1100190
sample estimates:
rate ratio.A_1st 
         1.02673 

[1] "A_1st vs H_2nd"

	Comparison of Poisson rates

data:  colSums(score)[cmb[, i]] time base: 1
count1.A_1st = 1306, expected count1 = 1759.5, p-value < 2.2e-16
alternative hypothesis: true rate ratio is not equal to 1
95 percent confidence interval:
 0.5507167 0.6322103
sample estimates:
rate ratio.A_1st 
       0.5901491 

[1] "A_2nd vs H_2nd"

	Comparison of Poisson rates

data:  colSums(score)[cmb[, i]] time base: 1
count1.A_2nd = 1272, expected count1 = 1742.5, p-value < 2.2e-16
alternative hypothesis: true rate ratio is not equal to 1
95 percent confidence interval:
 0.5360617 0.6161042
sample estimates:
rate ratio.A_2nd 
       0.5747854 

# python3 data preprocessing
import re, glob, os
score = re.compile("\>\d+&#8211;\d+")
game = ["qualifying", "R16", "QuaterFinal", "SemiFinal"]
w0 = open("out.txt", "w")

for y in range(1992, 2018):
  f = open(str(y), "rU")
  res = []
  g0 = []
  scoreflag = -1
  leg = False
  table = False
  for g in f:
    g0 += [g]
    tmp = score.findall(g)
    if "1st leg" in g:
      leg = True
      table = True
      scoreflag = -1
    if len(tmp) > 0 and leg and table:
      scoreflag += 1
      if len(tmp) > 1:
        tmp = tmp[-1]
      tmp = [re.sub("&#8211;", "-", tmp[0])]
      tmp = [re.sub(">", "", tmp[0])]
      if scoreflag % 3 < 2:
        res += [re.sub("<[^>]*>", "", g0[-2]).strip(), tmp[0]]
      elif scoreflag % 3 == 2:
        res += [tmp[0]]
    if '</table' in g:
      leg = False
      table = False
  N = [int(len(res)/5)-14, 8, 4, 2]
  lab = sum(map(lambda x, y: [x]*y, game, N), [])
  for i in range(int(len(res)/5)):
    line = [str(y), lab[i], res[j*i], res[j*i+1], res[j*i+2], res[j*i+3], res[j*i+4]]
    print(line)
    w0.write("\t".join(line) + "\n")

w0.close()

# R 解析
# R analysis
# ひだりの列が初戦のホームチーム
# スコアはチームの並びを保ったまま記載される
# 1列目:初戦がホーム
# 2列目:初戦がアウェイ
# 3列目:2戦目がアウェイ
# 4列目:2戦目がホーム

library(igraph)
library(rstan)
dat <- read.delim("out.txt", stringsAsFactors=FALSE, header=FALSE)


gr <- graph_from_edgelist(as.matrix(dat[,c(3,5)]), directed=FALSE)
adj <- as.matrix(as_adjacency_matrix(gr))


sb <- subset(dat, V2=="R16")
sapply(split(sb[c(3,5)], sb$V1), unlist)


score <- t(sapply(sapply(apply(dat[,6:7], 1, paste, collapse="-"), strsplit, "-"), as.numeric, USE.NAMES=FALSE))
colnames(score) <- c("H_1st", "A_1st", "A_2nd", "H_2nd")


goal <- cbind(c(score[,c(1,4)]), c(score[,2:3]))
Max <- max(goal)
tab <- mapply(function(z) table(factor(goal[,z], 0:Max)), seq(ncol(goal)))
me <- colMeans(goal)
d <- mapply(dpois, list(0:Max), me)


# ゴール数ヒストグラム
cols <- c("blue", "green")
b0 <- barplot(t(tab), beside=TRUE, col=cols, xlab="Goal")
legend("topright", legend=sprintf("%s: %.3f", c("Home", "Away"), me), pch=15, col=cols, title="Expected goal", cex=1.5)
legend("right", legend=sprintf("%d games, %d goals", nrow(goal), sum(goal)), bty="n", cex=1.5)
title(sprintf("Home and Away goals of Champions League %d - %d", min(dat$V1), max(dat$V1)))
for(i in seq(ncol(tab))){
  lines(colMeans(b0), d[,i]*sum(tab[,i]), col=cols[i], lwd=3)
}

poisson.test(colSums(goal))


# 1st 2nd leg のHA すべて区別する
cmb <- combn(4, 2)
colMeans(score)
mat <- diag(0, ncol(score))
rownames(mat) <- colnames(mat) <- colnames(score)
for(i in 1:ncol(cmb)){
  res <- poisson.test(colSums(score)[ cmb[,i] ])
  mat[ cmb[1,i], cmb[2,i] ] <- res$p.value
  print(sprintf("%s vs %s", colnames(score)[cmb[1,i]], colnames(score)[cmb[2,i]]))
  print(res)
}

tab <- mapply(function(z) table(factor(score[,z], 0:Max)), colnames(score))
me <- colMeans(score)

# ゴール数ヒストグラム
cols <- c("blue", "green", "pink", "orange")
b0 <- barplot(t(tab), beside=TRUE, col=cols, xlab="Goal")
legend("topright", legend=sprintf("%s: %.3f", colnames(tab), me), pch=15, col=cols, title="Expected goal", cex=1.5)
title(sprintf("Home and Away goals of Champions League %d - %d", min(dat$V1), max(dat$V1)))
for(i in seq(ncol(tab))){
  lines(colMeans(b0), d[,i]*sum(tab[,i]), col=cols[i], lwd=3)
}

2018-05-29

MikuHatsune2018-05-29

ギブスサンプリング

Metropolis-Hastings サンプリングをやったので、ギブスサンプリングをやってみる。

ある変数X=¥{x_1,x_2,¥dots,x_{i-1},x_i,x_{i+1},¥dots,x_m¥} について、i 番目を取り除いたX_{-i}

X_i|X_{-i} ¥sim x_1,x_2,¥dots,x_{i-1},~~~x_{i+1},¥dots,x_mi 番目が抜けている)

で順次サンプリングして、その値を入れなおしてまたサンプリングする、を1¥dots m について行う。

 

結局、あるひとつi 以外のすべてを固定して、ひとつサンプリングすることになる。多次元正規分布の場合は、こちらを参考に

¥mu_{i|-i}¥leftarrow ¥mu_{i}+¥Sigma_{-i-i}^{-1}(X_{-i}-¥mu_{-i})

¥Sigma_{i|-i}¥leftarrow ¥Sigma_{ii}-¥Sigma_{-i-i}^{-1}¥Sigma_{-ii}

とする。ただし¥underbrace{¥Sigma_{-i-i}^{-1}}_{1¥times (m-1)}=¥underbrace{¥Sigma-ii}_{1¥times(m-1)}¥underbrace{¥Sigma_{-i-i}}_{(m-1)¥times(m-1)} である。なので¥mu_{i|-i}¥Sigma_{i|-i} はそれぞれ1¥times 1 の大きさというかスカラーで、各X_i はひとつずつサンプリングされる。

スクリプトはここからパクる。

 

初期値が真の分布からは遠いところにしておいて、収束するまでにちょっと時間がかかるようにする。

Metropolis-Hastings(MH)では、真の分布に近づくまでに少し時間がかかる。サンプリングされた分布は、xy 同時にサンプリングしているので、移動方向は斜めである。

Gibbs (GB)は、各X_i でサンプリングしているので、動き方は各x, y 方向にジグザグである。詳細釣り合い条件や移動確率が1であることはここでは取り上げないが、動きやすいのでMH より収束に向かいやすい。

この2次元の例では、一発目から真の分布に取り込まれている様子がわかる。

f:id:MikuHatsune:20180529171307p:image

 

多次元に拡張するが、図示することを考えて3次元にしてみる。MH もGB もひだり下のマゼンタのところからサンプリングを開始したが、HM は最初の50点の時点ではまだ収束していないが、GB は7点目くらいからもう収束している。

f:id:MikuHatsune:20180529171355p:image

MH の移動方向は斜めだが、GB では各点がサンプリングされるまでに、xyz でそれぞれX_i をサンプリングしているので、カクカクになっている。

# 2D
# 初期値
x0 <- c(-5, 2)
iter <- 1500 # 繰り返し回数
x <- matrix(0, iter, 2) # ランダムウォークの座標
dratio <- rep(0, iter) # 確率密度比
out <- rep(0, iter) # 確率密度比にしたがって確率的に動いたか、動かなかったかの記録
x[1,] <- x0

# 2次元正規分布の相関
rho <- 0.7
sig <- matrix(c(1, rho, rho, 1), 2)
mu <- c(4, 5) # 2次元正規分布の各パラメータの平均

# Metrololis-Hastings
for(i in 2:iter){
  walk <- runif(2, -0.5, 0.5) # 歩く量
  v <- x[i-1,] + walk # 次の候補点
  dratio[i-1] <- dmvnorm(v, mean=mu, sigma=sig)/dmvnorm(x[i-1,], mean=mu, sigma=sig) # 確率密度比
  out[i-1] <- rbinom(1, 1, min(1, dratio[i-1])) # 動いたか、動かなかったか
  x[i,] <- x[i-1,] + out[i-1]*walk
}

kd1 <- kde2d(x[,1], x[,2], c(bandwidth.nrd(x[,1]), bandwidth.nrd(x[,2])), n=1000)
cols <- jet.colors(iter)
plot(x, type="p", pch=16, cex=0.6, col=2, xlab="x1", ylab="x2", main="Metropolis-Hastings sampling", xlim=xl, ylim=yl)
for(i in 2:iter){
  segments(x[i-1,1], x[i-1,2], x[i,1], x[i,2], col=cols[i])
}
points(x[1,1], x[1,2], pch="★", col=6)
points(mu[1], mu[2], pch="★", col=5)
contour(kd1, add=TRUE, col=1)

# Gibbs
X <- matrix(x0, nr=1)
for(j in 2:iter){
  x <- X[j-1,]
  for(i in seq(mu)){
    s <- sig[-i, i] %*% solve(sig[-i, -i])   # Σ_ab Σ_bb ^ -1
    # (PRML 2.81) μ_a|b = μ_a + Σ_ab Σ_bb ^ -1 (x_b - μ_b)
    mu_a_b <- mu[i] + s %*% (x[-i] - mu[-i])
    # (PRML 2.82) Σ_a|b = Σ_aa - Σ_ab Σ_bb ^ -1 Σ_ba
    sig_a_b <- sig[i, i] - s %*% sig[i, -i]
   # [Gibbs] x_a 〜 p(x_a|x_{-a}) = N(μ_a|b, Σ_a|b)
    x[i] <- rnorm(1, mu_a_b, sqrt(sig_a_b))
  }
  X <- rbind(X, x)
}
colnames(X) <- sprintf("V%d", seq(ncol(X)))

kd2 <- kde2d(X[,1], X[,2], c(bandwidth.nrd(X[,1]), bandwidth.nrd(X[,2])), n=1000)
cols <- jet.colors(iter)
plot(X, type="p", pch=16, cex=0.6, col=2, xlab="x1", ylab="x2", main="Gibbs sampling", xlim=xl, ylim=yl)
#lines(x, col=cols)
for(i in 2:iter){
  for(j in seq(ncol(X))){
    y <- rbind(replace(X[i-1,], 0:(j-1), X[i, 0:(j-1)]), replace(X[i-1,], 1:j, X[i, 1:j]))
    segments(y[1,1], y[1,2], y[2,1], y[2,2], col=cols[i])
  }
}
points(X[1,1], X[1,2], pch="★", col=6)
points(mu[1], mu[2], pch="★", col=5)
contour(kd2, add=TRUE, col=1)
# 3D
# 初期値
x0 <- c(-5, 2, -4)

iter <- 1500 # 繰り返し回数
x <- matrix(0, iter, length(x0)) # ランダムウォークの座標
dratio <- rep(0, iter) # 確率密度比
out <- rep(0, iter) # 確率密度比にしたがって確率的に動いたか、動かなかったかの記録
x[1,] <- x0
# 3次元の多次元正規分布を適当に作る
rho <- c(0.7, 0.6, 0.9)
sig <- diag(0, length(x0))
sig[lower.tri(sig)] <- rho
sig <- sig + t(sig)
diag(sig) <- rep(1, length(x0))
mu <- c(4, 5, 3)

# Metropolis-Hastings
for(i in 2:iter){
  walk <- runif(length(x0), -0.5, 0.5) # 歩く量
  v <- x[i-1,] + walk # 次の候補点
  dratio[i-1] <- dmvnorm(v, mean=mu, sigma=sig)/dmvnorm(x[i-1,], mean=mu, sigma=sig) # 確率密度比
  out[i-1] <- rbinom(1, 1, min(1, dratio[i-1])) # 動いたか、動かなかったか
  x[i,] <- x[i-1,] + out[i-1]*walk
}
colnames(x) <- sprintf("V%d", seq(ncol(x)))

cols <- jet.colors(iter)
plot3d(x, type="p", pch=16, size=0.01)
points3d(head(x, 50), size=5)
title3d("Metropolis-Hastings sampling", line=3)
for(i in 2:iter){
  segments3d(x[(i-1):i,], col=cols[i])
}
spheres3d(x[1,], col=6, radius=0.5)

# Gibbs
X <- matrix(x0, nr=1)
for(j in 2:iter){
  x <- X[j-1,]
  for(i in seq(mu)){
    s <- sig[-i, i] %*% solve(sig[-i, -i])   # Σ_ab Σ_bb ^ -1
    # (PRML 2.81) μ_a|b = μ_a + Σ_ab Σ_bb ^ -1 (x_b - μ_b)
    mu_a_b <- mu[i] + s %*% (x[-i] - mu[-i])
    # (PRML 2.82) Σ_a|b = Σ_aa - Σ_ab Σ_bb ^ -1 Σ_ba
    sig_a_b <- sig[i, i] - s %*% sig[i, -i]
   # [Gibbs] x_a 〜 p(x_a|x_{-a}) = N(μ_a|b, Σ_a|b)
    x[i] <- rnorm(1, mu_a_b, sqrt(sig_a_b))
  }
  X <- rbind(X, x)
}
colnames(X) <- sprintf("V%d", seq(ncol(X)))
# ギザギザの描画
Y <- X[1, , drop=FALSE]
for(i in 2:nrow(X)){
  for(j in 1:ncol(X)){
    y <- rbind(replace(X[i-1,], 0:(j-1), X[i, 0:(j-1)]), replace(X[i-1,], 1:j, X[i, 1:j]))
    Y <- rbind(Y, y)
  }
}

cols <- jet.colors(iter)
plot3d(X, size=0.01)
points3d(head(X, 50), size=5)
title3d("Gibbs sampling", line=3)
lines3d(Y, col=rep(cols, each=3))
spheres3d(x[1,], col=6, radius=0.5)