【備忘録】シーズン全体の得点の分布をシミュレーションしてみた

こんにちは

サッカーデータ革命という本を読んでいたのですが、その中に試合ごとの得点数がシーズン全体でどのように分布するかは、ポアソン分布で説明できるという話がありました。
※書籍の画像は載せられませんが、似たようなことをされているサイトがありましたので、その画像を引用させていただきます。

こちらの画像のように、実際の得点数(Actual)の分布と、ポアソン分布(Poisson)で予測される分布が非常に近い値になるとのことです。

これを読んでバスケにも当てはめてみたくなったので、実際にシミュレーションしてみました。
結果としてはサッカーほど近い値にはならず、面白い結論にはならなかったのですが、備忘録やポートフォリオとして投稿しようと思います。

シミュレーションの流れ

まず最初に、どのようにシミュレーションするかを考える必要があります。
サッカーの場合は、ゴールが「試行回数が多く、成功確率が小さい」という特徴があり、こういった特徴はポアソン分布として捉えられるため、ポアソン分布でシミュレーションしています。
ただ、バスケの場合は、成功確率が小さいとは言えず、しかもシュートのタイプごとに得点と成功率が異なるため、サッカーとは別の形でシミュレーションする必要があります。
そこで今回は下記の流れでシミュレーションを行いました。

① 試合ごとに何回シュートを行うかを指数分布でランダムに算出
② n回目ののシュートがどのシュートタイプ(3Pシュート・2Pシュート・FT)になるかを多項分布でランダムに選択
③ ②で選択されたシュートが成功するどうかをベルヌーイ分布でランダムに選択し、得点に換算する
④ ②と③を①の回数だけ繰り返して、1試合での得点数を算出する
⑤ ④を想定される試合数だけ繰り返して試合ごとの得点数の分布をシミュレーションする

※指数分布：次に何かが起こるまでに経過する時間が従う分布。ある期間に平均して何回起こるかが分かれば分布が特定できる
※多項分布：複数の選択肢がある中で、各選択肢が何回起こったかを表す確率分布。
※ベルヌーイ分布：「成功か失敗か」のように２種類のみの結果しか得られないような実験、試行の結果を0と1で表した分布。

噛み砕いて説明すると、
・試合ごとに何回アテンプトするか(①, ④)
・試合内の各シュートは3P、2P、FTのどれになるか(②)
・打ったシュートは成功するか失敗するか(③)
それぞれを確率分布に則ってシミュレーションすれば、実際の得点の分布と比較できるという形です。

では実際にやってみましょう

結果

まずシミュレーションを回します。
(コードに興味のない方は飛ばしていただいて構いません)

# import libraries --------------------------------------------------------

library(tidyverse)
library(furrr)

# load data ---------------------------------------------------------------

source("../B_data.R")

data_base <- game_summary %>%
 inner_join(
   game %>%
     filter(Season == "2020-21", EventId == 2)
   , by = "ScheduleKey"
 ) %>%
 filter(OT1 == 0) #延長除外。予測値の分布考える際に邪魔になるので
 
f3ga <- sum(data_base$F3GA)
f2ga <- sum(data_base$F2GA)
fta <- sum(data_base$FTA)
f3g_per <- sum(data_base$F3GM) / sum(data_base$F3GA)
f2g_per <- sum(data_base$F2GM) / sum(data_base$F2GA)
ft_per <- sum(data_base$FTM) / sum(data_base$FTA)
game_team_cnt <- nrow(data_base)
atmpt_num <- (f3ga + f2ga + fta) / game_team_cnt #40分間の平均アテンプト回数
exp_num <- 250 #指数分布のnの値
pts_dist <- rep(NA, game_team_cnt) #予測された得点の分布

base_df <- data.frame(
 rowname = rep(seq(exp_num), each = 1),
 pattern = rep(c("f3ga", "f2ga", "fta"), each = exp_num)
) %>%
 mutate(
   atmpt_pt = case_when(
     pattern == "f3ga" ~ 3
     ,pattern == "f2ga" ~ 2
     ,pattern == "fta" ~ 1
   )
 )
flag_prob <- c(f3ga = f3ga, f2ga = f2ga, fta = fta) / sum(f3ga, f2ga, fta)
pts_prob <- c(f3ga = f3g_per, f2ga = f2g_per, fta = ft_per)

score <- function() {
 exp_vec <- rexp(n = exp_num, rate = atmpt_num) # 前のシュートとの時間の間隔を作成して格納
 atmpt_pred <- sum(cumsum(exp_vec) < 1)
 t <- as.tibble(
   rownames_to_column(
     as.data.frame(
       t(rmultinom(n = exp_num, size = 1, prob = c(f3ga, f2ga, fta)))
     )
   )
 ) %>%
   mutate(
     rowname = as.numeric(rowname)
   ) %>%
   rename(
     ,f3ga = V1
     ,f2ga = V2
     ,fta = V3
   ) %>%
   pivot_longer(
     -rowname
     ,names_to = "pattern"
     ,values_to = "flag"
   ) %>%
   filter(
     flag == 1
   )
 
 
 base_df %>%
   filter(rowname <= atmpt_pred) %>%
   inner_join(t, by = c("rowname", "pattern")) %>%
   #    filter(runif(n()) <= flag_prob[pattern]) %>% # flag相当
   mutate(pts = atmpt_pt * rbinom(n(), size = 1, prob = pts_prob[pattern])) %>% 
   pull(pts) %>%
   sum()
}

# 同試合数でシミュレーションver
set.seed(1031)
pts_dist <- map_dbl(seq(game_team_cnt), ~score())

pts_dist_tibble <- as_tibble(pts_dist) %>%
 mutate(
   type = "simulation"
 ) %>%
 rename(PTS = value)


#可視化
p_same_num <- data_base %>%
 select(PTS) %>%
 mutate(
   type = "actual"
 ) %>%
 bind_rows(pts_dist_tibble) %>%
 ggplot(aes(x = PTS, fill = type))+
 geom_histogram(position = "dodge", bins = 10)+
 labs(x = "得点数", y = "件数", title = "得点数の実際の分布とシミュレーションの比較(seed = 1031)")+
 theme_classic()
plot(p_same_num)

実際にシミュレーションした結果は下記になります。

両者とも山型の分布になっていますが、実際の分布のほうが中央に寄っていて、シミュレーションのほうが裾が広くなっているように見えます。サッカーのようにめちゃくちゃ合致している感じはありませんね。
(ちなみに得点の幅を10分割して、それぞれに当てはまる件数を算出してχ二乗検定を行ったところp値は0.009868と小さい値になりました)。

ただこれはシミュレーション1回分でたまたまこうなってしまったのかもしれません。例えば10回実施したら近づくかもです。
ということで、シミュレーションを10回分回してみました。

# 10倍の試合数でシミュレーションver
plan(multiprocess(workers = 6))
plan()
pts_dist <- future_map_dbl(seq(game_team_cnt * 10), ~ score(), .options = future_options(seed = 1031L)) 


pts_dist_tibble <- as_tibble(pts_dist) %>%
 mutate(
   type = "simulation"
 ) %>%
 rename(PTS = value)


#可視化
p_10times_num <- data_base %>%
 select(PTS) %>%
 mutate(
   type = "actual"
 ) %>%
 bind_rows(pts_dist_tibble) %>%
 ggplot(aes(x = PTS, fill = type))+
 geom_density(alpha=0.5)+
 labs(x = "得点数", y = "確率密度", title = "得点数の実際の分布とシミュレーションの比較(seed = 1031)", caption = "シミュレーションを10シーズン分行った")+
 theme_classic()
plot(p_10times_num)

結果はこちらです。

グラフの形式は変わっていますが、こちらでも変わらず実際の分布のほうが中央に寄っていて、シミュレーションのほうが裾が広くなっていますね。
そうなると、今回のシミュレーションでは捉えきれていない要素がありそうです。

ここからは仮説ですが、バスケのガベージタイムとクラッチタイムが加味しきれていないのが要因かもしれません。

ある試合でガベージタイムになったら主要メンバーを下げるため、その後の得点ペースがゆっくりになって極端な得点がなくなったり、クラッチタイムでも同じく試合のペースがゆっくりになって極端な得点がなくなったりするんじゃないかと思っています。こういった要素はサッカーは少なそうなので、サッカーの場合はポアソン分布に近似するのかなと考えています。(終盤まで得点を取るインセンティブが働くと考えられるので)

まあ、チームごとのスタッツを使わず、リーグ全体の平均アテンプト回数、アテンプト回数の内訳、それぞれの成功率のみを使ってある程度近い数値を出せるのは、それはそれで良い結果かもしれません。

ということで今回は以上です。
ご質問・ご指摘などありましたらコメントいただければと思います。
(敬称略)

Tweets by b__s__k__t

今回のデータやコードの格納場所

basketball_data_analysis/point_predict_distribution at master 揃 pauze-pauze/basketball_data_analysis