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 の行を見ればよい。初戦で負けると 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 をしたけどあまり意味がない。

　

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

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

R アニメ統計数学

MikuHatsune2018-06-12

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

C93 で買いそびれていた。
　
808 のサークルにアンケート調査をして、新刊を落としたサークルと落とさなかったサークルで、制作にとりかかった時期の分布を出している。
ここで、落とさなかったサークルは平均64日前から、落としたサークルは平均48日前から制作に取り組んでいる、と考察しているが、ヒストグラムの分布を見るにニ峰性で、要約統計量として平均を出すのも、中央値や最頻値を出すのもよろしくない。
落としたサークルの平均48日が、ヒストグラムのみぎから3番目の40-59日のビンに入るので、まあ、悪くないとして、落としたことのないサークルの平均64日が、最も頻度の少ない60-79日のビンに収まってしまうことに違和感がバリバリである。

　
ここで、数値から離れて、コミケで同人誌を制作するというモチベーションを考えてみると、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

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

以下の行列は制作開始前の平均日数である。行はコミケ当選後に制作を始めるかコミケ当選前から制作を始めるサークルかで、列は新刊を落としたことがあるか、ないか、である。コミケ当選後では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-06-11

就学年数が長いから近視になるのか近視になるから就学年数が長いのか

医学統計遺伝医学

読んだ。
BMJ. 2018 Jun 6;361:k2022.
観察研究では、どちらも関連があるが、この論文の結論から言えば就学年数が長いから近視になり、近視だから就学年数が長い、ということにはならない。これを近視に関する多型でメンデリアンランダマイゼーションした場合（Nat Genet. 2016 Jul;48(7):709-17.）と、就学年数に関する多型でメンデリアンランダマイゼーションした場合（Nature. 2016 May 26;533(7604):539-42.）で解析し、因果関係を推定している。
　
メンデリアンランダマイゼーションを双方向に行って因果関係を推定するのは、
例えばCRP とBMI はどちらも高いことが多いが、これはCRP の多型でランダマイズしてもBMI は増加しないが、BMI に関係するFTO という多型でランダマイズするとCRP が増加する、ということでやっている。
Int J Obes (Lond). 2011 Feb;35(2):300-8.

2018-05-31

行列演算だから速いだろうと油断してはいけない

R Rを使いこなす医学

回帰分析の最小二乗法による係数の推定は
$b=(X^TX)^{-1}X^Ty$
で求められる。ここで、標本数 $n$ 、パラメータ数 $p$ とすると
$\underbrace{b}_{p\times 1}=(\underbrace{X^TX}_{p\times p})^{-1}\underbrace{X^T}_{p \times n}\underbrace{y}_{n\times 1}$
である。
普通に演算すると、 $(X^TX)^{-1}$ に $X^T$ をかけて $y$ をかける、という順番だが、 $(X^TX)^{-1}(X^Ty)$ として先に $X^T$ と $y$ をかけて $(X^Ty)$ としてから $(X^TX)^{-1}$ にかけるのが高速である。
　
20180607追記
スーパーポスドクが計算量について教えてくれた。

$l\times m$ 行列と $m\times n$ 行列の掛け算では、 $l\times n$ 個の成分を計算することになり、各成分について $m$ 回の掛け算と $m-1$ 回の足し算が必要になる。

$p\times p$ 行列 $A$ 、 $p\times n$ 行列 $B$ 、 $n\times 1$ 行列 $C$ の掛け算を行う場合、

$(AB)C$ と計算すると、 $AB$ は $p\times n$ 行列だから、
掛け算は $p\times p \times n + n\times p\times 1 = p(np+n)$ 回
足し算は $(p-1)\times p\times n + (n-1)\times p\times 1 = p(np-1)$ 回
　
$A(BC)$ と計算すると、 $BC$ は $p\times 1$ 行列だから、
掛け算は $n\times p\times 1 + p\times p\times 1 = p(n+p)$ 回
足し算は $(n-1)\times p\times 1 + (p-1)\times p\times 1 = p(n+p-2)$ 回

前者だと3次の項が出てくるが、後者だと2次まで。

n <- 1000
p <- 100
X <- matrix(runif(n*p), n, p)
y <- matrix(runif(n), n, 1)
iter <- 3000

前からかけていく場合は

system.time(replicate(iter, solve(t(X) %*% X) %*% t(X) %*% y))

   ユーザ   システム       経過  
    64.515     54.919     10.156

先に $(X^Ty)$ を行う場合は

system.time(replicate(iter, solve(t(X) %*% X) %*% (t(X) %*% y)))

   ユーザ   システム       経過  
    41.333     37.563      6.672

2018-05-30

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

R 統計数学

チャンピオンズリーグを見ていた。
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倍くらいホームのほうが点を取りやすい。

　
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試合合計で本当に勝負がかかっているときがもっとも点をとる、ようである。

個別にポアソン検定をすると、
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

ギブスサンプリング

R 機械学習 rstan 統計 Rを使いこなす数学

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_m$ （ $i$ 番目が抜けている）
で順次サンプリングして、その値を入れなおしてまたサンプリングする、を $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次元の例では、一発目から真の分布に取り込まれている様子がわかる。

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

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)

2018-05-25

ものすごいわかりやすかったKullback-Leibler divergence

R 機械学習統計数学

情報量を $-logQ(x)$ と置くと、
単調に減少する
ふたつの分布が独立ならば、加算的に扱える
というので、 $-logQ(x)$ と書いておく。
　
符号化してその平均長を考えると
$\sum-Q(x)logQ(x)$ はエントロピーである。
さてここで、 $Q(x)$ はわからないけど、理論的にこれならば、上のエントロピーが成り立つ。しかし、何度も言うが $Q(x)$ はわからないので、ひとまず $P(x)$ を考える。これの平均長は
$-\sum xlogP(x)$
だが、 $x\sim Q(x)$ （実際にデータとして得ている $x$ は、実態はわからないけれども $Q(x)$ からサンプリングされている）と書けるので、結局 $-\sum Q(x)logP(x)$ である。
このふたつの、（真の分布 $Q(x)$ のときの平均長）-（真の分布がわからないけれどもとりあえず $P(x)$ とおいた平均長）を計算すると
$D(Q, P)=-\sum Q(x)log\frac{P(x)}{Q(x)}$
となる。
これは常に0 以上で、 ${P(x)}{Q(x)}=X$ とすれば、 $y=X-1$ と $y=logX$ のふたつの関数を考えれば成り立ち、0 になるときは恒等的に $P(x)=Q(x)$ のときである。