パワプロ6で一番効率のよい変化球レベルの上げ方

パワポケ1はやり込んだが当時はN64を持っておらず、パワプロ6 は友達の家で遊びまくった。
当時の投手サクセスの変化球レベルアップのシステムは、経験値ではなく練習コマンドを何回かやると運がよければ変化球レベルが上がる。
変化球レベルが低く、かつ球種が少ないと上がりやすいが、たくさん球種を持つと全然上がらなかった。
当時はネットなど発達しておらず、全然情報がなかったが、内部的にはこんなアルゴリズムらしい。

いまとなってはゲームをする元気もなく、TAS動画を見ているが、パワプロのプレイ動画を上げる思い出ホイホイの投稿者がこんな動画を上げていた。
www.youtube.com
www.youtube.com

というわけで変化球レベルの効率のよい上げ方と、必要とされる練習回数を推定というかシミュレーションする。
サイトを参考に練習後における最大乱数と、実際にサンプリングされる乱数を一様分布から取ってくる。
投法（オーバーかサイド/アンダーか）とシンカーの関係はやってみてもいいが、単純化してシンカー以外の4方向をレベル上げしたらどうなるかを検討する。
5方向かつシンカーを考慮しても0-7レベル分の8段階とシンカーかどうかで計算量が16倍されるので暇があれば検討してもよい。
パワプロ５サクセス変化球講座基本編
パワプロ６サクセス変化球講座

シンカー以外の4方向は何からレベルを上げていっても数学的には対称性があるので、4方向を0から7まで上げるやり方は0から7の玉を4回重複を許して取り出す組み合わせなとなり 通りある。これらをすべてvertex として、1ポイント変化球レベルが上がったときの推移の仕方をグラフで考える。
このグラフについて変化球のパターン（vertex もしくはnode）とレベルが1ポイント上がったときの推移はエッジとなる。ここである変化球パターンからに推移するときのエッジについて、何回か練習をして溜まった変化球習得ポイントと、上記乱数を比較して変化球レベルが1ポイント上がるが、ここで変化球レベルが上がるのに何回（何週）練習が必要になるかは、確率分布になる。これを としておくと、各エッジで簡単にレベルが上がる（練習回数が少なくて済む）場合となかなかレベルが上がらない（練習回数がたくさん必要）場合とがあるので、これは重みになる。分布でもっていてもよいが、グラフに実装するには練習回数分布の平均値でも最頻値でもいいが としておく。
これで変化球のレベルの上がり方パターンのグラフを作成する。

グラフができたのでshortest_pathsで最適な（必要練習回数が最も少ない）レベルの上げ方のパスを取得する。
いでにクラスター検出もしてみる。レベルがあがって球種も増えるほど強そうなクラスターになる。
ここでは総変化量が21となる(7,7,7,0) と(6,5,5,5) を比較のために最短経路とノードを描出した。クラスター数は自動で決定され、この2つの変化球パターンとしては別グループに属するようだ。(6,5,5,5) のほうが変化球種類が多く複雑なのでおそらくこちらのグループのほうが高度というか投手力が強いグループと思われる。

Rとigraphを使ったネットワーク解析と可視化

グラフの描出にいくつかlayout 関数を使ってみたが、sugiyama というのが見栄えがよかった。これは階層グラフに対していい感じになるアルゴリズムらしく、総変化量と球種が段階的に増えていくので最適と思った。
https://igraph.org/r/html/1.2.6/layout_with_sugiyama.html
https://orsj.org/wp-content/corsj/or63-1/or63_1_20.pdf

総変化量が21のパターンについてどの変化量の組み合わせが最も効率的かを検討する。このパターンは、(7,7,7,0), (7,7,6,1), ..., と11通りある。総変化量21に達するまでに必要な週数の分布はこのようになる。

必要練習回数の最頻値のみを考慮すれば、(7,6,6,2) がもっとも早く99週で到達すると見込まれる。
一方で、すべての組み合わせについて分布を比較してみると、確率は引き分けがだいたい5%くらいあって、一方のパターンが勝つ（練習回数が少なくて済む）確率はせいぜい55%くらいなのでやる気の増減やあかつき大学一軍補正の50%の確率で+1、みたいなもので容易にブレる程度の勝率と考えられるので、有意にレベルが上がりやすい、とは言えないと思われる。(7,6,5,3) がすべての組み合わせで勝負して勝ち越すことが多いので、単なる分布でみれば(7,6,6,2) が平均（最頻値）がもっとも低いが、(7,6,5,3) のほうが分布の裾が中心に寄っているのでその影響で対比較をしたら(7,6,5,3) のほうが勝ちやすい（はやくレベルが上がりやすい）のだろうか。

別に総変化量21 のみにこだわらずすべての総変化量で検討すればよいと思ったのでやってみたが、一度にプロットするのは数が多すぎてできなかったので
20〜24については以下の通り。

すべてを無理やりプロット（boxplotにした）すると以下の通り。動画の後編の2:30くらいで(5,5,5,2) から(5,5,5,3) にレベルを上げるのに8週程度、と言っているのでだいたいよかろう。
もともとの解析リンクとは週数が微妙に異なっているので、もしかしたら変化球ポイントの(i) 現在の総変化量が11以上の場合はさらに2倍、(ii) 新しい球種を覚える場合はさらに1.5倍+10、の効果が重複する、のところで+10 の扱いが間違っているのかもしれない。

本当のサクセスはやる気によるポイントの上下や、リリーフ選択のタイミングなどでもっとばらつくだろうが、ここでは簡単に絶好調で、サクセス最初からリリーフ選択、ということにして条件を減らしている。
やる気があれば上記のシンカーと同様に、いつのタイミングでやる気を回復させるか、させないかなど検討できる。

# http://pawapuro56success.web.fc2.com/5henkakyuu.html
# http://pawapuro56success.web.fc2.com/6henkakyuu.html
# 2:31
# https://www.nemotos.net/igraph-tutorial/NetSciX_2016_Workshop_ja.html
yakuki <- 2
throw <- 0 # オーバースローは0 サイドかアンダーなら1

# 2回生（1年目）4月1週から、月4週、勝ち進めば最大4回生（3年目）11月1週目まで練習の機会がある
# 1年目11月2週にリリーフを希望する
Nweek <- c(list(4:12), list(1:12), list(1:11))
Week <- matrix(0, 4, length(unlist(Nweek)), dimnames=list(1:4, unlist(mapply(function(z) paste0(z, "_", Nweek[[z]]), seq(Nweek)))))

rensyu <- function(yaruki=0, throw=0, sinker=0, akatsuki=1, relief=1){
  # やる気 +4, +2, 0, -2, -4
  # オーバースロー 0 サイド・アンダー 1
  # オーバースローがシンカー練習をすると -1
  # あかつき一軍は +2, パワフル大学は50% の確率で -1
  # リリーフ希望だと50% の確率で +1
  ifelse(sinker == 1, 5 + 2*yaruki - 1 + throw, 5 + 2*yaruki) + akatsuki + sample(0:relief, size=1)
}

# 球種とレベル
henkakyu <- list(over=rep(0, 4), sinker=0)

# 最大乱数 = (1)基本乱数 + (2)球種補正 + (3)変化量補正
# (1)基本乱数
# 現在の総変化量が11以上の場合はさらに2倍
# 新しい球種を覚える場合はさらに1.5倍+10
# (2)球種補正
kyusyuhosei <- c("0"=0, "1"=0, "2"=0, "3"=15, "4"=32, "5"=55)
# (3)変化量補正
levelhosei <- c("0"=0, "1"=10, "2"=20, "3"=30, "4"=40, "5"=75, "6"=90)

# 実際のシミュレーションに必要なのはここから
# 変化球のレベル・球種パターン
e <- expand.grid(rep(list(0:7), 4))
pattern <- e[e[,1]>=e[,2] & e[,2]>=e[,3] & e[,3]>=e[,4],]
v <- do.call(rbind, mapply(function(z) cbind(pattern, z), 0:7, SIMPLIFY=FALSE))
# vt <- t(v)
vt <- t(pattern)

# ある変化球パターンのときから1レベルあげるときの最大乱数
M <- rep(0, nrow(e))
for(i in 1:nrow(e)){
  e0 <- vt[, e[i, 1]] # 前
  e1 <- vt[, e[i, 2]] # 後
  idx <- which(e0 != e1)
  ifnew <- (e0[idx] == 0)
  s <- sum(e0)
  maxransu <- c(
    kihon=(5*s) * ((s > 10) + 1) * ifelse(ifnew, 1.5, 1) + 10*ifnew,
    kyusyu=kyusyuhosei[sum(e0 > 0) + 1],
    level=levelhosei[e0[idx]+1]
  )
  M[i] <- sum(maxransu)
}

# シミュレーションする
maxtry <- 100
iter <- 1000
success <- matrix(0, nrow(e), iter)
for(i in 1:nrow(e)){
  for(j in 1:iter){
    r <- replicate(maxtry, rensyu())
    R <- sample(0:M[i], replace=TRUE, size=maxtry)
    success[i, j] <- which(cumsum(r) >= R)[1]
  }
}
# 平均（最頻値だが）的なレベルアップの乱数（週数）
med <- apply(success, 1, quantile, 0.5, type=1)

# あるパターンから1レベルあげるときの推移をグラフにする
library(igraph)
vmat <- matrix(0, ncol(vt), ncol(vt))
for(i in 1:ncol(vt)){
  tmp <- vt[, i]
  for(j in 1:5){
    idx <- which(apply(vt == replace(tmp, j, tmp[j]+1), 2, all))
    if(length(idx) > 0){
      vmat[i, idx] <- 1
    }
  }
}
diag(vmat) <- 0
lcols <- c("pink", rev(heat.colors(max(colSums(vt))))) # 総変化量
kcols <- c("pink", "darkcyan", "black", "darkcyan", "black") # 球種
shape <- c("circle", "circle", "circle", "square", "square")
ecols <- topo.colors(max(med))

g <- graph_from_adjacency_matrix(vmat)
layoutsugiyama <- layout_with_sugiyama(g)
e <- get.edgelist(g)

#V(g)$label <- NA
V(g)$label <- colSums(vt)
V(g)$label.cex <- 0.5
V(g)$color <- lcols[colSums(vt)+1]
V(g)$color[1] <- "pink"
V(g)$frame.color <- kcols[colSums(vt > 0) + 1]
V(g)$frame.width <- 1.5
V(g)$size <- 5
V(g)$shape <- c("circle", "square")[(colSums(vt >= 1) >= 4) + 1]
V(g)$shape <- shape[colSums(vt != 0)+1]
#E(g)$color <- gray(0.8)
E(g)$color <- ecols[med]
E(g)$width <- med/1.5
E(g)$arrow.size <- 0.3
E(g)$weight <- med

pch <- c(16, 21, 21, 22, 22)
par(mar=c(0, 0.2, 0.4, 2.5), xpd=TRUE)
plot(g, layout=layoutsugiyama, vertex.label.font=2, vertex.label.color="black")
legend("bottomleft", legend=c(0:4), horiz=FALSE, pch=pch, col=kcols, pt.cex=2, title="球種", pt.lwd=3)
legend("topleft", legend=seq(ecols), pch=15, pt.cex=3, col=ecols, title="", bty="n", xpd=TRUE, x.intersp=-0.3, inset=0)
legend("topleft", legend=NA, title="1レベル上がるのに平均的にかかる週数", bty="n")
legend("right", legend=seq(lcols)-1, pch=15, pt.cex=3, col=lcols, title="総変化量", bty="n", xpd=TRUE, x.intersp=0.05, inset=-0.10, adj=c(0.5, 0.5), xjust=-1, title.adj=-2)


i <- 316 # 5,5,5,5
i <- 119 # 7,7,6,0
i <- 317 # 6,5,5,5
i <- 120 # 7,7,7,0
spath <- shortest_paths(g, from=1, output="both", weights=med)
es <- mapply(function(z) as.numeric(as.character(z)), spath$epath) # 最短エッジ

E(g)$color[es[[i]]] <- "black"
V(g)$frame.width[spath$vpath[[i]]] <- 3
V(g)$size[spath$vpath[[i]]] <- 7
plot(g, layout=layoutsugiyama)

# ceb <- cluster_edge_betweenness(g, bridge=FALSE) # これは警告
ceb <- cluster_fast_greedy(as_undirected(g))
dendPlot(ceb, mode="hclust")

par(mar=c(0, 0.2, 0.4, 2.5), xpd=TRUE)
plot(ceb, g, layout=layoutsugiyama, edge.color=E(g)$color, col=V(g)$color, mark.border=NA, vertex.label.font=2)
legend("bottomleft", legend=c(0:4), horiz=FALSE, pch=pch, col=kcols, pt.cex=2, title="球種", pt.lwd=3)
legend("topleft", legend=seq(ecols), pch=15, pt.cex=3, col=ecols, title="", bty="n", xpd=TRUE, x.intersp=-0.3, inset=0)
legend("topleft", legend=NA, title="1レベル上がるのに平均的にかかる週数", bty="n")
legend("right", legend=seq(lcols)-1, pch=15, pt.cex=3, col=lcols, title="総変化量", bty="n", xpd=TRUE, x.intersp=0.05, inset=-0.10, adj=c(0.5, 0.5), xjust=-1, title.adj=-2)

# 総変化量21のノードをテキストで指示する
x <- c("(6,5,5,5)"=317, "(7,7,7,0)"=120)
for(i in seq(x)){
  xy <- norm_coords(layoutsugiyama$layout)[x[i],]
  text(xy[1], xy[2], names(x)[i], adj=c(1.3, 1.3), font=2)
}


# 総変化量の全パターンで解析
res <- NULL
ll <- 20:24 # 描出する量が多いので一部分にする
ll <- 1:28  # 全部描出する。barplot 
#for(l in 1:max(colSums(vt))){
for(l in ll){
  i21 <- which(colSums(vt) == l)
  y <- mapply(function(w) mapply(function(z) success[z,], es[[w]]), i21, SIMPLIFY=FALSE)
  y21 <- mapply(rowSums, y)
  colnames(y21) <- paste0("(", apply(vt[,i21, drop=FALSE], 2, paste0, collapse=","), ")")
  res <- cbind(res, y21)
}
lidx <- unlist(mapply(rep, seq(ll), table(match(colSums(vt), ll))))
vcols <- rainbow(length(ll))
ymed <- apply(res, 2, quantile, 0.5, type=1)

library(vioplot)
par(mar=c(4, 7, 1, 1))
vioplot(res, horizontal=TRUE, side="right", col=vcols[lidx], axes=FALSE, xaxt="n", yaxt="n", border=NA, lineCol=NA, rectCol=NA, xlim=c(2, length(lidx)-1))
axis(1, lwd=0, lwd.ticks=1)
#axis(2, lwd=0, lwd.ticks=1, at=1:ncol(res), labels=colnames(res), las=1)
abline(h=1:ncol(res), lty=3)
mtext("総変化量に達するまでにかかる練習回数（週数）", side=1, line=2.5)
mtext("変化球のパターンと総変化量", side=2, line=6.0)
text(ymed, 1:ncol(res), ymed)
box()

sw <- strwidth(colnames(res))[1]/2
sw0 <- 1.10
sh <- strheight(colnames(res))[1]/2
sh0 <- 1.50
xat <- par("usr")[1] - (par("mgp")[2] * diff(grconvertX(c(0, 1), from="lines", to="user")))
xat0 <- xat - sw
for(l in seq(ll)){
  y01 <- which(lidx == l)
  y01 <- y01[c(1, length(y01))]
  xpoly <- xat0 + c(-1, 1, 1, -1)*sw*sw0
  ypoly <- y01[c(1, 1, 2, 2)] + c(-1, -1, 1, 1)*sh*sh0
  polygon(xpoly, ypoly, xpd=TRUE, border=NA, lwd=2, col=adjustcolor(vcols[l], 0.4))
  text(xpoly[1], mean(ypoly), ll[l], xpd=TRUE, pos=2)
}
axis(2, lwd=0, lwd.ticks=1, at=1:ncol(res), labels=colnames(res), las=1) # 透過性の上から描出したほうが見栄えがよい


resmed <- tapply(apply(res, 2, median), lidx, quantile, 0.5, type=1)
par(mar=c(4, 7, 1, 1))
boxplot(res, horizontal=TRUE, col=vcols[lidx], axes=FALSE, border=NA, xlim=c(10, length(lidx)-10))
mtext("総変化量に達するまでにかかる練習回数（週数）", side=1, line=2.5)
mtext("変化球のパターンと総変化量", side=2, line=6.0)
axis(1, lwd=0, lwd.ticks=1)
#axis(2, lwd=0, lwd.ticks=1, at=1:ncol(res), labels=colnames(res), las=1)
abline(h=cumsum(rle(lidx)$lengths)+0.5, lty=3)
text(resmed, tapply(1:ncol(res), lidx, mean), resmed)
box()
# このあと上のpolygon を実行する


l <- 21
i21 <- which(colSums(vt) == l)
y <- mapply(function(w) mapply(function(z) success[z,], es[[w]]), i21, SIMPLIFY=FALSE)
y21 <- mapply(rowSums, y)
colnames(y21) <- paste0("(", apply(vt[,i21], 2, paste0, collapse=","), ")")
ymed <- apply(y21, 2, median)
vcols <- rainbow(length(i21))
names(vcols) <- colnames(y21)

library(vioplot)
par(mar=c(4, 6, 1, 1))
vioplot(y21, horizontal=TRUE, side="right", col=vcols)
axis(1, lwd=0, lwd.ticks=1)
axis(2, lwd=0, lwd.ticks=1, at=1:ncol(y21), labels=colnames(y21), las=1)
abline(h=1:ncol(y21), lty=3)
mtext(sprintf("総変化量%dに達するまでにかかる練習回数（週数）", l), side=1, line=2.5)
mtext("変化球のパターン", side=2, line=4.5)
text(ymed, 1:ncol(y21), ymed, pos=1)

wlmat <- matrix(0, length(i21), length(i21), dimnames=replicate(2, colnames(y21), simplify=FALSE))
for(i in 1:nrow(wlmat)){
  for(j in 1:ncol(wlmat)){
    wlmat[i, j] <- mean(y21[,i] < y21[,j])
  }
}

# 総当り結果
w <- graph_from_adjacency_matrix(wlmat > t(wlmat))
we <- get.edgelist(w)
wlayout <- layout_in_circle(w)
win <- rowSums(as.matrix(as_adjacency_matrix(w)))
V(w)$label <- paste0(colnames(y21), "\n", win)
V(w)$color <- vcols
V(w)$label.color <- "black"
V(w)$size <- 4+1.5*win
V(w)$frame.color <- "white"
E(w)$color <- vcols[we[,1]]
E(w)$width <- 3
E(w)$arrow.size <- 1.0
par(mar=c(0, 0, 0, 0))
plot(w, layout=wlayout, vertex.label.font=2)

w <- tapply(ymed, mapply(function(z) sum(eval(parse(text=z))), paste0("c", names(ymed))), c)
wx <- unlist(mapply(rep, as.numeric(names(w)), sapply(w, length)))
plot(wx, unlist(w), pch=16, las=1, xlab="総変化量", ylab="総変化量に達するまでにかかる練習回数（週数）")