読んだ。
A comparison of single-cell trajectory inference methods: towards more accurate and robust tools.
single cell のRNAseq などからデータを取得して、細胞分化系統樹を解析するのに多種多様な手法やパッケージが出ている。
59手法を試してみて、分化系統樹にどのような仮定をおくかをフローチャート形式にして、どういう場合にどういう手法を使うのがよいか、まで提言している。
とりあえずSlingshot, monocle, TSCAN, cellTree あたりを使っておけばよさそうである。
vim 8.1 にしたら \ の後は / か ? か & でなければなりませんって文句言われるし"Vim: Caught deadly signal ABRT" って言われて突然落ちる
vim 8.0 でvimrc をいい感じにしていたのに、ppa (https://launchpad.net/~jonathonf/+archive/ubuntu/vim) を
sudo add-apt-repository ppa:jonathonf/vim
でvim-gnome をインストールすると 8.1 がインストールされるようになっている。
そうすると、いままで 8.0 ではエラーが出ていなかったのに、
\ の後は / か ? か & でなければなりません
って起動時に毎回文句言われるようになった。これは
vim -N
で回避できるらしいが、肝心のスクリプトを書いているときに
Vim: Caught deadly signal ABRT
と言われて突然落ちる。python を書いているときに頻発する。
パッチの問題とか言われているようだが、パッチの当て方がよくわからないので vim 8.0 がapt-get できるようにしたいが、ppa をいくらいじっても 8.1 しか出てこない。
こちらのppa (https://launchpad.net/~laurent-boulard/+archive/ubuntu/vim) を使うと
sudo add-apt-repository ppa:laurent-boulard/vim
sudo apt-cache policy vim-gnome
vim-gnome: インストールされているバージョン: 2:8.0.1520-1~xenial~lboulard+1 候補: 2:8.0.1520-1~xenial~lboulard+1 バージョンテーブル: *** 2:8.0.1520-1~xenial~lboulard+1 500 500 http://ppa.launchpad.net/laurent-boulard/vim/ubuntu xenial/main amd64 Packages 500 http://ppa.launchpad.net/laurent-boulard/vim/ubuntu xenial/main i386 Packages 100 /var/lib/dpkg/status 2:7.4.1689-3ubuntu1.2 500 500 http://jp.archive.ubuntu.com/ubuntu xenial-updates/main amd64 Packages 500 http://security.ubuntu.com/ubuntu xenial-security/main amd64 Packages 2:7.4.1689-3ubuntu1 500 500 http://jp.archive.ubuntu.com/ubuntu xenial/main amd64 Packages
ようやく出るようになった。
いやでもやっぱりpython で落ちる。
決勝トーナメントに向けて初戦が大事というが初戦はどれくらい大事なのか
2018FIFAワールドカップが始まった。日本の決勝トーナメント進出は初戦に勝利できるかにかかっている! とかなんとかよく聞くが、実際に初戦に勝利するのはどれだけ大事かは定量的に言われていない気がする。
2002年と2010年に決勝トーナメントに進出したので8年周期のジンクスがある! みたいに言っていたTV があったが何を言っているのかさっぱりわからなかった。
というわけで、32チームを4チーム8グループに分けて、上位2チームが決勝トーナメントに進出するような仕様になった1998年から2014年まで5大会分の延べ160チームの勝敗表から、初戦に勝利すると決勝トーナメント進出にどれだけ有利になるかを調べる。
初戦に勝つのを、決勝トーナメントに進出するのを とすると、初戦に勝って決勝トーナメントに進出するのは である。
地味に数えればよい。行は初戦での勝点、つまり上から負け、引き分け、勝ち、であり、列は決勝トーナメントに進出したからの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 をしたけどあまり意味がない。
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 "
なぜ人は新刊を落とすのか
こんなツイートを見かけた。
808 のサークルにアンケート調査をして、新刊を落としたサークルと落とさなかったサークルで、制作にとりかかった時期の分布を出している。
ここで、落とさなかったサークルは平均64日前から、落としたサークルは平均48日前から制作に取り組んでいる、と考察しているが、ヒストグラムの分布を見るにニ峰性で、要約統計量として平均を出すのも、中央値や最頻値を出すのもよろしくない。
落としたサークルの平均48日が、ヒストグラムのみぎから3番目の40-59日のビンに入るので、まあ、悪くないとして、落としたことのないサークルの平均64日が、最も頻度の少ない60-79日のビンに収まってしまうことに違和感がバリバリである。
ここで、数値から離れて、コミケで同人誌を制作するというモチベーションを考えてみると、60-70日前というのは、夏冬ともにコミケ当選の時期と被っており、当選前から制作しようと考えるか、当選後に制作しようと考えるか、という仮定をおいてみる。
すると、このヒストグラムは、新刊を落としたことのある/ない、に加えて、コミケ当選(60-70日程度前)の前/後で制作を開始する、が混合した分布と考えられる。ビンの高さは左右非対称のように見えるので、正規分布ではなくガンマ分布を想定した。
ガンマ分布は形パラメータをふたつとり、 と書くと、コミケ当選後に制作を始めるサークルの割合をとすると、コミケ当選前から制作を始めるサークルの割合は となり、混合分布は となる。これをヒストグラムの累積分布から、落としたことのあるサークル、落としたことのないサークル、のパラメータ各2つずつ、と ゴリ押しで推定する。
推定した結果、、コミケ当選後から制作するサークルのガンマ分布パラメータふたつ、コミケ当選前から制作するサークルのガンマ分布パラメータふたつ、の順で並べると以下の通りである。
[[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)
就学年数が長いから近視になるのか近視になるから就学年数が長いのか
読んだ。
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.
行列演算だから速いだろうと油断してはいけない
回帰分析の最小二乗法による係数の推定は
で求められる。ここで、標本数、パラメータ数 とすると
である。
普通に演算すると、 に をかけて をかける、という順番だが、 として先に と をかけて としてから にかけるのが高速である。
20180607追記
スーパーポスドクが計算量について教えてくれた。
行列と行列の掛け算では、個の成分を計算することになり、各成分について回の掛け算と回の足し算が必要になる。
行列、行列、行列の掛け算を行う場合、
と計算すると、は行列だから、
掛け算は 回
足し算は 回
と計算すると、は行列だから、
掛け算は 回
足し算は 回
前者だと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
先に を行う場合は
system.time(replicate(iter, solve(t(X) %*% X) %*% (t(X) %*% y)))
ユーザ システム 経過 41.333 37.563 6.672
ホーム&アウェー方式のアウェーゴールは妥当なのか
チャンピオンズリーグを見ていた。
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, "–", 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+–\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("–", "-", 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) }