1. 베이지안 신뢰구간 1

빈도주의 통계학에서 신뢰구간(Confidence Interval)이라고 얘기하는 개념이 베이지안에서는 신뢰구간(Credible Interval)이다. 적절한 용어가 통계학회에서도 제시하고 있지 않아 (혹은 제시하고 있는데 홍보가 덜되서 모르거나… 어쨌든…) 용어는 신뢰구간을 사용하지만, 실무에서 사용하는 의미는 동일하다.

2. 올해(2017년) 최고 수위타자 타율은?

전반기를 지나고 있는 현재시점(2017년 6월) 최고타자 수위타자 타율은 어떻게 될까? 이를 빈도주의 통계 및 베이즈 통계를 활용하여 95% 신뢰수준으로 구해본다.

2.1. 환경설정 및 데이터 준비

위키 KBO 수위타자 순위를 긁어온다. 80년대부터 최근까지 수위타자 타율의 변화를 시각화한다.

# 0. 환경설정 -------------------------------------------
# list_of_packages <- c("Lahman", "ggplot2", "tidyverse", "ggthemes", "extrafont",
#                       "tibble", "rvest", "stringr", "extrafont")
# new_packages <- list_of_packages[!(list_of_packages %in% installed.packages()[,"Package"])]
# if(length(new_packages)) install.packages(new_packages)
# 
# sapply(list_of_packages, require, character.only = TRUE)

# 1. 대한민국 최고 타자 데이터 ------------------------------------
# devtools::install_github('stefano-meschiari/latex2exp')
Sys.setlocale("LC_ALL", "English")
[1] "LC_COLLATE=English_United States.1252;LC_CTYPE=English_United States.1252;LC_MONETARY=English_United States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252"
hitter_url <- "https://ko.wikipedia.org/wiki/KBO_%EB%A6%AC%EA%B7%B8_%ED%83%80%EC%9C%A8_%EA%B4%80%EB%A0%A8_%EA%B8%B0%EB%A1%9D"
hitter_xml <- read_html(hitter_url)

hitter_tbl <- html_nodes(hitter_xml, xpath='//*[@id="mw-content-text"]/div/table[3]')
hitter_tbl <- html_table(hitter_tbl, fill=TRUE)[[1]]

Sys.setlocale("LC_ALL", "Korean")
[1] "LC_COLLATE=Korean_Korea.949;LC_CTYPE=Korean_Korea.949;LC_MONETARY=Korean_Korea.949;LC_NUMERIC=C;LC_TIME=Korean_Korea.949"
# 2. 최고 타자 타율 시각화 ------------------------------------

hitter_tbl %>% 
  mutate(시대 = cut(연도, breaks = c(0, 1990, 2000, 2010, 2020), labels=c("80년", "90년", "00년", "10년"))) %>% 
  ggplot(data = ., aes(x=타율, color=시대)) +
    geom_histogram(aes(y=..density..), binwidth=0.005, colour="gray", fill="darkgray") +
    geom_density(alpha=.1, aes(fill=시대)) +
    scale_x_continuous(limits=c(0.3, 0.45)) +
    theme_bw(base_family="NanumGothic") +
    theme(legend.position = "top") +
    scale_color_brewer(palette="Dark2") +
    labs(x="타율", y="밀도 (Density)") +
    geom_vline(aes(xintercept=mean(타율)),
             color="blue", linetype="dashed", size=1)

2.2. 최고 수위타자 사전분포

최고 수위타자에 대한 사전분포를 추정한다. 타석에 들어서서 안타를 치는 것은 이항분포가 되니, 이에 적합한 사전분포는 베타분포로 필요한 것이지만, 베타분포에 대한 \(\alpha, \beta\) 모수를 추정하면 된다. 경험적 베이즈 방법을 활용하여 사전분포 베타분포에 대한 모수를 추정한다. 적률법을 활용하여 추정한다. 데이터가 적은 경우 수렴이 안되는 상황이 발생할 수 있으니 그점을 유념한다.

# 3. 베타분포 ------------------------------------
## 3.1. 최고타자 사전분포 모수추정 ------------------------
hitter_avg_decade <- hitter_tbl %>% 
  mutate(시대 = cut(연도, breaks = c(0, 1990, 2000, 2010, 2020), labels=c("80년", "90년", "00년", "10년"))) %>% 
  dplyr::filter(시대 %in% c("10년"))

## 3.1. 적률법 ------------------------
alpha_beta <- MASS::fitdistr(hitter_avg_decade$"타율", dbeta,
                    start = list(shape1 = 560, shape2 = 970))

alpha0 <- alpha_beta$estimate["shape1"]
beta0 <- alpha_beta$estimate["shape2"]

## 3.2. 최고타자 사전분포 시각화 ------------------------
gg_x <- seq(0.3, 0.45, .001)
beta_dens <- dbeta(gg_x, alpha0, beta0) 

beta_df <- data.frame(gg_x, beta_dens) 

beta_df %>% 
  ggplot(data = ., aes(x=gg_x, y=beta_dens)) +
    geom_line() +
    scale_x_continuous(limits=c(0.3, 0.45)) +
    theme_bw(base_family="NanumGothic") +
    theme(legend.position = "top") +
    scale_color_brewer(palette="Dark2") +
    labs(x="타율", y="밀도 (Density)") +
    geom_vline(aes(xintercept=alpha0/(alpha0+beta0)),
               color="blue", linetype="dashed", size=1)

2.3. 현재 수위타자 경쟁

사전분포에 대한 선행작업이 완료되었다면, 다음으로 현재 수위타자 경쟁을 하고 있는 데이터를 긁어온다. 다음 스포츠 야구 타자 통계웹사이트에서 25명 데이터를 긁어와서 이중 수위타자가 될 가능성이 높은 10명만 추려본다.

## 3.4. 최고타자 타율 데이터 가져오기 ------------------------
baseball_url <- "http://score.sports.media.daum.net/record/baseball/kbo/brnk.daum"
baseball_html <- xml2::read_html(baseball_url)

Sys.setlocale("LC_ALL", "English")
[1] "LC_COLLATE=English_United States.1252;LC_CTYPE=English_United States.1252;LC_MONETARY=English_United States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252"
hitter_table <- rvest::html_nodes(x=baseball_html, xpath='//*[@id="table1"]')
baseball_hitter <- rvest::html_table(hitter_table)[[1]]
Sys.setlocale("LC_ALL", "Korean")
[1] "LC_COLLATE=Korean_Korea.949;LC_CTYPE=Korean_Korea.949;LC_MONETARY=Korean_Korea.949;LC_NUMERIC=C;LC_TIME=Korean_Korea.949"
## 3.5. 최고타자 타율 신뢰구간 ------------------------
baseball_bayes_df <- baseball_hitter %>% 
  mutate(선수 = str_extract_all(`선수 (팀)`, '^[^\\s]*'),
          팀 = str_extract_all(`선수 (팀)`,  "\\([^()]+\\)")) %>% 
  dplyr::select(선수, 팀, 타수, 안타, 타율) %>% 
  dplyr::filter(row_number() <=10 ) %>% 
  tbl_df %>% unnest()

2.4. 베이즈 타율 95% 신뢰구간

베이지안으로 추정한 타율은 사전분포에서 추정한 \(\alpha, \beta\)가 있다면 타수와 안타 데이터를 바탕으로 추정이 가능하다. 또한 사후분포의 \(\alpha, \beta\)를 갱신하면 95% 신뢰구간도 수월하게 계산할 수 있다. 빈도주의 기법과 비교하여 베이지안 기법의 차이점도 함께 비교해본다.

baseball_bayes_df <- baseball_bayes_df %>% 
  mutate(alpha1 = 안타 + alpha0,
         beta1 = 타수 - 안타 + beta0,
         베이즈_타율 =  round((안타 + alpha0) / (타수 + alpha1 + beta0), 3), 
         low  = qbeta(.025, alpha1, beta1),
         high = qbeta(.975, alpha1, beta1))

frequentist_df <- baseball_bayes_df %>%
  group_by(선수) %>%
  do(tidy(binom.test(.$안타, .$타수))) %>%
  select(선수, 타율=estimate, low = conf.low, high = conf.high) %>%
  mutate(method = "빈도주의 신뢰구간")

bayes_df <- baseball_bayes_df %>%
  select(선수, 타율=베이즈_타율, low, high) %>%
  mutate(method = "베이지안 신뢰구간")

comparison_df <- bind_rows(frequentist_df, bayes_df) %>% 
  ungroup()

DT::datatable(comparison_df)

2.4. 베이즈 타율 95% 신뢰구간 시각화

빈도주의 방법을 통한 95% 신뢰구간을 시각화하면 우선 신뢰구간이 상당히 넓은 점과 함께, 꿈의 타율 4할을 넘을 수 있는 선수가 다수 발견된다. 반면에 95%신뢰구간 상한으로 베이지안 추정을 통해 갈음해 보면 4할을 넘는 선수가 없다. 참고로 점선은 최고 수위타자 2010년대 평균 타율이다.

## 3.6. 최고타자 타율 신뢰구간 시각화 ------------------------

comparison_df %>% 
  mutate(선수 = reorder(선수, 타율)) %>%
  ggplot(aes(타율, 선수, color = method, group = method)) +
  geom_point(size=3.5) +
  geom_errorbarh(aes(xmin = low, xmax = high)) +
  geom_vline(xintercept = alpha0 / (alpha0 + beta0), color = "red", lty = 2) +
  labs(x="", y="", color = "") +
  theme(legend.position = "top")


  1. Understanding credible intervals (using baseball statistics)