聰明不如鈍筆
총명불여둔필
assignment KindeR

R로 깔끔하게 몬티 홀 문제 시뮬레이션 하기(feat. map_df())


미국 TV 게임 쇼 '렛츠 메이크 어 딜(Let's Make a Deal)'에서 유래한 '몬티 홀 문제'라는 퍼즐이 있습니다.


문 세 개 중에 하나를 선택해 문 뒤에 있는 선물을 가져갈 수 있는 퀴즈쇼에 참가한 당신.


맨 위에 있는 그림처럼 문 하나 뒤에만 자동차가 있고 나머지 문 뒤에는 염소가 있다. 사회자는 문 뒤에 무엇이 있는지 미리 알고 있다.


만약 이 참가자가 1번 문을 선택하면 사회자 몬티 홀 씨는 자동차가 없는 (이 그림에서는 2번) 문을 임의로 하나 열어 선택지를 줄여준다.


그리고 나서 1번 대신 3번을 선택하겠냐고 묻는다.


원래 선택했던 문 번호를 유지하는 것과 문 번호를 바꾸는 것 가운데 어느 쪽이 참가자가 자동차를 가져가는 데 유리할까?


사람들은 대부분 번호를 바꾸지 않는 쪽을 선택합니다.


사회자가 염소가 있는 문 하나를 알려줬기 때문에 자동차를 가져갈 확률이 $\frac{1}{3}$에서 $\frac{1}{2}$로 올라갔다고 생각하기 때문입니다.


실제로는 번호를 바꾸는 쪽이 유리합니다.


그냥 조금만 더 생각해 보면 이유를 짐작할 수 있습니다.


참가자가 문 세 개 가운데 하나를 골랐을 때 뒤에 자동차가 있을 확률은 $\frac{1}{3}$입니다.


이건 나머지 문 두 개 뒤에 자동차가 있을 확률이 $\frac{2}{3}$라는 뜻입니다.


그런데 사회자가 염소가 있는 문 하나를 골라줍니다.


그러면 나머지 문 하나가 이 $\frac{2}{3}$ 확률을 전부 가져가게 됩니다.



조금 더 쉽게 이해할 수 있도록 표를 하나 그려볼까요?


참가자가 전부 1번 문을 선택한 상황에서 현재 번호를 유지했을 때(스테이)와 바꿨을 때(스위치) 결과는 아래 표처럼 나타납니다.


 1번  2번  3번  사회자  스테이  스위치
 염소  염소  자동차  2번  염소(1번)  자동차(3번)
 염소  자동차  염소  3번  염소(1번)  자동차(2번)
 자동차  염소  염소  2·3번  자동차(1번)  염소(2·3번)


스테이했을 때는 자동차가 있는 문을 택했을 확률이 마지막 한 번 = $\frac{1}{3}$이지만 스위치 했을 때는 $\frac{2}{3}$로 올라갑니다.


한 번 R를 가지고 간단한 시뮬레이션을 진행해 정말 이런 결과가 나오는지 한번 알아보도록 하겠습니다.


언제나 그렇듯 제일 먼저 할 일은 tidyverse 패키지를 (설치하고) 불러오기.

#install.packages('tidyverse')
library('tidyverse')
## -- Attaching packages ---------------------------------------------------------------------------------- tidyverse 1.3.0 --
## √ ggplot2 3.3.0     √ purrr   0.3.4
## √ tibble  3.0.1     √ dplyr   0.8.5
## √ tidyr   1.1.0     √ stringr 1.4.0
## √ readr   1.3.1     √ forcats 0.5.0
## -- Conflicts ------------------------------------------------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()


이 문제를 해결하려면 임의로 번호를 고르는 과정을 반복해야 합니다.


R에서 이 기능을 담당하는 함수는 sample()입니다.


예를 들어 1~10 가운데 숫자 하나를 고르라는 코드는 아래처럼 쓸 수 있습니다.

sample(1:10, 1)
## [1] 10


또 벡터에서 특정 숫자만 제외하고 싶을 때는 대괄호([ ]) 뒤에 빼기(-) 부호를 써서 나타내면 됩니다.


예를 들어 1~10 가운데 1, 3, 5를 제외한 나머지 숫자를 나타내라고 할 때는 이렇게 쓰면 됩니다.

c(1:10)[-c(1, 3, 5)]
## [1]  2  4  6  7  8  9 10


노파심에 말씀드리면 여기서 c() 함수는 묶는다(concatenate)는 뜻입니다.


그리고 R에서 대괄호는 원래 몇 번째 원소를 출력하라는 뜻입니다.


만약 c(1:10)[c(1, 3, 5)]라고 입력했다면 그냥 1, 3, 5를 출력했을 겁니다.


이때 빼기를 넣으면 이 순서는 빼고 출력하라는 뜻입니다.


따라서 1~10 가운데 첫 번째, 세 번째, 다섯 번째 = 1, 3, 5를 제외하고 숫자 두 개를 임의로 추출하라는 내용은 이렇게 정리할 수 있습니다.

sample((1:10)[-c(1, 3, 5)], 2)
## [1] 2 7


이제 기본은 익히셨을 테니 1~3번 문 가운데 자동차가 있는 문 하나를 지정하겠습니다.

car <- sample(1:3, 1)
car
## [1] 3


이어서 참가자가 처음으로 선택하는 문 번호도 하나 임의로 고릅니다.

initial_pick <- sample(1:3, 1)
initial_pick
## [1] 2


이 경우에는 사회자가 고를 수 있는 문은 1번밖에 없습니다.


2번은 참가자가 선택했고 3번 뒤에는 자동차가 있기 때문.


단, 경우에 따라서는 참가자가 이미 자동차가 뒤에 있는 문을 골랐을 수도 있습니다.


이를 일반화해서 표현하자면 참가자가 골랐거나 자동차가 뒤에 있는 문을 빼고 하나를 사회자에게 고르라고 하면 됩니다.


이를 코드로 나타내면 아래처럼 쓸 수 있습니다.

monty_hall <- sample(c(1:3)[-c(initial_pick, car)], 1)
monty_hall
## [1] 1


그리고 참가자가 고르지도 않았고 사회자가 열지도 않은 마지막 문이 하나 있습니다.


이 문은 아래처럼 지정하면 됩니다.

the_other <- c(1:3)[-c(initial_pick, monty_hall)]
the_other
## [1] 3


우리는 스테이를 해야 할지 아니면 스위치를 해야 할지 궁금합니다. 


만약 현재 사례처럼 맨 마지막 문 뒤에 자동차가 있을 때는 문을 바꿔야 합니다.


이 내용을 코드로 표현할 때는 맨 마지막 문 번호와 자동차가 있는 문 번호가 같으면 스위치, 아닐 때는 스테이를 표현하도록 쓰면 됩니다.


이번에는 물론 스위치가 정답입니다.

result <- ifelse(the_other==car, 'switch', 'stay')
result
## [1] "switch"


여기 사용한 ifelse() 함수는 마이크로소프트(MS) 엑셀에서 쓰는 if() 함수와 같은 구실을 합니다.


그러니까 ifelse(조건, 조건이 맞을 때, 조건에 맞지 않을 때) 형태로 쓰시면 됩니다.


이제 지금까지 만든 내용을 switch_door라는 함수(function)로 만들어 보겠습니다.


프로그래밍 과정에서 함수는 반복적으로 사용한 코드 덩어리라고 생각하시면 쉽습니다.


R에서 함수는 function(...){} 형태로 정의합니다.


예를 들어 어떤 입력 값 x에 1을 더하는 함수 add()를 만들고 싶다면 아래처럼 쓸 수 있습니다.

add <- function(x){
  x + 1
}


그다음 add(1)이라고 입력하면 2를 출력합니다.

add(1)
## [1] 2


우리는 위에서 이미 틀을 다 짰기 때문에 이를 함수로 바꾸는 게 별로 어렵지 않습니다.

switch_door <- function(try){
  car <- sample(1:3, 1)
  initial_pick <- sample(1:3, 1)
  monty_hall <- sample(c(1:3)[-c(initial_pick, car)], 1)
  the_other <- c(1:3)[-c(initial_pick, monty_hall)]
  tibble(try=try,
         result=ifelse(the_other==car, 1, 0))
}


이 함수는 try라는 변수를 입력 받아 열을 하나 만드는 걸 제외하면 지금까지 우리가 논의한 걸 그냥 주욱 정리했을 뿐입니다.


위에서는 'switch'와 'stay' 형태로 시뮬레이션 결과를 정리했지만 이번에는 1과 0을 쓴다는 것도 차이점입니다.


아, tibble()은 tidyverse 생태계에서 사용하는 데이터 프레임, 그러니까 행과 열로 된 데이터 형태라고 생각하시면 됩니다.


이제 함수를 만들었기 때문에 위에 있는 코드를 전부 칠 필요 없이 이 함수만 호출하면 됩니다.

switch_door(1)
## # A tibble: 1 x 2
##     try result
##   <dbl>  <dbl>
## 1     1      1


result가 1이 나왔으니 이번에도 문 번호를 바꾸는 편이 유리합니다.


이제 이 작업을 여러 번 반복해 보도록 하겠습니다.


tidyverse() 생태계에서는 purrr 패키지에 들어 있는 map_*() 함수를 활용해 반복 작업을 처리합니다.


purrr 패키지 사용법이 낯선 분은 'R에서 깔끔하게 반복 작업 처리하기(feat. purrr)' 포스트가 도움이 될 수 있습니다.


map 함수는 기본적으로 map(반복할 대상, 반복할 함수) 형태로 씁니다.


예를 들어 숫자 1~5에 위에서 만든 add() 함수를 적용하고 싶으면 이렇게 쓰면 됩니다. 

map(1:5, add)
## [[1]]
## [1] 2
## 
## [[2]]
## [1] 3
## 
## [[3]]
## [1] 4
## 
## [[4]]
## [1] 5
## 
## [[5]]
## [1] 6


물론 리스트 형태뿐만 아니라 데이터 프레임 형태도 반복 작업을 할 수 있습니다.


먼저 add() 함수를 살짝 손보겠습니다.

add <- function(x){
  tibble(x=x, y=x+1)
}


이어서 map_df() 함수를 쓰면 데이터 프레임 형태로 결과가 나옵니다.

map_df(1:5, add)
## # A tibble: 5 x 2
##       x     y
##   <int> <dbl>
## 1     1     2
## 2     2     3
## 3     3     4
## 4     4     5
## 5     5     6


기본 원리를 이해하셨을 테니 이제 switch_door() 함수를 100번 반복해 보겠습니다.

map_df(1:100, switch_door)
## # A tibble: 113 x 2
##      try result
##    <int>  <dbl>
##  1     1      1
##  2     2      1
##  3     2      0
##  4     3      0
##  5     4      1
##  6     5      1
##  7     6      1
##  8     6      0
##  9     7      0
## 10     8      1
## # ... with 103 more rows


이상합니다. 분명 우리는 100번을 반복하고 싶었는데 이 113번 반복했습니다.


어디에 이상이 있는지 알아보는 차원에서 함수에 사용한 모든 변수를 출력해 보도록 하겠습니다.


그냥 tibble() 안에 변수를 전부 넣으면 됩니다.

switch_door <- function(try){
  car <- sample(1:3, 1)
  initial_pick <- sample(1:3, 1)
  monty_hall <- sample(c(1:3)[-c(initial_pick, car)], 1)
  the_other <- c(1:3)[-c(initial_pick, monty_hall)]
  tibble(try=try,
         car=car,
         initial_pick=initial_pick,
         monty_hall=monty_hall,
         the_other=the_other,
         result=ifelse(the_other==car, 1, 0))
}


그리고 위에서 쓴 것과 똑같은 방식으로 100번을 반복해 봅니다.

map_df(1:100, switch_door)
## # A tibble: 111 x 6
##      try   car initial_pick monty_hall the_other result
##    <int> <int>        <int>      <int>     <int>  <dbl>
##  1     1     2            3          1         2      1
##  2     2     2            2          1         3      0
##  3     3     2            3          1         2      1
##  4     4     3            1          2         3      1
##  5     5     3            3          2         1      0
##  6     6     2            3          1         2      1
##  7     7     3            2          1         3      1
##  8     8     3            3          1         2      0
##  9     9     2            1          3         2      1
## 10    10     1            3          1         2      0
## # ... with 101 more rows


열 번째 행을 보니 문제를 알 수 있습니다. 자동차가 1번 문 뒤에 있는데 사회자가 1번 문을 열었습니다.


자동차가 있는 문 번호 그리고 참가자가 고른 번호는 제외해야 하는데 sample() 함수 작동 과정에서 문제가 있었던 것.  


그래서 monty_hall 변수 부분을 손질합니다.


length() 함수로 monty_hall에 들어 있는 숫자가 몇 개인지 파악한 다음 숫자가 하나밖에 없을 때는 그냥 그 숫자를 출력하고 두 개일 때만 하나를 추출하도록 하는 방식입니다.


이 내용을 포함해 switch_door() 함수를 수정합니다.

switch_door <- function(try){
  car <- sample(1:3, 1)
  initial_pick <- sample(1:3, 1)
  monty_hall <- c(1:3)[-c(car, initial_pick)]
  monty_hall <- ifelse(length(monty_hall)==1, monty_hall, sample(c(monty_hall), 1))
  the_other <- c(1:3)[-c(initial_pick, monty_hall)]
  tibble(try=try,
         car=car,
         initial_pick=initial_pick,
         monty_hall=monty_hall,
         the_other=the_other,
         result=ifelse(the_other==car, 1, 0))
}


그다음 다시 반복 작업을 진행해 봅니다.

map_df(1:100, switch_door)
## # A tibble: 100 x 6
##      try   car initial_pick monty_hall the_other result
##    <int> <int>        <int>      <int>     <int>  <dbl>
##  1     1     1            3          2         1      1
##  2     2     1            1          2         3      0
##  3     3     3            2          1         3      1
##  4     4     3            1          2         3      1
##  5     5     1            3          2         1      1
##  6     6     3            1          2         3      1
##  7     7     2            2          3         1      0
##  8     8     3            3          1         2      0
##  9     9     3            1          2         3      1
## 10    10     2            2          1         3      0
## # ... with 90 more rows


이번에는 100개가 제대로 나왔습니다.


그래도 여전히 사회자가 문을 잘못 연 경우가 있을지 몰라 다시 확인해 봅니다.

map_df(1:100, switch_door) %>%
  filter(monty_hall!=car & monty_hall!=initial_pick)
## # A tibble: 100 x 6
##      try   car initial_pick monty_hall the_other result
##    <int> <int>        <int>      <int>     <int>  <dbl>
##  1     1     3            1          2         3      1
##  2     2     1            3          2         1      1
##  3     3     1            1          3         2      0
##  4     4     2            1          3         2      1
##  5     5     1            2          3         1      1
##  6     6     3            2          1         3      1
##  7     7     2            2          3         1      0
##  8     8     3            2          1         3      1
##  9     9     2            3          1         2      1
## 10    10     1            1          2         3      0
## # ... with 90 more rows


아닙니다. 다행스럽게도 사회자가 전부 규칙에 따라 문을 골랐습니다.


이 코드에는 파이프 기호(%>%)와 filter() 함수가 등장했습니다.


이는 전부 tidyverse 생태계에 들어 있는 dplyr 패키지를 통해 활용할 수 있는 기능입니다.


아래서 쓰게 될 summarise() 함수 역시 마찬가지.


아직 dplyr 활용법이 낯선 분이 계시다면 '최대한 친절하게 쓴 R로 데이터 뽑아내기(feat. dplyr)' 패키지가 도움이 될 수 있습니다.


100번을 반복한 다음에 result가 몇 번이나 나왔는지 평균을 내면 문을 바꾸는 게 옳은지 아닌지 알 수 있겠죠?


이를 알려주는 코드는 이렇게 쓸 수 있습니다.

map_df(1:100, switch_door) %>%
  summarise(win_pct=mean(result))
## # A tibble: 1 x 1
##   win_pct
##     <dbl>
## 1    0.67


이 포스트 처음에 이야기했던 것처럼 $\frac{2}{3}$에 가까운 값이 나왔습니다.


한 번 이 실험을 반복하면 어떤 추이를 나타내는지 알아볼까요?


이런 실험을 진행할 때는 누적 합을 계산하는 cumsum() 함수를 활용하면 좋습니다.

 

예를 들어 1부터 10까지 누적 합을 구하고 싶을 때는 아래처럼 쓰면 됩니다. 

crossing(try=1:10) %>%
  mutate(cumsum=cumsum(try))
## # A tibble: 10 x 2
##      try cumsum
##    <int>  <int>
##  1     1      1
##  2     2      3
##  3     3      6
##  4     4     10
##  5     5     15
##  6     6     21
##  7     7     28
##  8     8     36
##  9     9     45
## 10    10     55


우리는 result를 계속 더한 다음에 이를 반복 횟수로 나누면 문을 바꿨을 때 승률이 어떻게 변하는지 확인할 수 있을 겁니다.


이런 변화를 알아볼 때는 숫자만 나열하는 것보다는 그래프로 그리는 쪽이 더 보기 편할 겁니다.


이 내용을 코드로 정리하면 아래와 같습니다.

map_df(1:10000, switch_door) %>%
  mutate(sum=cumsum(result), win_pct=sum/try) %>%
  ggplot(aes(x=try, y=win_pct)) +
  geom_hline(yintercept=2/3) +
  geom_line() +
  scale_y_continuous(limits=c(0, 1))


실험을 반복할수록 가로로 그은 기준 선(=$\frac{2}{3}$)으로 수렴한다는 사실을 알 수 있습니다.


사람은 본능적으로 변화를 두려워하게 마련이지만 너무 변하지 않으려고 하다가 눈 앞에 있는 기회를 놓치는 건 아닌지 돌아볼 필요도 있겠습니다.


아, 혹시 위에 있는 코드에서 그래프를 그리는 부분이 이해가 가지 않는다는 분이 계시면 '최대한 친절하게 쓴 R로 그래프 그리기(feat. ggplot2)' 포스트가 도움이 될 수 있습니다.


그럼 모두들 Happy Tidyversing, Happy Charting!


댓글,

KindeR | 카테고리 다른 글 더 보기