奨励会時の昇級が分かったので追記した。
藤井棋士が最後の王座戦に勝利し、全タイトル戦を制覇し八冠となった。
ここで、八冠になるまでに藤井〇〇という称号がいくつかあるわけだが、それぞれを個別の藤井棋士とみなすと、称号保持の期間がどれくらいあるか日数が出る。
というわけで生存曲線を作ってみる。
日本棋院とwikiを参考に、藤井棋士の段位や称号を得た時期をまとめる。アマチュア時代の級はよくわからなかったので奨励会に入るまでがタダの棋士、奨励会時代は六級から一級まで（このシステムはよくわかっていない）、一段からの昇段日をまとめた。
全部で32人の藤井〇〇という称号があるようだった。
ただし、実際の生存解析では各個人は独立である。藤井竜王・名人は藤井竜王、藤井名人であることを含んでいるし、藤井N冠はN-1冠すべてを含んでいるので各人は独立ではない。
直近のタイトル戦では12月7日の竜王戦まで4連敗するまでは八冠のタイトルを保持できるから、2023年12月7日まで八冠と仮定した。
中央生存日数で296日は八冠でいられるようだ。
MSTの信頼区間の計算はBrookmeyer & Crowley の手法を使った。
https://www.jstor.org/stable/2530286
Kaplan-Meier Survival Estimates (Survival Curves) - StatsDirect
統計学入門−第11章
library(survival) library(prodlim) txt <- " class,title,start,end 出生,ama,2002-07-19,2012-09-21 六級,ama,2012-09-22,2012-11-09 五級,ama,2012-11-10,2013-05-04 四級,ama,2013-05-05,2013-06-01 三級,ama,2013-06-02,2013-09-14 二級,ama,2013-09-15,2014-04-04 一級,ama,2014-04-05,2014-06-20 一段,dan,2014-06-21,2015-02-27 二段,dan,2015-02-28,2015-10-17 三段,dan,2015-10-18,2016-09-30 四段,dan,2016-10-01,2018-01-31 五段,dan,2018-02-01,2018-02-16 六段,dan,2018-02-17,2018-05-17 七段,dan,2018-05-18,2020-08-19 八段,dan,2020-08-20,2021-07-02 九段,dan,2021-07-03, 棋聖,title,2020-07-16, 王位,title,2020-08-20, 叡王,title,2021-09-13, 竜王,title,2021-11-13, 王将,title,2022-02-12, 棋王,title,2023-03-19, 名人,title,2023-06-01, 竜王・名人,title,2023-06-01, 王座,title,2023-10-11, 二冠,kan,2020-08-20,2021-09-12 三冠,kan,2021-09-13,2021-11-12 四冠,kan,2021-11-13,2022-02-11 五冠,kan,2022-02-12,2023-03-18 六冠,kan,2023-03-19,2023-05-31 七冠,kan,2023-06-01,2023-10-10 八冠,kan,2023-10-11, " sysd <- Sys.Date() # sysd <- as.Date("2023-12-07") dat <- read.csv(text=txt, stringsAsFactors=FALSE) dat$end <- ifelse(dat$end == "", as.character(sysd), dat$end) dat$days <- as.numeric(as.Date(dat$end) - as.Date(dat$start)) dat$event <- (dat$end != sysd) + 0 fit <- survfit(Surv(days, event=event, type="right") ~ 1, data=dat) sfit <- prodlim(Hist(days, event=event) ~ 1, data=dat) ss <- summary(sfit) s1 <- ss$table[ss$table[, "n.event"] == 0, ] s2 <- merge(subset(dat[dat$days%in%s1[, "time"],], title=="title"), s1, by.x="days", by.y="time") # MST の信頼区間の計算 g <- switch((is.null(fit$strata))+1, unlist(mapply(rep, names(fit$strata), fit$strata)), rep(1, length(fit$n.event))) er <- tapply(fit$surv*fit$std.err, g, c) p <- tapply(fit$surv, g, c) t0 <- tapply(fit$time, g, c) i0 <- lapply(p, function(z) tail(which(z > pai), 1)) i1 <- lapply(p, function(z) head(which(z < pai), 1)) j0 <- mapply(function(z0, z1) z0[z1], t0, i0, SIMPLIFY=FALSE) j1 <- mapply(function(z0, z1) z0[z1], t0, i1, SIMPLIFY=FALSE) mst <- mapply(function(z0, z1, z2, z3){z0[z1]+(z0[z2]-z0[z1])*(z3[z1]-pai)/(z3[z1]-z3[z2])}, t0, i0, i1, p) cis <- do.call(rbind, list( mapply(function(z0, z1, z2, z3){z0[tail(which(z1+2*((z0>z2)-0.5)*qnorm(0.975)*z3 > pai & z0 < z2), 1)+1]}, t0, p, mst, er), mst, mapply(function(z0, z1, z2, z3){z0[head(which(z1+2*((z0>z2)-0.5)*qnorm(0.975)*z3 < pai & z0 > z2), 1)-1]}, t0, p, mst, er) ) ) half_time <- sfit$time[which(sfit$surv < 0.5)[1]] xl <- c(0, 365*5) axis1.at <- at.t <- floor(seq(0, ceiling(max(xl/365)), by=0.5)*365) axis1.labels <- sprintf(ifelse(at.t%%365 == 0, "%.f", "%.1f"), at.t/365) par(mar=c(4, 4, 1, 1), las=1, cex.lab=1.5) plot(sfit, lwd=3, xlim=xl, xlab="タイトル保持期間（年）", ylab="", axis1.at=axis1.at, axis1.labels=axis1.labels, atrisk.at=at.t, marktime=TRUE, background=TRUE, background.horizontal=NA, atrisk.labels="藤井棋士", confint.col="yellow") mtext("生存確率", 2, line=2, las=3, cex=1.5) p <- par()$usr # segments(half_time, 0.5, y1=p[3], lty=3) # segments(half_time, 0.5, x1=-100, lty=3) # text(half_time, p[3], sprintf("%d日", half_time), cex=2, adj=c(-0.1, -0.5)) arrows(cis[2, 1], 0.5, cis[1, 1], angle=90, length=0.1, col=4, lwd=3) arrows(cis[2, 1], 0.5, cis[3, 1], angle=90, length=0.1, col=4, lwd=3) segments(cis[2, 1], 0.5, y1=p[3], lty=3) txt <- sprintf("%.f日 [95%s CI: %.f ｰ %.f]", cis[2, 1], "%", cis[1, 1], cis[3, 1]) text(half_time, p[3], txt, cex=2, adj=c(0, -0.5)) txt <- sprintf("%s年%d月%d日時点", format(sysd, "%Y"), as.numeric(format(sysd, "%m")), as.numeric(format(sysd, "%d"))) legend("topright", legend=txt, bty="n", cex=2) # 生存確率が同じタイトルを同一x軸上にtext する xd <- 100 r <- rle(s2[, "surv"]) xx <- unlist(mapply(rep, s2[cumsum(r$lengths), "days"], r$lengths)) dup <- duplicated(s2[, "days"]) | duplicated(s2[, "days"], fromLast=TRUE) dif <- c(xd, diff(s2[, "days"])) tcex <- 1.2 for(i in seq(nrow(s2))){ x0 <- s2[i, "days"] y0 <- s2[i, "surv"] x1 <- xx[i] + xd if(i != 1){ # 王位と棋聖のように近いタイトルを上下離す flag <- 2*((!(dif < xd/2 & dif > 0))-0.5)[i] y1 <- y0 + (p[4]-p[3])*(x1-x0)/(p[2]-p[1])*flag } else { y1 <- y0 } if(!dup[i]){ txt <- sprintf("%s （%d日）", s2[i, "class"], s2[i, "days"]) segments(x0, y0, x1, y1, lty=3, lwd=2) text(x1, y1, txt, pos=4, cex=tcex) } else { if(duplicated(s2[, "days"])[i]){ # 竜王・名人がかぶる txt <- sprintf("%s （%d日）", paste(s2[dup, "class"], collapse="/"), s2[i, "days"]) segments(x0, y0, x1, y1, lty=3, lwd=2) text(x1, y1, txt, pos=4, cex=tcex) } } } d0 <- list(as.POSIXct(dat$start), as.POSIXct(dat$end)) interval <- list(seq(as.POSIXlt(format(min(d0[[1]]), "%Y-01-01")), max(d0[[2]]), by="year"), seq(as.POSIXlt(format(min(d0[[1]]), "%Y-01-01")), max(d0[[2]]), by="month") ) interval[[2]] <- setdiff(interval[[2]], interval[[1]]) y <- 1:nrow(dat) xl <- as.numeric(c(as.POSIXlt("2012-04-01"), max(d0[[2]]))) yl <- c(1, nrow(dat)) cols <- c(ama="red", dan="yellow", title="cyan", kan="orange") par(mar=c(2.5, 6, 2.5, 1.5)) plot(d0[[1]], y, type="n", xlim=xl, ylim=yl, xaxt="n", yaxt="n", xlab="", ylab="", bty="[") axis(2, at=y, labels=dat$class, las=1) axis(1, at=interval[[1]], labels=format(interval[[1]], "%Y")) axis(3, at=interval[[1]], labels=format(interval[[1]], "%Y")) abline(v=interval[[1]]) # abline(v=interval[[2]], lty=3) # 毎月 for(i in seq(nrow(dat))){ #segments(x0=d0[[1]][i], y0=i, x1=d0[[2]][i], lwd=20) px <- c(d0[[1]][i], d0[[2]][i])[c(1,2,2,1)] py <- i + c(-1, -1, 1, 1)*0.5 polygon(px, py, col=cols[dat$title[i]], border=NA) text(px[2], i, dat$day[i], pos=4, xpd=TRUE) } abline(h=c(0, y)+0.5, lty=3)
