聰明不如鈍筆
총명불여둔필
한겨레 "서울시내 초등학교 신입생이 달랑 18명 '출산의 늪'"
Q. 교육통계서비스에 따르면 지난해(2019년) 기준으로 서울 지역 초등학교 학급별 학생 숫자는 평균 25명입니다. 그렇다면 한 반에 생일이 같은 학생이 존재할 확률은 얼마일까요?
A. 정답은 56.9%입니다.
예, 정말입니다. 얼핏 생각하면 (올해 같은 윤년을 제외하면) 1년은 365일이라 적어도 (365일 절반인) 183명은 있어야 생일이 같은 사람이 나올 것 같은 느낌이 드는 게 사실.
그런데 실제로는 23명만 모여도 생일이 같은 사람이 있을 확률(50.7%)이 절반을 넘어갑니다.
이렇게 '생일 문제' 그러니까 '사람이 임의로 모였을 때 그 중에 생일이 같은 두 명이 존재할 확률을 구하는 문제'는 직관과 어긋나기 때문에 '생일 역설(paradox)'이라고 부르기도 합니다.
The Pudding 'The Birthday Paradox Experiment'
이 블로그에서도 2015년 '이정희 그리고 생일 역설(문제)' 포스트를 통해 이 개념을 소개한 적이 있습니다.
2012년 제19대 국회의원 선거를 앞두고 당시 통합진보당 비례대표 경선 과정에서 온라인 투표자 가운데 이름은 다른데 주민번호 뒷자리가 똑같거나 엇비슷한 이들이 여럿 나왔습니다.
이에 대해 이정희 공동대표는 "우리나라 주민번호 체계상 동일한 지역에서 출생신고를 한 사람 20명만 모이면 그중 한 쌍 이상은 뒷번호 7개가 일치할 확률이 대단히 높다"고 말했습니다.
이 해명은 결국 거짓말로 밝혀졌지만 수학적으로 따지면 이 대표 이야기가 아주 허튼 소리는 아니라는 게 이 포스트 내용이었습니다.
수학적으로 '사람 n명이 임의로 모였을 때 모든 사람이 생일이 다를 확률' $\overline{p}(n)$은 아래처럼 계산할 수 있습니다.
$$\begin{align}\overline{p}(n) &= 1\times(1-\frac{1}{365})\times(1-\frac{2}{365})\times\cdots\times(1-\frac{n-1}{365})\\ &= \frac{365\times364\times363\times\cdots\times(365-n+1)}{365^n}\\ &= \frac{365!}{365^n(365-n)!} = \frac{n!\cdot\binom{365}{n}}{365^n} = \frac{\ _{365} \mathrm{P} _n}{365^n}\end{align}$$
따라서 1(=100%)에서 $\overline{p}(n)$을 빼면 '사람이 n명이 임의로 모였을 때 그 중에 생일이 같은 두 명이 존재할 확률' ${p}(n)$을 계산할 수 있습니다.
$$\begin{align}{p}(n) &= 1 - \frac{\ _{365} \mathrm{P} _n}{365^n}\end{align}$$
n=23으로 ${p}(n)$을 계산하면 0.5072972 ≒ 50.7%가 나옵니다. 같은 방식으로 n=25일 때 ${p}(n)$을 계산하면 0.5686997을 얻을 수 있습니다. 맨 처음에 56.9%가 나온 이유입니다.
이렇게 공식을 통해서만 생일 역설 문제에 접근할 수 있는 건 아닙니다. 한번 R를 가지고 생일 문제에 접근하는 모의실험(시뮬레이션)을 진행해 볼까요?
언제나 그렇듯 제일 먼저 불러 올 패지키는 tidyverse입니다.
library('tidyverse')## -- Attaching packages -------------------------------------------------------------------------------- tidyverse 1.2.1 --
## √ ggplot2 3.2.1 √ purrr 0.3.2 ## √ tibble 2.1.3 √ dplyr 0.8.3 ## √ tidyr 0.8.3 √ stringr 1.4.0 ## √ readr 1.3.1 √ forcats 0.3.0
## -- Conflicts ----------------------------------------------------------------------------------- tidyverse_conflicts() -- ## x dplyr::filter() masks stats::filter() ## x dplyr::lag() masks stats::lag()
우리가 처음 쓸 함수는 sample()입니다. sample은 이름 그대로 모집단에서 표본을 추출하는 구실을 합니다.
예를 들어 1부터 10까지 숫자 가운데 다섯 개를 고르되 복원을 허용한다고 명령하고 싶을 때는 아래처럼 쓰면 됩니다.
sample(10, 5, replace=TRUE)## [1] 7 1 10 2 6
우리는 1년 365일 가운데 한 반 학생 25명 생일을 고르고 싶으니까 이렇게 쓰면 되겠죠?
sample(365, 25, replace=TRUE)## [1] 5 216 294 269 336 4 211 87 13 326 11 277 49 73 19 338 241 320 322 ## [20] 198 116 332 350 340 19
위에 있는 결과에는 19가 두 번 나오는데 얼핏 봐서는 겹치는 자료가 있는지 알기가 쉽지 않습니다.
그럴 때는 duplicated() 함수 도움을 받으면 됩니다.
아래 코드에서는 1이 두 번 나오고 duplicated() 함수는 두 번째가 중복이 맞다(TRUE)고 알려주고 있습니다.
c(1, 1, 2, 3, 4) %>%
duplicated()## [1] FALSE TRUE FALSE FALSE FALSE
노파심에 말씀드리면 여기 등장하는 '%>%'는 파이프 기호입니다.
현실 세계에서 파이프가 한 곳에서 다른 곳으로 액체나 기체를 보내주는 것처럼 코딩에서 파이프도 함수 사이에서 자료를 전달하는 구실을 맡습니다.
참고로 위에 있는 코드는 아래와 100% 똑같은 뜻입니다.
duplicated(c(1, 1, 2, 3, 4))파이프를 써서 코드를 쓰는 방법에 익숙하지 않으시다면 '최대한 친절하게 쓴 R로 데이터 뽑아내기(feat. dplyr)' 포스트가 도움이 될 수 있습니다.
다시 본론으로 돌아가겠습니다. 우리는 그저 생일이 겹치는 사람이 있는지 아닌지만 알고 싶을 뿐이니까 몇 번째가 겹치든 관계가 없습니다.
어디에든 겹치는 원소가 있는지 아닌지만 알고 싶은 거죠. 이럴 때는 any() 함수를 쓰시면 그만입니다.
c(1, 1, 2, 3, 4) %>%
duplicated() %>%
any()## [1] TRUE
이 코드를 조금만 바꾸면 25명이 모여 있을 때 생일이 같은 사람이 있는지 없는지 알 수 있습니다.
sample(365, 25, replace=TRUE) %>%
duplicated() %>%
any()## [1] TRUE
이번에는 '우연히' 생일이 같은 사람이 있었습니다. 거꾸로 없을 때도 있습니다.
어떤 비율로 이런 사람이 있는지 없는지 알아보려면 위에 쓴 코드를 아주 많이 반복해 보면 됩니다. 그게 바로 시뮬레이션입니다.
R에서 '깔끔하게' 반복 작업을 하고 싶을 때는 map_*() 함수를 쓰면 됩니다.
아래 코드는 위에서 한 번만 했던 작업을 100번 반복하라는 뜻입니다.
map_dbl(1:100, ~sample(365, 25, replace=TRUE) %>%
duplicated() %>%
any())## [1] 1 0 1 1 1 0 1 0 1 0 1 1 1 1 0 0 0 1 1 0 1 1 0 1 0 0 1 0 0 1 0 0 0 1 1 1 1 ## [38] 1 0 0 1 1 1 1 0 0 0 0 0 0 1 1 0 1 0 1 1 0 0 1 1 0 1 1 1 1 1 0 1 1 1 1 1 1 ## [75] 0 0 1 0 1 1 1 1 0 1 1 0 1 1 0 0 0 1 1 0 1 1 0 1 0 0
이 과정이 잘 이해가가 가지 않으시다면 'R에서 깔끔하게 반복 작업 처리하기(feat. purrr)' 포스트를 읽어 보셔도 좋습니다.
이번에도 숫자가 따로 따로 나와서 어느 정도 비율이 되는지 짐작하기가 쉽지 않습니다.
이럴 때는 sum() 함수로 그냥 이 값을 더해주면 그만입니다.
map_dbl(1:100, ~sample(365, 25, replace=TRUE) %>%
duplicated() %>%
any()) %>%
sum()## [1] 63
100번 실행을 진행해 63번 생일이 겹치는 사람이 나왔으니까 이때 확률은 63%입니다.
우리는 정답(56.9%)을 알고 있는데 이 숫자는 너무 큽니다.
10만 번 시뮬레이션을 돌리면 어떤 결과가 나올까요?
map_dbl(1:100000, ~sample(365, 25, replace=TRUE) %>%
duplicated() %>%
any()) %>%
sum()## [1] 56783
10만 번 가운데 5만6783번이면 56.8%입니다. 이 정도면 꽤 믿을 만한 결과라고 할 수 있지 않을까요?
사람이 25명일 때 문제를 해결하는 방법을 찾았으니 이제 사람 숫자별로 같은 계산을 반복하면 그만입니다.
우리는 crossing() 함수를 써서 계산 결과를 확장할 겁니다.
'R로 스도쿠 퍼즐을 풀어보자(feat. function(function(x)))' 포스트에서도 소개해 드렸던 crossing() 함수는 R 기본 함수인 expand.grid()처럼 결합 가능한 모든 조합을 만드는 구실을 합니다.
예를 들어 숫자 1~5와 알파벳 a~e를 서로 짝짓게 만들려면 아래처럼 쓰면 됩니다.
crossing(a=1:5, b=letters[1:5])## # A tibble: 25 x 2 ## a b ## <int> <chr> ## 1 1 a ## 2 1 b ## 3 1 c ## 4 1 d ## 5 1 e ## 6 2 a ## 7 2 b ## 8 2 c ## 9 2 d ## 10 2 e ## # ... with 15 more rows
일단 우리는 사람이 2명~366명일 때 이 계산을 100번씩 반복해 볼 겁니다.
이걸 R에게 명령하려면 아래처럼 쓰면 됩니다
crossing(saram=2:366, try=1:100) %>%
mutate(birthday=map(saram, ~sample(365, .x, replace=TRUE)))## # A tibble: 36,500 x 3 ## saram try birthday ## <int> <int> <list> ## 1 2 1 <int [2]> ## 2 2 2 <int [2]> ## 3 2 3 <int [2]> ## 4 2 4 <int [2]> ## 5 2 5 <int [2]> ## 6 2 6 <int [2]> ## 7 2 7 <int [2]> ## 8 2 8 <int [2]> ## 9 2 9 <int [2]> ## 10 2 10 <int [2]> ## # ... with 36,490 more rows
여기 등장하는 mutate() 함수는 계산 결과를 바탕으로 열을 추가하라는 뜻입니다.
mutate() 함수를 어떻게 쓰는지 모르시겠다면 '최대한 친절하게 쓴 R로 데이터 뽑아내기(feat. dplyr)' 포스트를 읽고 오셔도 좋습니다.
마찬가지로 map()을 쓴 이유가 궁금하시다면 'R에서 깔끔하게 반복 작업 처리하기(feat. purrr)' 포스트가 도움이 될 수 있습니다.
자, 다시, 위에서 보이는 birthday 열에 정수(integer)가 두 개씩 있는 건 사람이 두 명이라 날짜를 두 개씩 추출했기 때문입니다.
이렇게 데이터가 중첩(nested) 상태일 때는 unnest() 함수로 내용을 확인할 수 있습니다.
crossing(saram=2:366, try=1:100) %>%
mutate(birthday=map(saram, ~sample(365, .x, replace=TRUE))) %>%
unnest(birthday)## # A tibble: 6,716,000 x 3 ## saram try birthday ## <int> <int> <int> ## 1 2 1 113 ## 2 2 1 68 ## 3 2 2 294 ## 4 2 2 121 ## 5 2 3 137 ## 6 2 3 356 ## 7 2 4 353 ## 8 2 4 308 ## 9 2 5 296 ## 10 2 5 101 ## # ... with 6,715,990 more rows
그러니까 사람이 두 명일 때 첫 번째 실험은 113과 68을 뽑았습니다. 그리고 두 번째로 294와 121을 뽑았습니다.
현재 이런 식으로 각 사람 숫자에 따라 100번씩 실험을 진행한 상태입니다.
25명만 가지고 했을 때처럼 생일(birthday)에 겹치는 게 있는지 없는지 알려주면 더 좋겠죠?
열을 하나 더 추가하겠습니다.
crossing(saram=2:366, try=1:100) %>%
mutate(birthday=map(saram, ~sample(365, .x, replace=TRUE)),
duplicated=map_lgl(birthday, ~.x %>% duplicated() %>% any()))## # A tibble: 36,500 x 4 ## saram try birthday duplicated ## <int> <int> <list> <lgl> ## 1 2 1 <int [2]> FALSE ## 2 2 2 <int [2]> FALSE ## 3 2 3 <int [2]> FALSE ## 4 2 4 <int [2]> FALSE ## 5 2 5 <int [2]> FALSE ## 6 2 6 <int [2]> FALSE ## 7 2 7 <int [2]> FALSE ## 8 2 8 <int [2]> FALSE ## 9 2 9 <int [2]> FALSE ## 10 2 10 <int [2]> FALSE ## # ... with 36,490 more rows
이제 사람 숫자 별로 그룹을 지어서 몇 개나 겹치는지 알아보면 그만입니다.
아까는 합계를 구했는데 TRUE는 1, FALSE는 0값을 가지고 있기 때문에 평균(mean)을 내도 결과를 확인할 수 있습니다.
crossing(saram=2:366, try=1:100) %>%
mutate(birthday=map(saram, ~sample(365, .x, replace=TRUE)),
duplicated=map_lgl(birthday, ~.x %>% duplicated() %>% any())) %>%
group_by(saram) %>%
summarise(probability=mean(duplicated))## # A tibble: 365 x 2 ## saram probability ## <int> <dbl> ## 1 2 0 ## 2 3 0.01 ## 3 4 0.01 ## 4 5 0.03 ## 5 6 0.02 ## 6 7 0.06 ## 7 8 0.07 ## 8 9 0.08 ## 9 10 0.12 ## 10 11 0.14 ## # ... with 355 more rows
group_by(), summarise() 콤비를 이해하기 어렵다고 생각하시는 분도 '최대한 친절하게 쓴 R로 데이터 뽑아내기(feat. dplyr)' 포스트가 도움이 될 수 있습니다.
계속 앞쪽만 보고 있으니 결과가 잘 나왔는지 아닌지 헷갈리기도 합니다. tail() 함수로 뒤쪽을 한 번 볼까요?
crossing(saram=2:366, try=1:100) %>%
mutate(birthday=map(saram, ~sample(365, .x, replace=TRUE)),
duplicated=map_lgl(birthday, ~.x %>% duplicated() %>% any())) %>%
group_by(saram) %>%
summarise(probability=mean(duplicated)) %>%
tail()## # A tibble: 6 x 2 ## saram probability ## <int> <dbl> ## 1 361 1 ## 2 362 1 ## 3 363 1 ## 4 364 1 ## 5 365 1 ## 6 366 1
그래프를 한 번 그려보겠습니다. 뒤 쪽은 어차피 1이 나올 테니 75명까지만 한번 그려볼까요?
crossing(saram=2:75, try=1:100) %>%
mutate(birthday=map(saram, ~sample(365, .x, replace=TRUE)),
duplicated=map_lgl(birthday, ~.x %>% duplicated() %>% any())) %>%
group_by(saram) %>%
summarise(probability=mean(duplicated)) %>%
ggplot(aes(x=saram, y=probability)) +
geom_line()이번에는 ggplot()라는 녀석이 등장했습니다. 이런 그래프를 그리는 법을 잘 모르시겠다면 '최대한 친절하게 쓴 R로 그래프 그리기(feat. ggplot2)' 포스트를 읽어 보시면 도움이 될 겁니다.
직접 이 코드를 실행해 보신 분은 아시겠지만 이 정도만 그리는 데도 생각보다 시간이 오래 걸립니다.
68명만 있어도 이미 생일이 같은 사람이 있을 확률이 99.9%입니다.
2명부터 68명까지만 1만 번씩 시뮬레이션을 진행하고 나서 그래프를 그려보겠습니다.
crossing(saram=2:68, try=1:10000) %>%
mutate(birthday=map(saram, ~sample(365, .x, replace=TRUE)),
duplicated=map_lgl(birthday, ~.x %>% duplicated() %>% any())) %>%
group_by(saram) %>%
summarise(probability=mean(duplicated)) %>%
ggplot(aes(x=saram, y=probability)) +
geom_line()확실히 시뮬레이션 횟수가 늘어나니 그래프가 훨씬 매끈합니다.
여기서 치사한(?) 말씀을 하나 드리겠습니다.
우리는 이렇게 열심히 시뮬레이션을 진행하고 있지만 사실 R은 사람 숫자에 따라 생일이 같은 사람이 있을 확률을 알려주는 pbirthday() 함수를 마련해 두고 있습니다.
위에서 사람이 68명만 있어도 99.9%라고 말씀드린 것 역시 이 함수를 가지고 계산한 결과입니다.
pbirthday(68)## [1] 0.9987264
그렇다면 이 시뮬레이션 결과가 실제 계산값과 얼마나 차이가 나는지도 알아볼 수 있겠죠?
아래는 공식으로 계산한 값 위에 시뮬레이션 결과를 그리라는 코드입니다.
crossing(saram=2:68, try=1:10000) %>%
mutate(birthday=map(saram, ~sample(365, .x, replace=TRUE)),
duplicated=map_lgl(birthday, ~.x %>% duplicated() %>% any())) %>%
group_by(saram) %>%
summarise(probability=mean(duplicated)) %>%
mutate(actual=map_dbl(saram, pbirthday)) %>%
ggplot(aes(x=saram, y=probability)) +
geom_line(aes(y=actual), color='#00bfa5', lwd=1.25) +
geom_line()이렇게 dplyr에 들어 있는 crossing()을 바탕으로 tidyverse 생태계에 들어 있는 각종 패키지를 조립하면 시뮬레이션도 어렵지 않게 진행할 수 있습니다.
그러면 여러분 모두 Happy Tidyversing!
댓글,
