R을 이용한 하프 이닝 시뮬레이션 2

마르코프 체인 시뮬레이션

마르코프 체인 모델을 여러 번 시뮬레이션하여 2016년 야구의 반 이닝 동안 득점 분포를 얻을 수 있습니다. 첫 번째 단계는 모든 가능한 상태 간 전이에서 득점된 점수를 나타내는 행렬을 구성하는 것입니다. 여기서 $ N_{runners}$는 상태의 주자 수를, $O$는 아웃 수를 나타냅니다. 이닝에서 이미 타석에 들어선 모든 선수는 베이스에 있거나, 아웃되었거나, 득점했기 때문에 타격 플레이에서 득점 수는 다음과 같습니다.

$$runs = (N_{runners}^{(b)} + O^{(b)} + 1 ) - (N_{runners}^{(a)} + O^{(a)})$$

요컨대, 득점은 플레이 전 주자와 아웃의 합에서 플레이 후 주자와 아웃의 합을 뺀 값에 1을 더한 것입니다. 예를 들어, 1루와 2루에 주자가 있고 1 아웃인 상황에서 플레이 후 2루에 주자가 있고 2 아웃인 경우, 득점 수는 다음과 같습니다.

$$runs = (2 + 1 + 1) - (1 + 2) = 1$$

우리는 상태를 입력으로 받아 주자와 아웃의 수를 합산하여 반환하는 새로운 함수 num_havent_scored()를 정의합니다. 그런 다음 이 함수를 모든 가능한 상태에 적용(map_int() 함수를 사용)하고, 해당 합계를 벡터 runners_out에 저장합니다.

num_havent_scored <- function(s) {
  s |>
    str_split("") |>
    pluck(1) |>
    as.numeric() |>
    sum(na.rm = TRUE)
}

runners_out <- T_matrix |>
  row.names() |>
  set_names() |>
  map_int(num_havent_scored)

 

outer() 함수와 뺄셈 연산을 사용하여 모든 가능한 상태 쌍에 대해 득점 계산을 수행하고, 결과 행렬은 R_runs 행렬에 저장됩니다. R_runs 행렬을 살펴보면 몇 가지 음수 값과 이상하게 큰 양수 값이 있음을 알 수 있습니다. 하지만 이는 예를 들어 "000 0" 상태에서 "000 2" 상태로의 전이와 같은 전이가 불가능하기 때문에 문제가 되지 않습니다. 이 행렬을 정사각형으로 만들기 위해 cbind() 함수를 사용하여 이 득점 행렬에 0으로 이루어진 추가 열을 추가합니다.

R_runs <- outer(
  runners_out + 1, 
  runners_out, 
  FUN = "-"
) |>
  cbind("3" = rep(0, 24))

 

이제 우리는 새로운 함수인 simulate_half_inning()을 사용하여 야구의 반 이닝을 시뮬레이션할 준비가 되었습니다. 입력값은 확률 전이 행렬 P, 득점 행렬 R, 그리고 시작 상태 s(1부터 24 사이의 정수)입니다. 출력값은 해당 반 이닝에서 득점한 점수입니다.

simulate_half_inning <- function(P, R, start = 1) {
  s <- start
  path <- NULL
  runs <- 0
  while (s < 25) {
    s_new <- sample(1:25, size = 1, prob = P[s, ])
    path <- c(path, s_new)
    runs <- runs + R[s, s_new]
    s <- s_new
  }
  runs
}

 

이 시뮬레이션에는 두 가지 핵심 문장이 있습니다. 현재 상태가 s인 경우, sample() 함수는 전이 행렬 P의 s행을 사용하여 새로운 상태를 시뮬레이션합니다; 새로운 상태는 s_new로 표시됩니다. 이닝에서 득점한 총점수는 득점 행렬 R의 s행과 s_new열에 있는 값을 사용하여 업데이트됩니다.

 

map_int() 함수를 사용하면 많은 수의 야구 반 이닝을 시뮬레이션할 수 있습니다. 아래 코드에서는 주자 없고 아웃 카운트 없는 상태(상태 1)에서 시작하여 10,000개의 반 이닝을 시뮬레이션하고, 득점을 simulated_runs 벡터에 수집합니다. set.seed() 함수는 난수 시드를 설정하여 독자가 이 코드를 실행함으로써 이 특정 시뮬레이션의 결과를 재현할 수 있게 합니다.

set.seed(111653)
simulated_runs <- 1:10000 |>
  map_int(~simulate_half_inning(T_matrix, R_runs))

 

반 이닝에서 득점 가능한 점수를 찾기 위해, 우리는 table() 함수를 사용하여 simulated_runs의 값들을 표로 만듭니다.

table(simulated_runs)

 

우리의 10,000번 시뮬레이션에서, 5점 이상 득점한 반 이닝은 50 + 34 + 10 + 2 + 2 = 98회였습니다. 따라서 5점 이상 득점할 확률은 98 / 10,000 = 0.0098입니다. 이 계산은 sum() 함수를 사용하여 확인할 수 있습니다.

sum(simulated_runs >= 5) / 10000

 

simulated_runs에 mean() 함수를 적용하여 득점한 평균 점수를 계산합니다.

mean(simulated_runs)

 

10,000번의 반 이닝 동안, 평균 0.477점이 득점되었습니다.

 

다양한 주자와 아웃 상황에서의 득점 가능성을 이해하기 위해, 다른 시작 상태에 대해 이 시뮬레이션 절차를 반복할 수 있습니다. 우리는 상태 j에서 시작하여 득점한 평균 점수를 계산하는 runs_j() 함수를 작성합니다. map_int() 함수를 사용하여, 가능한 모든 시작 상태 1부터 24까지에 대해 runs_j() 함수를 적용합니다. 출력은 mean_run_value 열에 저장된 평균 득점의 벡터입니다. 이 값들은 아래에 시뮬레이션된 기대 득점 행렬로 표시됩니다.

runs_j <- function(j) {
  1:10000 |>
    map_int(~simulate_half_inning(T_matrix, R_runs, j)) |>
    mean()
}

erm_2016_mc <- tibble(
  state = row.names(T_matrix), 
  mean_run_value = map_dbl(1:24, runs_j)
) |>
  mutate(
    bases = str_sub(state, 1, 3),
    outs_ct = as.numeric(str_sub(state, 5, 5))
  ) |>
  select(-state)
  
erm_2016_mc |>
  pivot_wider(names_from = outs_ct, values_from = mean_run_value)