0619-2잡소리MCMC

 1. 복권 사기?

 

 

얼마전에 복권 1등이 50명이나 나와서 이거 사기 아니냐라는 이야기가 있었는데요. 이걸 또 검증한다고 기사가 나오고. 참...

 

 

 

검증 기사 보니까 복권 평균 1억장 정도 팔린다고 합니다. 한 사람이 여러 장 사는 것이겠지요. 그리고 1등에 당첨될 확률이 약 900만분의 1이라고 합니다.

 

그럼 평균 m=n*p=약 12명 정도 됩니다. 이건 공식입니다.

 

개인이 복권 1장 샀을 때는 베르누이 시행이라고 합니다. 개인이 여러장 샀을 때는 이항분포라고 합니다. 이때 평균이 m=n*p가 됩니다.

 

그럼 1억 장이 팔렸다고 하면 정확한 분포는 1등 당첨 개수 Y=B(n, p)=B(1억, 1/900만) 이렇게 됩니다. 그럼 이 분포에서 이 당첨 개수 Y가 50보다 클 확률을 구해야 합니다. 그런데 이것 계산하기 쉽지 않죠.

 

 

좀 더 간단히 계산하는 방법이 없을까요. 만약 n이 1억 정도 된다고 합시다. 그럼 이 이항분포는 포아송 분포로 생각해도 됩니다. 즉 평균이 m인 포아송 분포가 됩니다. 그럼 이 포아송 분포에서 50보다 클 확률은 좀 더 간단히 계산해 낼 수 있습니다.

 

R 프로그램에서

 

1-ppois(50, 12)=0

 

 

즉 평균이 12인 포아송 분포에서 50명 이상의 1등 당첨자가 나올 확률은 0입니다. 그럼 이거 사기일까요?

 

아직은 모릅니다. 통계 프로그램에서 0이라고 나왔다고 해서 이게 0이 아닙니다. 소숫점 몇 자리에서 짤린 숫자입니다. 즉 0.0000000000000000000000039489 이럴 수 있다는 것입니다.

 

그럼 나올 확률이 엄청나게 낮아도 이게 처음 시도에 나올 수도 있고, 아니면 언제가는 한번 쯤 튀어 나올 수 있는거든요.

 

예를 들어 복권 1장 쌌을 때 1등 당첨 확률이 1/900만 이잖아요. 사실상 나올 확률이 0 이거든요. 그런데도 1등 당첨될 것 기대하고 사죠.

 

일반인 수준에서 더 간단하게 구할 수는 없을까요. 포아송 분포하면 뭔가 어려워 보이잖아요. 이건 일정 기간, 면적 안에서 사건이 일어날 횟수가 포아송 분포를 가집니다. 손흥민 선수가 한 게임에서 몇 골을 넣는가 하는 것도 포아송 분포를 합니다.

 

정규분포를 이용해도 됩니다. 즉 n이 엄청 크고 p가 매우 작으면 이건 평균이 m=n*p인 정규분포로 갑니다. 이 분포에서 50보다 클 확률을 구하면 됩니다.

 

 

그래서 사실상 확률이 0이지만 이렇다고 사기라고 확신을 할 수 없습니다.

 

그러나 50명 이상이 계속 나오거나 아니면 뛰엄뛰엄 나오더라도 많이 나오면 이건 무슨 문제가 있다고 판단해도 됩니다. 즉 과학의 문제가 아니라 사람들의 판단의 문제입니다.

 

조금 더 과학적으로 시도해보면 지금까지 뽑힌 1등 숫자와 정규분포와의 적합도를 비교해보면 됩니다. Q-Q plot을 그린다든지 아니면 K-S 검증을 한다든지 하면 된다는 것이죠.

 

 

문제가 있다고 해서 이게 꼭 내부의 조작질이라고 이야기할 수는 없습니다. 숫자를 뽑아내는 랜덤 메카니즘에 문제가 있을 수 있거든요.

 

 

 

2. MCMC

 

딴지 자유계시판에서 원주율 파이를 통계학으로 풀 수 있다는 계시물이 나와서. 저도 이게 뭔 소리인가 하다가 옛날에 책에서 본 기억이 살아나네요.

 

 

알파고가 나왔을 때 쓴 알고리즘이 MCMC라고 해서 한때 화제가 된 적이 있는데요. MCMC는 마코브 체인 몬테칼르로 알고리즘을 말합니다. 마코브체인 이론이 있고, 또 몬테칼르로 알고리즘이 따로 있습니다. 둘 다 분리해서 공부하면 그리 어려운 분야는 아닌데 이게 결합되는 것을 이해 하는 것이 힘듭니다.

 

 

원주율 파이 구하는 방법이 이 몬테칼르로 알고리즘의 대표적인 예입니다.

 

아래 그림처럼 길이가 1인 정사각형을 그려 놓고 거기에 내접하는 원을 생각합니다.

그런 다음 -1에서 1 사이의 랜덤넘버 짝을 구해 냅니다. 이렇게요

 

(x, y)=(-0.25, 0.54), (0.72, 0.31), (0.1, -0.34), 이렇게 (x, y) 짝 값을 랜덤하게 만 개 정도 뽑아냅니다. 그럼 이 점들이 원 안에 들어가는 경우가 있고, 아니면 원밖에 있는 경우가 있을 겁니다. 그럼 원 안에 있는 점의 개수의 비율을 구하면 원의 면적이 됩니다. 그럼 원의 면적을 알고 원의 원주율 파이도 구할 수 있죠. 점의 갯수 세는 것도 손으로는 안 되겠죠. 이것 프로그램도 간단합니다. 

 

 

이와 같이 확률적 랜덤 넘버를 구해 작업을 하는 것을 몬테칼르로 알고리즘이라 합니다.

 

 

또 다른 대표적인 예를 볼까요.

 

윤석렬이가 예능에 나와 적분 구하는 것을 이야기한 모양입니다. 저는 자세히 안봐서요.

 

아마 몬테칼르로 방법으로 적분하는 것을 이야기하려고 한 모양입니다. 제가 여기서 한번 쓴 적이 있는 윤석렬 아버지가 한국 통계학의 대부입니다. 아마 아버지 통해 예능에서 강의할 만한 것을 찾은 것 같습니다. 아버지가 이걸 아는게 아니라 주변에 자기 제자 교수들에게 물어 봤겠죠.

 

 

우리가 학교에서 적분 정의에 의해 구하는 법을 배웠죠. x축을 잘게 잘라서 얇은 사각형을 만들어 적분 면적을 구하는 것이죠.

 

 

그럼 이런 생각을 해 볼 수가 있습니다. 등간격으로 짜르지 말고 [0,1] 사이에 랜덤넘버를 구해 함수값에 집어 넣고 합을 구하면 되지 않을까 하는 생각을 할 수 있습니다.

 

네 이렇게 해도 됩니다. 여기서 합을 구하는 것이 아니라 평균을 구해야 합니다. 뽑은 랜덤 넘 숫자가 많으면 함수값을 더하면 엄청나게 큰 값이 나오거든요.

 

예를 들어 x제곱을 0에서 1까지 적분을 한다고 하죠. 그럼 계산하면 1/3이 나옵니다.

 

이걸 R에서 몬테칼르로 방법으로 하면 명령어 하나면 됩니다.

 

                                mean( ( runif(n, 0, 1))^2 )

 

==> 수정했습니다

 

이렇게 하면 됩니다.

 

n=10이면 0.207, n=100이면 0.3876, n=10000이면 0.3272, n=1000000이면 0.3330

 

이렇게 됩니다.

 

여기서 이런 것을 발견할 수 있습니다.

 

n이 커질수록 정답에 수렴한다. 현실세계에서는 어느 정도 수준에서 끊어야 한다. 그래서 정확한 값은 없습니다.

 

보면 정답 0.3333을 중심으로 변동이 심하게 일어납니다. 특히 초반에서요.

 

 

그래서 이런 문제가 생깁니다. n이 커질수록 정답에 수렴하는데 가능하면 빠르게 수렴하는 방법이 없을까. 그리고 정답 근처에서 변동이 심하게 일어나는데 이 변동을 줄일 수가 없을까 하는 문제입니다. 이런 문제를 variance reduction 문제라고 합니다.

 

 

이와 같이 하는 적분 방법을 몬테칼르로 방식에 의한 적분이라고 합니다.

 

 

그럼 여기서 조금 더 개념을 확장해 볼까요. 기댓값을 한번 구해 보죠

.

 

E[X제곱] = ∫x제곱*f(x)dx

                 =∫x제곱dP(조금 고상하게 표현하면, 여기서 P는 확률메저라고 합니 다)

 

 

입니다.

 

여기서 f(x)는 정규분포 등 확률분포입니다. 즉 x가 확률분포 f(x)에서 왔다는 것이죠. 그럼 이것도 적분 비슷합니다.

 

 

정규분포에서 랜덤 넘버 한 십만개 정도 뽑아낸 다음 이 십만개 값을 x제곱에 대입한 다음 평균값을 구하면 이 기댓값 E[X제곱]을 구할 수 있습니다.

 

 

그럼 이 정규분포에서 랜덤 넘버를 어떻게 구할 수 있을까요.

 

쉽지 않습니다. 정규분포는 아주 특이한 방법을 씁니다. 사실 균등분포 흔히 우리가 랜덤넘버라고 하는 것도 쉽지 않습니다. 그래서 유사, 가짜 랜덤넘버라고 합니다. psuedo 랜덤넘버라고 하죠. 진짜 랜덤넘버를 뽑을 수만 있다면 이건 세계적인 토픽감이고 뭔가 상을 하나 받을만 합니다.

 

딴지 게시판에 누가 댓글로 이걸 이야기 한 사람이 있는데 자연현상에서 랜덤 넘버를 뽑는 방법을 이야기한 것을 보았습니다. 사실 자기가 진짜 랜덤 넘버 뽑는 방법을 찾아다고 주장하는 학자들이 꽤 있습니다.

 

 

예를 들어 한 만 페이지 되는 한글 파일이 있다고 해보죠. 그럼 각 줄에 첫 철자의 자음만을 뽑아서 랜덤 넘버를 만들 수 있을 겁니다. ㄱ=1, ㄴ=2, 이렇게 코딩한 다음 첫째 줄 첫 문장이 강물이면 ㄱ이니까 1이 추출되고 이런식으로 랜덤 넘버를 뽑는다는 것이죠.

 

 

그런데 균등하게 뽑히지 않겠죠. 우리나라 단어에 있는 자음이 나오는 분포가 균등하다고 믿기에는 좀 그렇죠.

 

그래서 정리하면 이런 문제들이 있습니다.

 

확률분포에서 어떻게 랜덤하게 숫자를 뽑을 수 있을까?

 

또 가능하면 빠르게 수렴하는 방법이 없을까. 기상 예측 프로그램에서는 어마어마한 작업이 필요합니다.

 

그리고 가능하면 변동을 줄일 수가 없을까 하는 문제들이 있습니다.

 

 

여기서 rejection 방법, importance 방법 이런 것이 있는데 나중에 MCMC 공부하는데 중요한 역할을 합니다.

 

 

사실 제가 이걸 공부한 것은 지도교수가 independence 강의를 할 때입니다. 이런 강의는 교수가 아무 주제나 잡아 강의하는 것입니다. 관심 있는 대학원 생만 듣는 것이고요. 제 경우는 2-3명 정도 들었고, 이 주제도 교수가 선정한 몇 개 주제 중 하나였습니다. 볼쯔만 머신 강의하는 중에 간단히 나오는 내용입니다. 본인이 다 공부를 해야 했습니다. 즉 학기에 교과서 같은 것을 가지고 체계적으로 공부를 한 것이 아닙니다.

 

그럴 수 밖에 없는 것이 그 당시만 해도 이 분야가 엄청 햇갈리는 분야였습니다.

 

앞에서 이야기한 컴퓨터학과나 통계학과에서 배우는 수치해석 문제도 있고, 획기적인 논문이 metropolis 알고리즘인데(저널 논문 쓴 3사람 중의 한 사람입니다) 이거 온도 t, 에너지 E가 나오는 등 열역학 이야기이고, 그리고 MCMC의 특수한 케이스인 Gibbs sampler도 이미지 프로세싱에서 나온 이야기이고.

 

또 통계학은 MLE 한계점을 극복하려고 여기저기에 이상한 알고리즘들이 많이 나왔거든요.

 

그래서 그 당시만 해도 일관적이고 체계적인 설명이 없었습니다. 지금이야 수 많은 책들이 있지만요.

 

그래서 모듈 같은 것을 사가지고 바로 응용할 수도 있고, 이론을 이해 못해도 기계적으로 알고리즘을 이해해서 프로그램 해도 되지만 제대로 이해해서 공부하시면 상당한 시간이 걸릴겁니다.