グラフ構成の原則の学習記録

【学習動機】

"The Elements of Graphing Data"の第2章では, グラフを構成する際の原則が, グラフの例を挙げながら紹介されている.

https://www.amazon.co.jp/Elements-Graphing-Data-William-Cleveland/dp/0963488414

そこの原則を表す文をメモしておく.

【学習内容】

用語の定義

Fig 2.1とFig 2.2 を引用. これらの用語を知っておくと, 可視化ライブラリの関数の引数の意味が分かってきてよい気がする.

f:id:moratoriamuo271:20190118125451p:plain

f:id:moratoriamuo271:20190118125516p:plain

Clear Vision

  • Make the data stand out. Avoid superfluity. (データを目立たせ, 過剰は避けろ.)

  • Use visually prominent graphical elements to show the data. (視覚的に目立ったグラフの要素を使え.)

  • Use a pair of scale lines for each variable. Make the data rectangle slightly smaller than the scale-line rectangle. Tick marks should point outward. (上下左右に目盛線を.データの箱は目盛線の箱より僅かに小さくすべき. 目盛は外側に.)

  • Do not clutter the interior of the scale-line rectangle. (目盛線の箱の内側を散らかすな.)

  • Do not overdo the number of tick marks. (目盛を使いすぎるな.3~10個で十分.)

  • Use a reference line when there is an important value that must be seen across the entire graph, but do not let the line interfere with the data. (グラフ全体を横切らなければならない重要な値があるなら基準線を使え.ただし,データの邪魔はするな.)

  • Do not allow data labels in the interior of the scale-line rectangle to interfere with the quantative data or to clutter the graph. (量的データを邪魔したり,グラフを散らかすようなデータラベルを目盛の箱の内側に置くな.)

  • Avoid putting notes and keys inside the scale-line rectangle. Put a key outside, and put notes in the caption or in the next. (注釈や凡例を目盛線の箱の内側に置くのは避けよ. 凡例は外側に. 注釈は説明文か文中に.)

  • Overlapping plotting symbols must be visually distinguishable. (重なり合うプロットは見分けがつくようにすべし.)

  • Superposed data sets must be readily visually assembled. (重ね合わせるデータセットは視覚的に分かりやすく整理すべし.)

  • Visual clarity must be preserved under reduction and reproduction. (複製や縮小を経ても視覚的明確さが損なわれないようにすべし.)

Clear Understanding

  • Put major conclusions into graphical form. Make captions comprehensive and informative. (主要な結論をグラフの形式に落とし込むこと.説明文は包括的で有益に.)

  • Error bars should be clearly explained.(エラーバーは明確に説明されるべし.)

  • When logarithms of a variable are graphed, the scale label should correspond to the tick mark labels. (対数グラフを書くときは, スケールラベルを目盛ラベルと対応させるべし.)

  • Proofread graphs. (グラフは校正せよ.)

  • Strive for clarity. (今までの原則, 全部やれ.)

【学習予定】

色調が与える印象などについても学んでいきたいところ.

Banking to 45°の学習記録

 

【学習動機】

ggplot2の学習を進めていて, グラフの作り方は多少分かったが, グラフをどう作るべきかは分からなかったため, それが分かるかもしれないと思い, “The Elements of Graphing Data”を読むことにした (和訳が出ているのだが, 入手困難.).

www.amazon.co.jp

これを読み進めていて, banking to 45°という概念に出会った. この概念について言及している日本語文献があまりなく(奥村先生のグラフの描き方に関するページ

https://oku.edu.mie-u.ac.jp/~okumura/stat/090503b.html

くらいしか見つからなかった.), これについて, ggplot2を利用しつつ, 記録しておこうと思った.

【学習記録】

アスペクト比はグラフから変化率を見極めるのに重要で, 最適なアスペクト比を求めるのに, banking to 45°という方法がある.

今回使用するライブラリは以下のもの.

library(ggplot2)
library(ggfortify)
library(ggthemes)
library(tidyverse)

アスペクト比

アスペクト比は, 縦の長さ/横の長さ. 我々は曲線の局所的な線の方向から, 風速とオゾンの変化率に関する情報をデコードしている. グラフのアスペクト比が異なると, この方向が変わってしまい, 読み取れる情報が変わりうる.

太陽の黒点データ

太陽の黒点データを見てみる.

head(sunspot.year)
## [1]  5 11 16 23 36 58
autoplot(sunspot.year) + theme(aspect.ratio=1)

autoplot(sunspot.year) + theme(aspect.ratio=0.055)

両方とも約11年周期のサイクルがあることには気付けるが, 上のグラフでは黒点の増加の方が減少よりも急速であることに気付きにくい. さらに注意深くみると, ピークが高いサイクルではその傾向が顕著である一方, ピークが低いサイクルではその傾向が見えないということまで気付くことができる.

Banking to 45°

局所的な線の方向, 傾きの絶対値を45°に中心化させると, 線の方向を読み取るのに最適になるようだ. その中心化がbanking to 45°と呼ばれている. (和訳は分からないし, 思い付かない. 関係ないが, ななめ45°という芸人のトリオがいたことをこれで思い出した.)

具体的な計算式についてはここでは省略する.

ggthemesパッケージに最適なアスペクト比を計算してくれる関数がある.

ggthemesのドキュメントから引用(https://cran.r-project.org/web/packages/ggthemes/ggthemes.pdf)

Calculate the optimal aspect ratio of a line graph by banking the slopes to 45 degrees as suggested by W.S. Cleveland. This maximizes the ability to visually differentiate differences in slope. This function will calculate the optimal aspect ratio for a line plot using any of the methods described in Herr and Argwala (2006). In their review of the methods they suggest using median absolute slope banking (’ms’), which produces aspect ratios which are generally the median of the various methods provided here.

このドキュメントの, Exampleで, coord_fixed(ratio = ratio)で比率を変えているが, bank_slopesの産出はアスペクト比を出していると思われるので, theme(aspect.ratio=ratio)とすべきなのではないかと思った. 以下ではそのようにしている.

x <- seq_along(sunspot.year)
y <- as.numeric(sunspot.year)
ratio <- bank_slopes(x, y)
m <- ggplot(data.frame(x = x, y = y), aes(x = x, y = y)) +
  geom_line() + theme_bw()
m

m + theme(aspect.ratio = ratio)

co2のトレンドデータ

co2のデータのトレンド部分について. co2をdecomposeで分解する.

co2d <- decompose(co2)
plot(co2d)

x <- seq_along(co2d$trend)
y <- as.numeric(co2d$trend)
df <-data.frame(x = x, y = y) %>% na.omit()
ratio <- bank_slopes(df$x, df$y)
m <- ggplot(df, aes(x = x, y = y)) +
  geom_line() + theme_bw()
m + theme(aspect.ratio = ratio)

m + theme(aspect.ratio = 0.055)

上のbanking to 45°によるアスペクト比のグラフからは, 曲線の凸性が見て取れるが, 下のグラフからそれを読み取るのは難しい.

EDAをしていく上で可視化は欠かせない作業だが, アスペクト比に気を付けながらグラフを作成していくことが重要のようだ.

【学習予定】

グラフのデザインについてもう少し学んでいきたいところ.

終わるまでは終わらないよ

 オ……オオ……クロノ…… ナ……ナツカシイ……。 イヤ……、アナタ方にとってハ イッシュンの事だったのデスネ。 シカシ、ワタシニとっては 400年ハながい時間デシタ……。 シカシ、クロウのかいアッテ森ハよみがえりマシタ。 ……サア、今夜ハ、 400年ブリのサイカイをいわおうではアリマセンカ。

 
 冒頭の引用は、クロノ・トリガーというゲームのロボのセリフである。この140文字程のセリフを読むだけで、どこかこみ上げてくるものがある。

 今日、松濤美術館で開催されている「終わりのむこうへ : 廃墟の美術史」という展覧会に行ってきた。廃墟趣味に関して、18世紀頃のポンペイ遺跡の発掘やグランドツアー(上流階級の子弟が知見を広げるために行われた旅行)の流行による興隆、20世紀頃のシュルレアリスムによる変容、日本への伝播とその後が理解できる良い展示だった。実際の廃墟に訪れてみたいものだと思った。

 廃墟の寂れた風景に、廃れゆく建築物に、人は何を思うのか。時の流れで老いて朽ちる事実に哀愁を覚えるのか。それとも、今にも崩れようとしながらも聳え立とうとする姿に感情が動かされるのか。動きのない景色のはずなのに、ただ眺めているだけで、何かの目的で建築され使用されていたであろう遠い過去と、最期を迎えようとしているどこか近い未来の双方向への流れを感じさせられる。そういう点が廃墟の魅力の一つなのかもしれない。

 モノは、時間をかけ、エネルギーをかけて、作られる。この事実に自分が気付いたのは高校生になってからだった。高校の選択授業で美術を選んで油絵を何日もかけて描くことを経験するまで、絵というモノは上手い人がささっと描いてしまうものだと勘違いをしていた。よくできているモノは、それを作るのに(その着想に至るための経験等を含めて)、ほぼ間違いなく時間やエネルギーが費やされている。モノをただ見ているだけでは見えない、その奥行きを、積み重ねを感じ取れるかどうかは、大切なことのように思う。原価などといった一次元の情報や表層的で短期的な視野に囚われてしまって、コストや有用性を勘定するのは、もったいない。見える目を養っていかなければならない。

バドミントン戦略シミュレーションとの衝動遊戯:第3ゲーム

 

はじめに

「バドミントン戦略シミュレーションとの衝動遊戯:第2ゲーム」

moratoriamuo.hatenablog.com

で, 高速化を図ってRcppによる書き換えを試みたが上手くいかなかったことを書いたところ, ill-identifiedさんからコメントをいただき, コメントにしたがってRcppの書き換えたところ, 上手くいった. 以下ではRcppの書き換えと, それによる速度比較について書く.

Rcppによる書き換え

書き換えたコードが以下のものである.

#include <Rcpp.h>

using namespace Rcpp;

// [[Rcpp::export]]

NumericMatrix frallycpp(Function f1, Function f0, int ngame, double reffect){
  NumericMatrix mresult(6,0);
  NumericVector vresult(6);
  
  for(int game = 0; game < ngame; ++game){
    NumericVector vcrally(0);
    int p1point = 0;
    int p0point = 0;
    int player = 1;
    
    while(p1point < 21 && p0point < 21){
      int judge = 1;
      int shot = 0;
      int crally = 0;
      
      while(judge==1){
        double effect = (shot==0)? 0: reffect;
        shot = (player==1)? as<int>(f1(shot)): as<int>(f0(shot));
        double prob = (shot==0)? 0.9-effect: 0.7-effect;
        judge = (R::runif(0, 1) < prob)? 1: 0;
        player = (player==1)? 0: 1;
        crally = crally + 1;
      }
      vcrally.insert(vcrally.end(),crally);
      
      if(player==1){
        p1point = p1point + 1;
      }else{p0point = p0point + 1;}
    }
    vresult[0] = mean(vcrally);
    vresult[1] = var(vcrally);
    vresult[2] = max(vcrally);
    vresult[3] = p1point;
    vresult[4] = p0point;
    vresult[5] = (p1point>p0point)? 1: 0;
    
    mresult = cbind(mresult,vresult);
  }
  mresult = transpose(mresult);
  return mresult;
}  

前回詰まっていた2点のうち、一点は型の間違いによるもの、もう一点はRcppのrunif()とRのrunif()で返り値の型が違うことによるものだった. これらを修正し, 対戦する戦略関数とシミュレーション回数とrisky shotの効果を引数にとって、試合結果の行列を返す関数を作成した.

Rの関数についても同様に書き換えた.

frallyR <- function(fstr1,fstr0,ngame,reffect){
  #### vresultをrbindするmatrixの作成 ####
  mresult <- NULL
  vresult <- NULL
  
  #シミュレーション
  for(i in 1:ngame){
    #### game result 変数の初期値 ####
    vcrally <- c()    #rally count vector
    p1point <- 0    #player1のポイント
    p0point <- 0    #player0のポイント
    player <- 1    #first serverはplayer1
    
    #### rallyの定義とgameの定義 ####
    #rallyをjudge==0になるまでのwhileループとする.
    #gameをどちらかのplayerのpointが21になるまでのwhileループとする.
    
    while(p1point < 21 & p0point < 21){
      #ラリー前に初期化する
      judge <- 1    #shotの成功失敗の判定(1=>成功,0=>失敗)
      shot <- 0    #ラリー最初のサーブは成功率へのeffectは0なので相手のshotはsafetyとみなす.
      crally <- 0    #rally count
      
      
      while(judge==1){
        
        #手番のplayerの戦略関数をfstrに代入
        fstr <- ifelse(player==1,fstr1,fstr0)
        
        #相手のshotからeffectを決定
        effect <- ifelse(shot==0, 0, reffect)
        
        #自分のshotの決定  
        shot <- fstr(shot)
        prob <- ifelse(shot==0, 0.9-effect, 0.7-effect)
        
        #shotの成功判定
        judge <- ifelse(runif(1)<prob, 1, 0)
        
        #手番playerの交代
        player <- ifelse(player==1, 0, 1)
        
        #rally count追加
        crally <- crally + 1
      }
      
      #ラリーカウントベクトルの更新
      vcrally <- c(vcrally, crally)
      #ポイントの更新
      if(player==1){
        p1point <- p1point + 1 
      }else{p0point <- p0point + 1}
      
    }
    
    #試合結果(ラリー数平均,ラリー数分散,最大ラリー数p1point,p0point,gamewinplayer)
    vresult <- c(mean(vcrally),var(vcrally),max(vcrally),p1point,p0point,ifelse(p1point>p0point,1,0))
    mresult <- rbind(mresult,vresult)
  }
  return(mresult)
}

さらに, 速度比較の対象として並列計算するケースを用意するために, シミュレーションを繰り返さずに結果を返す関数も作成しておく.

fonerallyR <- function(fstr1,fstr0,reffect){

  vresult <- NULL
  
    #### game result 変数の初期値 ####
    vcrally <- c()    #rally count vector
    p1point <- 0    #player1のポイント
    p0point <- 0    #player0のポイント
    player <- 1    #first serverはplayer1
    
    
    #### rallyの定義とgameの定義 ####
    #rallyをjudge==0になるまでのwhileループとする.
    #gameをどちらかのplayerのpointが21になるまでのwhileループとする.
    
    while(p1point < 21 & p0point < 21){
      #ラリー前に初期化する
      judge <- 1    #shotの成功失敗の判定(1=>成功,0=>失敗)
      shot <- 0    #ラリー最初のサーブは成功率へのeffectは0なので相手のshotはsafetyとみなす.
      crally <- 0    #rally count
      
      
      while(judge==1){
        
        #手番のplayerの戦略関数をfstrに代入
        fstr <- ifelse(player==1,fstr1,fstr0)
        
        #相手のshotからeffectを決定
        effect <- ifelse(shot==0, 0, reffect)
        
        #自分のshotの決定  
        shot <- fstr(shot)
        prob <- ifelse(shot==0, 0.9-effect, 0.7-effect)
        
        #shotの成功判定
        judge <- ifelse(runif(1)<prob, 1, 0)
        
        #手番playerの交代
        player <- ifelse(player==1, 0, 1)
        
        #rally count追加
        crally <- crally + 1
      }
      
      
      #ラリーカウントベクトルの更新
      vcrally <- c(vcrally, crally)
      #ポイントの更新
      if(player==1){
        p1point <- p1point + 1 
      }else{p0point <- p0point + 1}
      
    }
    
    #試合結果(ラリー数平均,ラリー数分散,最大ラリー数p1point,p0point,gamewinplayer)
    vresult <- c(mean(vcrally),var(vcrally),max(vcrally),p1point,p0point,ifelse(p1point>p0point,1,0))
  return(vresult)
}

速度比較

必要なライブラリーの読み込みと戦略関数と上で定義した3つ関数の読み込みを行う.

library(Rcpp)
library(foreach)
library(doSNOW)
## Loading required package: iterators
## Loading required package: snow
source("strategy.R")
sourceCpp("frally.cpp")
source("frally.R")
source("fonerally.R")

単純にRで計算するケース, Rcppを使用して高速化するケース, foreachを使用して並列計算するケースの3つのケースで比較する. 設定としては, プレイヤー1がrisky戦略, プレイヤー0がsafety戦略をとり, シミュレーション回数は10000回, risky shotの効果の大きさを0.284とする.

単純にRで計算するケース

system.time(
  mr <- frallyR(fstr_allrisky,fstr_allsafety,10000,0.284)
)
##    user  system elapsed 
##   16.03    0.11   16.22

Rcppを使用して高速化するケース

system.time(
  mcpp <- frallycpp(fstr_allrisky,fstr_allsafety,10000,0.284)
)
##    user  system elapsed 
##   23.60    0.00   23.64

foreachを使用して並列計算するケース

cl <- makeCluster(4, type = "SOCK")
registerDoSNOW(cl)
system.time(
  mrp <- foreach(i=1:10000, .combine = "rbind") %dopar% fonerallyR(fstr_allrisky,fstr_allsafety,0.284)
)
##    user  system elapsed 
##    8.21    0.47   13.58
stopCluster(cl)

結果の確認

速度比較はできたが, 一応, 中で同様のことができているか確認する. 扱いやすいようにデータフレームの形式に変換する. mr以外の2つの行列についても同様.

dfr <- as.data.frame(mr)
colnames(dfr) <- c("mean","var","max","p1point","p2point","winp")
rownames(dfr) <- NULL
head(dfr)
##       mean       var max p1point p2point winp
## 1 3.175000  3.276282   9      19      21    0
## 2 3.973684 14.999289  19      17      21    0
## 3 4.121212  8.047348  13      21      12    1
## 4 3.684211  7.411095  11      21      17    1
## 5 3.390244  8.343902  12      20      21    0
## 6 2.825000  3.276282   9      19      21    0
dim(dfr)
## [1] 10000     6
colMeans(dfr[,c("winp","mean","var","p1point","p2point")])
##      winp      mean       var   p1point   p2point 
##  0.495700  3.328941  5.794026 18.440500 18.600500
head(dfcpp)
##       mean       var max p1point p2point winp
## 1 3.171429 11.969748  21      14      21    0
## 2 3.000000  3.600000   8      20      21    0
## 3 3.625000  5.596774   9      21      11    1
## 4 3.400000  4.041026  10      21      19    1
## 5 3.051282  5.418354  11      18      21    0
## 6 3.317073  5.671951   9      21      20    1
dim(dfcpp)
## [1] 10000     6
colMeans(dfcpp[,c("winp","mean","var","p1point","p2point")])
##      winp      mean       var   p1point   p2point 
##  0.495600  3.331092  5.813098 18.517200 18.557800
head(dfrp)
##       mean       var max p1point p2point winp
## 1 3.390244  4.693902   8      20      21    0
## 2 3.111111  3.473016   9      21      15    1
## 3 3.151515  5.570076  12      21      12    1
## 4 4.027778 14.027778  16      15      21    0
## 5 2.536585  2.004878   6      21      20    1
## 6 3.717949 10.628880  15      18      21    0
dim(dfrp)
## [1] 10000     6
colMeans(dfrp[,c("winp","mean","var","p1point","p2point")])
##      winp      mean       var   p1point   p2point 
##  0.488000  3.328949  5.789926 18.462100 18.596600

おそらく出力に問題はないと思われる.

おわりに

昨年学んだ並列計算で高速化できたのは良かったけど, Rcppで書き換えたのに、むしろ遅いってどういうことなの……. コードの書き方がおそらく非効率なのだろう. コードの最適化やC++について今後少しずつ学んでみようと思う.

<追記: 2019/01/02>

C++側で, Rで定義した関数を呼び出すところがどうしても遅くなるようで,  高速化をしたいのなら, それは避けるべきみたい. C++で定義し直してやるとよい. 

加えて, cbindはサイズが大きくなると遅くなるので, mresultのmatrixを先に用意して代入していく方がおそらく速くなる.

Twitter有識者(鍵アカウントなので名前の公開は避けます.)の方々のご指摘で気付きました. ありがとうございます.

<追記おわり>

 

 

Rcppとforeachという持ち札が増えたのは良かった. Rcppのコードで詰まっていたところについてのコメントは大変助かりました. ill-identifiedさん, ありがとうございました.

2018年の後日譚々

 2019年になったので、2018年を振り返る。2018年は、「きっと何者にもなれないお前たちに告げる」と「私、特別になりたいの」という二つの対立したセリフが心を貫いた一年間であった。片方は、『輪るピングドラム』のセリフで、このアニメは昨年初めて見たのだが、衝撃的だった。見終わったとき何が何だか分からなかったのだが何かがすごいということは感じた。そんな直感の理由を探るべく解釈を探して読み漁ったりした。村上春樹作品もいくつか読んだし、『銀河鉄道の夜』とその関連作品も漁ったりした。謎は謎のまま残っていてこれからも探索していきたい。もう片方は、『響け!ユーフォニアム』のセリフだ。『リズと青い鳥』を映画館で見た機会に見直して、改めて表現の仕方含め面白かった。ちなみに、『リズと青い鳥』も映画館で二回見た。

 この二つのセリフはとりわけ印象に残った。就職活動を通じて、自己存在の確認をする上での、相対する命題のようなものとして、頭のなかに存在していた。それらの思考の欠片は、その頃の記事に書いたように思う。思考はしてみたものの、何か明快な理解が得られたわけではない。きっと、これは死ぬまでずっと続く問いなのだろう。ただ一歩一歩足を進めることが重要であるし、意味のあることだと分かった一年であった。この3,4年間コツコツやっていたデータ分析でもバドミントンでも、けっして最高の結果が得られたわけではないないが、自分の成長が大きく見られた。自分の歩幅は天才や巨人のそれと比較して大きくないのは分かっているので、一歩一歩着実に進めていこうという認識・覚悟が確かなものになった。

 そんな覚悟のもとで、今年は何をしていこうか。統計周りの話は当然として、今年はとくにデザインやプレゼンについて、学んでいきたいところだ。圧倒的に知識や経験が不足している。学習以外にも対外的に、どこか積極的に外に出ていきたい。昨年は、サッカー観戦や六大学野球観戦に行ったり、TokyoRやJapanRなど勉強会に参加してみたり、技術書典に行ったりした。今年も、まだやっていないことに挑戦していきたい。あと、英語。昨年は英語のできなさでとても苦しんだ。一年通して改善していきたい。

 今年も良い年になりますように。

バドミントン戦略シミュレーションとの衝動遊戯:第2ゲーム

 

「バドミントン戦略シミュレーションとの衝動遊戯:第1ゲーム」で書いた内容の補足的な記事です。

moratoriamuo.hatenablog.com

 

戦略関数について

戦略関数のコードを以下に掲載する。6種類の戦略関数を定義して、strategy.Rという名前で保存している。自分で戦略関数を付け加えれば、シミュレーションすることが可能となる。

#戦略関数リスト
list_fstr <- c("fstr_allsafety", "fstr_allrisky", "fstr_tit4tat", "fstr_tit4tat_inv","fstr_trigger", "fstr_random")

#### fstr_allsafety ####
#全てsafety shotを打つ戦略関数
fstr_allsafety <- function(shot){
  shot <- 0
  return(shot)
}

#### fstr_allrisky ####
#全てrisky shotを打つ戦略関数
fstr_allrisky <- function(shot){
  shot <- 1
  return(shot)
}

#### fstr_tit4tat ####
#相手のshotと同じshotを打つ戦略関数
fstr_tit4tat <- function(shot){
  return(shot)
}

#### fstr_tit4tat_inv ####
#相手のshotと違うshotを打つ戦略関数
fstr_tit4tat_inv <- function(shot){
  shot <- ifelse(shot==1, 0, 1)
  return(shot)
}

#### fstr_trigger ####
fstr_trigger <- function(shot){
  shot <- ifelse(shot==1, 1, 0)
  return(shot)
}

#### fstr_random ####
fstr_random <- function(shot){
  shot <- sample(c(0,1),1)
  return(shot)
}

高速化について

条件分岐とループ処理がほとんどを占めるので、Rだとやはり遅い。 Rcppパッケージを使ってC++によって高速化したかったのだが、担当日までにコードの書き換えが上手くできなかった。 二か所上手くいかなかったc++のコードの途中作成物を以下に掲載する。

#include <Rcpp.h>

using namespace Rcpp;

// [[Rcpp::export]]

List frally(Function f1, Function f0){
  int ngame = 1000;
  double reffect = 0.284;
  
  NumericMatrix mresult(ngame, 6);
  NumericVector vresult(6);
  
  for(int game = 0; game < ngame; ++game){
    NumericVector vcrally(0);
    int p1point = 0;
    int p0point = 0;
    int player = 1;
    
    while(p1point < 21 && p0point < 21){
      int judge = 1;
      int shot = 0;
      int crally = 0;
      
      while(judge==1){
        Function fstr = (player==1)? f1: f0;
        double effect = (shot==0)? 0: reffect;
        shot = fstr(shot);   //ここでassigning to 'int' from incompatible type type 'SEXP'と怒られる
        double prob = (shot==0)? 0.9-effect: 0.7-effect;
        judge = (runif(1) < prob)? 1: 0;  //ここも怒られる
        player = (player==1)? 0: 1;
        crally = crally + 1;
      }
      vcrally.insert(vcrally.end(),crally);
      
      if(player==1){
        p1point = p1point + 1;
      }else{p0point = p0point + 1;}
    }
    //vresult <- c(mean(vcrally),var(vcrally),max(vcrally),p1point,p0point,ifelse(p1point>p0point,1,0))
    //mresult <- rbind(mresult,vresult)
  }
  return 0;
}  

Rの関数(定義した戦略関数)2つを引数にとって試合結果をリストとして出力する関数のコードを書こうとしたが、1ラリーをシミュレーションしているところ(while(judge==1){}のループ処理のところ)が上手くいかなくて躓いた。戦略関数をC++で書いて、includeする方法も試したが、同じように上手くいかなかった。 もう少し、C++とRcppについて学んでみようと思う。

Rの高速化として、並列化について、今年の秋に『Rによるハイパフォーマンスコンピューティング』を読んで学んだので、そちらを試しても良かった。10000回のシミュレーションについては各々独立した処理なので、そんなに苦労せず書けただろうが、結局やらなかった。

バドミントンついて

今回、シミュレーションをするために、バドミントンのかなりの部分を捨象して、モデル化を行った。プレー中、自分のショットの成功確率やその効果は、自分の練習量や疲労度合によって変わるし、相手のショットの種類も初めて戦う相手だと全く分からない。不確実な状況下で、プレーを行う。そのあたりにもまた奥深さがあるのだが、今回は紹介しきれなかった。紹介しきれなかったことや、バドミントンの上達方法やプレーについて考えたことは、また今度記事にする予定でいる。

バドミントン戦略シミュレーションとの衝動遊戯:第1ゲーム

 

こちらはみゅーもり Advent Calendar 2018の21日目の記事になります.

https://adventar.org/calendars/3429

21日目ということで, 21点のスポーツ, バドミントンを題材にします.

【ABSTRACT】

バドミントンは瞬間瞬間でリスクコントロールをしていく一面のあるスポーツである. そのリスクコントロールの一面に注目してモデル化を行い, シミュレーションによってショットの選択戦略を評価した. リスキーなショットはリスクに見合う効果が無いのならば打たない方が良いということと, 相応の効果がある場合は攻め時と守り時でリスクコントロールをしていくと良いということが確認できた.

【OBJECTIVE】

モチベーション

この記事を書くモチベーションは二つある.

一つは, バドミントンをあまりよく知らない人に, バドミントンの奥深さを知ってもらうことである. 今年はバドミントンの知名度が大きく上がった1年になったように思える。世界選手権で日本男子初の金メダルを獲得した桃田選手を初め、数多くの日本人選手が活躍した.『はねバド!』というバドミントン部の高校生が主役のアニメも放映された. これらをきっかけにバドミントンに少しでも興味を持った人に, もう少しバドミントンをよく知ってもらいたいという気持ちがある.

以下は男子ダブルスのラリーの動画である.

https://www.youtube.com/watch?v=PCO23zl9vaA

このような素早いラリー(スマッシュの初速のギネス記録は493km/hなんだとか)の中で, プレイヤーが意識的に・無意識に、緩急やコースの打ち分け,ショットの種類について, 戦略を立て予測・判断してプレーしている. このことが分かると観戦するときの楽しみ方が増えるかもしれない.

もう一つは, バドミントンの中級者に上がりたい初級者に, リスクコントロールの大切さを伝えることである. バドミントンを始めて一通りのショットが打てるようになると,  スマッシュなど一見,強いショットをたくさん打ちたくなる. 体勢が悪いところからでもそのようなショットを打ってしまいたくなる. 決まると楽しいからである.しかしながら,  そのようなショットは得てして成功率が低く,  多用してばかりいると試合でなかなか勝てない. そんな初級者に, 私が先輩に教わったり, 『ベイビーステップ』や, もこう先生の厨ポケ狩り講座, 遊戯王などから学んだリスクコントロールの面白さ, 強さを伝えたい. 良い体勢からは有効打が打ちやすいし, 悪い体勢からは打ちにくい. 悪い体勢で有効打が打てないわけではないが成功する確率は当然低くなる.  このへんのショット選択の駆け引きやリスクコントロールが, 勝敗の鍵を握っている. このリスクコントロール術は自分より格下の相手に負けにくくすることであると同時に,  自分より格上の相手に勝つ可能性を高めることでもある.  直感的にやっていって言語化してしにくいことを,  人に分かりやすく示すためにシミュレーションで確認する.

ルールの紹介とモデルのための簡単化

分からない人のためにも, 現行のバドミントンのルールの確認を行う. サービスを打ってラリーが始まると, 自分の打ったシャトルが相手のコートの中に落ちるか, 相手の打ったシャトルが自分のコートの外に落ちたら, 得点が1点入る. 得点をとったプレイヤーが次のサービスを打つ. ラリーを繰り返し, 21点を先にとった方が1ゲームを取ることができる. 20-20となった場合はそれ以降,点差が2点になるか,どちらかが30点をとるまでゲームが続く. 2ゲーム先取でその試合の勝者になる.  シミュレーションでは簡単化のため,  デュースはルールに入れず21点を先取したら1ゲーム取れることにしている. さらに, 1ゲーム取れば勝ちとする.

【METHODS】

モデルのセッティング

バドミントンのリスクコントロールの要素の側面を抽出したモデルによって, シミュレーションを行う.

一番基本的なモデルのイメージは以下のようなものである. 1ラリーの流れを示している. 

 

f:id:moratoriamuo271:20181220165841p:plain

ラリーのイメージ

 

まず, player1がサービスを打つ. このときsafety shotかrisky shotを選ぶ. safety shotは成功率が0.9と高いが, 次の相手のshotの成功率に影響を与えない(reffect=0). 一方, risky shotは成功率が0.7と低いが, 次の相手のshotの成功率を0.2低下させる(reffect=0.2). つまりリスクはあるが, 効果のあるshotということだ. shotを選んだら, 判定のフェイズに入る. 成功率 - reffect でshotが成功する(サービス時点では, 相手のshotに依存しないのでreffect=0). 成功したら, player0の手番になり, 同じことを行う. 失敗したら, player0の得点となり, 次は得点したplayer0のサービスでラリーが始まる. (念のために述べておくが, 現実では審判がシャトルのジャッジをするのは地面に落ちたときのみである.) このラリーをどちらかが21点とるまで繰り返すのが1ゲームである. 以下のコードがそれを実行するコードである.

#### strategy functionの定義(戦略関数は別で用意して読み込む) ####
source("strategy.R")

#### playerの戦略関数の決定 ####
list_fstr    #戦略関数リストから選択
fstr1 <- fstr_allsafety
fstr0 <- fstr_allrisky


#### game result 変数の初期値 ####
vcrally <- c()    #rally count vector
p1point <- 0    #player1のポイント
p0point <- 0    #player0のポイント
player <- 1    #first serverはplayer1


#### rallyの定義とgameの定義 ####
#rallyをjudge==0になるまでのwhileループとする.
#gameをどちらかのplayerのpointが21になるまでのwhileループとする.

while(p1point < 21 & p0point < 21){
  #ラリー前に初期化する
  judge <- 1    #shotの成功失敗の判定(1=>成功,0=>失敗)
  shot <- 0    #ラリー最初のサーブは成功率へのeffectは0なので相手のshotはsafetyとみなす.
  crally <- 0    #rally count
  
  
  while(judge==1){
    
    #手番のplayerの戦略関数をfstrに代入
    fstr <- ifelse(player==1,fstr1,fstr0)
    
    #相手のshotからeffectを決定(ここではrisky shotの効果を0.2としている)
    effect <- ifelse(shot==0, 0, 0.2)
    
    #自分のshotの決定  
    shot <- fstr(shot)
    prob <- ifelse(shot==0, 0.9-effect, 0.7-effect)
    
    #shotの成功判定
    judge <- ifelse(runif(1)<prob, 1, 0)
    
    #手番playerの交代
    player <- ifelse(player==1, 0, 1)
    
    #rally count追加
    crally <- crally + 1
  }
  
  
  #ラリーカウントベクトルの更新
  vcrally <- c(vcrally, crally)
  #ポイントの更新
  if(player==1){
    p1point <- p1point + 1 
  }else{p0point <- p0point + 1}
  
}

#試合結果(ラリー数平均,ラリー数分散,最大ラリー数p1point,p0point,gamewinplayer)
vresult<- c(mean(vcrally),var(vcrally),max(vcrally),p1point,p0point,ifelse(p1point>p0point,1,0))

これで, 最も基本的なモデルによる1ゲームが定義された. (strategy.Rは戦略関数を定義している. これについては次に中身を説明し, 別記事でコードを掲載する.) ゲームを10000回繰り返すことで各戦略の勝率(ゲーム取得率)を計算していく. 上のモデルでは, player1もplayer0も同じ能力を持っている. shotの成功率とreffectの大きさが同じである.

準備した戦略

戦略関数は一つ前の相手のshotを入力してsafety shotを打つかrisky shotを打つかを返す関数である(サービス時点では一つ前の相手のshotはsafetyと定義している.). 今回は, 以下の6種類を用意している.

## [1] "fstr_allsafety"   "fstr_allrisky"    "fstr_tit4tat"    
## [4] "fstr_tit4tat_inv" "fstr_trigger"     "fstr_random"
  • fstr_allsafety: 相手が何を打ってこようとsafety shotを打つ戦略. safety戦略と呼ぶ.
  • fstr_allrisky: 相手が何を打ってこようとrisky shotを打つ戦略. risky戦略と呼ぶ.
  • fstr_tit4tat: 一つ前の相手の打ってきたshotと同じshotで返す戦略. いわゆる, しっぺ返し戦略.
  • fstr_tit4tat_inv: 一つ前の相手の打ってきたshotと違うshotで返す戦略. 逆しっぺ返し戦略または天邪鬼戦略と呼ぶ.
  • fstr_trigger: 一度相手がrisky shotを打ってきたら, こちらはずっとrisky shotを打ち続ける戦略. いわゆる, トリガー戦略.
  • fstr_random: 相手が何を打ってこようと確率0.5でsafety, 確率0.5でriskyを選択する戦略. ランダム戦略と呼ぶ.

【RESULTS】

基本モデル(reffect=0.2の場合)

・safety vs safety

上段はplayer1の勝率, ラリー回数の平均と分散, player1の平均得点, player2の平均得点であり, 下段はラリー回数の最大値である.

##    winp    mean     var p1point p2point 
##  0.5016 10.0062 90.1300 18.5416 18.5224
## [1] 113

player1の勝率は 0.5016なので, 互角の勝負をしている. 互いにsafetyな戦略で戦っているのでラリー回数の平均は 10.0062 回, 最大ラリー回数は113回と長いラリーを続けている. ちなみに、2013年の世界選手権大会で, ラリー数108打の2分も続いたラリーがあった. 詳しくは以下の記事と動画を見てほしい.

https://rocketnews24.com/2014/03/25/425720/

・risky vs risky

上段はplayer1の勝率, ラリー回数の平均と分散, player1の平均得点, player2の平均得点であり, 下段はラリー回数の最大値である.

##    winp    mean     var p1point p2point 
##  0.4996  2.3988  2.2404 18.5287 18.5728
## [1] 19

player1の勝率は0.4996なので, 互角の勝負をしている. 互いにriskyな戦略で戦っているのでラリー回数の平均は2.3988回, 最大ラリー回数は19回と, safety vs safetyより断然短いラリーで終わっていることが分かる.

・safety vs risky

上段はplayer1の勝率, ラリー回数の平均と分散, player1の平均得点, player2の平均得点であり, 下段はラリー回数の最大値である.

##    winp    mean     var p1point p2point 
##  0.7746  3.7063  8.0190 20.1195 16.1272
## [1] 35

player1の勝率は0.7746であり, risky戦略の分の悪さが露骨に出ている. 現実でも, やはり相手が勝手にミスってくれると戦っていて楽である. ただし, この結果はrisky shotの成功率を下げる効果が0.2とあまり強くない(risky shotとsafety shotの成功率の差分と等しい)のが原因である. 後でこの効果をいじって, リベンジマッチを行う.

以降, player1にはランダム戦略のみで戦ってもらい, player0に他の戦略で戦ってもらって, 様子を見てみる.

・ランダム vs safety

上段はplayer1の勝率, ラリー回数の平均と分散, player1の平均得点, player2の平均得点であり, 下段はラリー回数の最大値である.

##    winp    mean     var p1point p2point 
##  0.2888  5.4007 21.6333 16.7775 19.8119
## [1] 62

player1の勝率は0.2888である. このセッティングではsafetyが圧倒的に強く,safetyには勝てない.

・ランダム vs risky

上段はplayer1の勝率, ラリー回数の平均と分散, player1の平均得点, player2の平均得点であり, 下段はラリー回数の最大値である.

##    winp    mean     var p1point p2point 
##  0.6191  2.9211  4.1788 19.2997 17.6817
## [1] 29

player1の勝率は0.6191である. やはり,このセッティングにおけるrisky戦略の弱さがうかがえる.

・ランダム vs しっぺ返し

上段はplayer1の勝率, ラリー回数の平均と分散, player1の平均得点, player2の平均得点であり, 下段はラリー回数の最大値である.

##    winp    mean     var p1point p2point 
##  0.4961  4.2780 10.7494 18.4860 18.5934
## [1] 40

player1の勝率は0.4961である. ほぼ互角の戦いをしている.

・ランダム vs 天邪鬼

上段はplayer1の勝率, ラリー回数の平均と分散, player1の平均得点, player2の平均得点であり, 下段はラリー回数の最大値である.

##    winp    mean     var p1point p2point 
##  0.5062  3.3518  7.2292 18.6457 18.5635
## [1] 35

player1の勝率は0.5062である. ほぼ互角であることが分かる.

・ランダム vs トリガー

上段はplayer1の勝率, ラリー回数の平均と分散, player1の平均得点, player2の平均得点であり, 下段はラリー回数の最大値である.

##    winp    mean     var p1point p2point 
##   0.495   4.278  10.731  18.532  18.526
## [1] 47

player1の勝率は0.495である. これまた,ほぼ互角か.

以上のシミュレーションではsafety shotを打ち続けるのが支配的な戦略となっている. これは, risky shotの効果がそのリスクに見合うものでなかったからだと分かる. したがって, 以下ではrisky shotの効果を上昇させて, safetyとriskyが釣り合うようにして, 戦略の強さを確認していく.

リベンジモデル(reffect=0.284の場合)

・リベンジ safety vs risky

safety vs riskyで, risky戦略は分の悪い結果となったが, risky shotの効果を0.2から0.284まで強化してみる.

上段はplayer1の勝率, ラリー回数の平均と分散, player1の平均得点, player2の平均得点であり, 下段はラリー回数の最大値である.

##    winp    mean     var p1point p2point 
##  0.5057  3.3521  5.8394 18.6062 18.5297
## [1] 30

player1の勝率は0.5057であり, 互角の勝負となっている. これでsafetyとriskyが釣り合った状態になった. これ以上, risky shotの効果を高めるとrisky shotが支配的な戦略となってしまうので, risky shotの効果を0.284にして他の戦略を見ていく.

・リベンジ ランダム vs しっぺ返し

上段はplayer1の勝率, ラリー回数の平均と分散, player1の平均得点, player2の平均得点であり, 下段はラリー回数の最大値である.

##    winp    mean     var p1point p2point 
##  0.5995  3.9948  8.7118 19.1948 17.7721
## [1] 39

player1の勝率は0.5995である. しっぺ返し戦略の方が分が悪いことが分かる. これは相手のrisky shotに対してrisky shotで挑んでしまうため, 相手のrisky shotが成功した場合, ミスしやすいという状況に陥ってしまう一方で, 相手がsafety shotを打っていても, こちらはrisky shotで攻めないため, 相手の安定性を崩せていないという解釈ができる.

・リベンジ ランダム vs 天邪鬼

上段はplayer1の勝率, ラリー回数の平均と分散, player1の平均得点, player2の平均得点であり, 下段はラリー回数の最大値である.

##    winp    mean     var p1point p2point 
##  0.4062  3.0051  5.0159 17.8443 19.1606
## [1] 29

player1の勝率は0.4062である. 天邪鬼戦略の方が強いことが分かる. これは相手のrisky shotが成功したら, safety shotで安全に返し, 相手がsafety shotを打って来たら, risky shotで攻めに行ってるのが功を奏しているのだろう. 攻め時と守り時のリスクコントロールがよくできていると解釈できる.

・リベンジ ランダム vs トリガー

上段はplayer1の勝率, ラリー回数の平均と分散, player1の平均得点, player2の平均得点であり, 下段はラリー回数の最大値である.

##    winp    mean     var p1point p2point 
##  0.6063  3.9975  8.7577 19.2244 17.7147
## [1] 39

player1の勝率は0.6063である. しっぺ返し戦略同様分が悪い. 解釈も同様である.

1ゲーム21点先取でない場合

現行のルールでは1ゲーム21点先取となっているのだが, 11点先取のルールが新ルールとして,今年提案されていたのだが,この案は却下された.11点だった場合,1ゲームを先取する場合,最適なリスクコントロール戦略は変わってくるだろう.おそらく,よりリスキーな戦略の方が有効になるのではないだろうか.

・11点制 リベンジ safety vs risky

risky shotの効果を0.284まで強化した状態で, safety戦略とrisky戦略で戦わせてみる. 予想では, risky戦略の方が勝率が高くなるのではないだろうか.

上段はplayer1の勝率, ラリー回数の平均と分散, player1の平均得点, player2の平均得点であり, 下段はラリー回数の最大値である.

##    winp    mean     var p1point p2point 
##  0.4982  3.3655  5.8964  9.2070  9.2109
## [1] 30

player1の勝率は0.4982であり, 予想と異なり, 21点のときと結果は変わっていない.

・271点制 リベンジ safety vs risky

逆に, 21点制より多い271点制を試してみる.

上段はplayer1の勝率, ラリー回数の平均と分散, player1の平均得点, player2の平均得点であり, 下段はラリー回数の最大値である.

##     winp     mean      var  p1point  p2point 
##   0.5105   3.3345   5.7968 262.5071 261.8643
## [1] 35

player1の勝率は0.5105であり, これまた, 21点のときと結果は変わらない.

以上, 実験してみたように, このセッティングでは, 点数の制度を変更しても勝率(ゲーム取得率)はあまり変わらなかった.

【CONCLUSIONS】

得られた結論

リスキーなショットは, リスクを背負うだけの強さがない限り(今回で言えば,reffectが0.2のケース)は使うべきでないことがよく分かった. 以下の動画のように, 上手く決まれば効果は絶大だが, 成功確率がそれに見合わないなら多用しない方が良いのである.

https://www.youtube.com/watch?v=eX9r9lTEbqw

一方でリスキーなショットにそれなりの効果がある場合(今回で言えば,reffectが0.284のケース)は, 天邪鬼戦略のように, 相手のリスキーなショットには安全なショットで守り, 相手の安全なショットにはリスキーなショットで攻めてリスクをコントロールするのが有効な戦略であることが確認された。

あることを失敗したら相手の得点になるという要素が濃いスポーツでは、このような戦略が最適になるのではないだろうか。逆に、あることを成功したら自分の得点になるという要素が濃いスポーツもある. この2つの要素の濃淡でリスク管理の質が変わってくると予想される.

バドミントンは相手のコートに返球し続ければ理論上負けない競技である. どんなに素晴らしいショットを打とうが芸術点は入らないので, 得点したとしても1点は1点であり, 自分からミスをして失点すると負ける競技なのである. だからといって, リスキーなショットが不要かというとそうではない. 相手にミスらせないと勝てないのである. 相手がミスをしないのならば, リスキーなショットを混ぜ合わせて得点しにいかないといけない. ミスを減らしたうえで そのような駆け引きができるようになるとバドミントンの強さが上がっていく. リスキーなショットは決まると楽しいので, つい使ってしまいがちだが, 強さを求めるのならば, リスキーなショットの使いどころは考えた方が良いだろう. ただし, 全く使わないのではなく, 有効な場面で使うことが重要である. また, 日々練習を繰り返し, 各ショットの成功率を上げることと, ショットの効果を高めることが重要である.

モデルの拡張可能性

戦略関数やモデルの設定を変更すれば容易に拡張可能である. 現在のモデルでは, 一つ前の相手のshotのみ記憶していて, それに対して反応を決めることしかできない. 一つ前相手のショットのみ記憶できるが、ショットの割合の記憶(最初は平坦な事前分布にしてベイズの公式で更新していく)するケースやショットの履歴を完全記憶できるケースに拡張でき, もっと複雑な戦略関数でシミュレーションすることができる. 現実のところは覚えているのはだいたいの割合だろう. 割合から予測を立てている気がする. 逆にそれを考えたうえで, 勝負をかけたいところで今まで打っていなかったshotを打つこともある. 割合を用いた戦略を立てるモデルにすれば, より現実に近いシミュレーションができるかもしれない.

 

戦略関数や高速化について

 別の記事に用意した。

 

moratoriamuo.hatenablog.com

【Acknowledgments】

人生初のAdvent Calendar参加でした. カレンダーの企画者のみゅーもり様, 共に参加している皆様, そしてここまで読んで下さった読者の皆様, ありがとうございました. ドンドンドーナツどーんと行こう!