Alias method Algorithm

"Alias method" 를 사용하여 probability ( 확률 ) 기반의 RandomValue ( Index ) 를 얻고자 하는데 관련 한글 문서가 없었다. 영문 위키를 찾아봤지만 생략된 내용이 많아 이해하기 어려웠다. 구글링을 하다 좋은 영문 포스팅을 찾게 되었고

포스팅 : https://www.keithschwarz.com/darts-dice-coins/
위키 : https://en.wikipedia.org/wiki/Alias_method

 

이를 통해 대략적인 원리를 이해하는데 도움을 받았다.

그래서 누군가에게 도움이 될까 싶어 이번 기회에 한글로 정리하고자 한다. 원문에서는 수식 및 공식이 많고 이미지가 많은데 최대한 정리하여 핵심 아이디어만 작성한다.

---

Alias method 를 설명하기 전에 확률 관련 설명부터 시작하는데 원문에서는 RandomValue ( 임의의 값 ) 을 얻는 확률의 이야기를 두 가지 예시로 시작하고 있다.

"동전 던지기 ( 앞 / 뒤 )" 와 "n 면을 가진 주사위를 던지기" 이다.
동전 던지기는 각 면이 1 / 2 확률을 가지고 있고,
주사위 던지기는 1 / n 의 확률을 가지게 된다.

만약, 우리가 동전 던지기를 하고 앞면 또는 뒷면 을 얻는 프로그램 ( 알고리즘 ) 을 작성한다면 어떻게 해야할까?

앞으로 설명되는 모든 실험 ( 코드 ) 에는 전제조건이 있다.

- [0, 1) [닫힌, 열린) 의 임의의 유리수를 얻는 함수는 비용이 O(1) 이다
- 위 함수는 균일 분포 ( uniformly-distributed ) 를 잘 따른다

그 외 다른 following operations 도 있지만 원리를 이해하는데 큰 의미는 없어서 생략한다.

 

예를 들어

Random.NextSingle() // 0 ~ 1 의 값을 얻는 dot. code

과 같은 함수로 [0, 1) 사이의 값을 얻고 0.5 보다 작으면 앞면, 크면 뒷면으로 구분하면 된다.
주사위도 같은 개념으로 RandomValue 를 얻고 n 을 곱하여 그 수를 floor ( 내림 ) 을 하면 원하는 면 ( Index ) 가 나온다.

 

N 이 4 인 경우 각 n ( 면 ) 의 확률 범위는 아래와 같다.
Index : [0, 1, 2, 3]
Range : [ 0 / 0.25 / 0.5 / 0.75 / 1 ]

예시) 0.1 인 경우는 0. 1 * 4 -> 0.4 ->  0 이 나온다
예시) 0.6 인 경우는 0.6 * 4 -> 2.4 -> 2 이 나온다

 

이것이 원문의 "Simulating a Fair Die" 이다.
공평한 확률을 가지는 Die ( 주사위 ) 이다.

---

다음은 "Simulating a Loaded Die with a Fair Die" 인데 주사위의 경우 동일한 숫자의 면을 여러 개 가지는 경우도 있다.
이러면 각 면이 가지는 확률이 동일하지 않는데 이런 경우의 RandomValue 는 어떻게 얻을 수 있을까?

예시는 아래와 같다.

Index : [0, 1, 2, 3]

Probability : [1/2, 1/3, 1/12, 1/12]
Index [0, 1, 2, 3] 의 모든 확률을 더하면 1 이 나온다.

결과부터 이야기하면 최소공배수 ( LCM ) 을 사용하면 된다.
* LCM ( Least Common Multiple )

 

쉽게 이야기하면 예시의 분모인 ( 2, 3, 12 ) 의 최소공배수인 "12" 로 분모를 모두 치환하면 확률을 [6/12, 4/12, 1/12, 1/12] 로 바꿀 수 있다.
이 상태에서 "Simulating a Fair Die" 와 같이 N 이 12 개인 주사위로 생각하고 RandomValue 를 얻은 뒤 n 을 곱하여 floor 를 하면 원하는 Index 가 나온다.
이 Index 가 각 면에 해당되진 않기에 적절한 MappingTable ( Array ) 를 경유해 원하는 Index 를 찾을 수 있다.
위 예시에서 MappingTable 은 아래와 같다.

총 길이가 12 개의 배열
MappingTable : [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 2, 3]
* Probability 의 분자 값 만큼 MappingTable 의 Item 수가 많아진다

원문에서는 Algorithm 마다 성능 분석 ( big O for Time / Memory ) 을 테이블로 작성했는데 해당 내용이 궁금하면 원문을 찾아가서 볼 수 있다.

다시 돌아와서 "Simulating a Fair Die" 는 계산이 빠르다 ( 계산 자체는 Array 를 Indexing 하면 끝 )
다만, MappingTable 의 크기가 LCM 의 값이기에 주어진 상황에 따라 Array 크기가 너무 커서 실제 사용이 불가능한 경우가 있다.

여기서 다른 개념이 필요해서 잠깐 주제가 변경되는데 "Simulating a Biased Coin" 이야기가 나온다.

"Biased Coin" 은 동전 던지기를 말하는데 이 때 동전의 앞면 / 뒷면이 나올 확률이 다른 경우이다.
* 위에서 1/2 확률로 동전 던지기의 결과를 구하는 방법을 작성했는데 이는 사실 N 이 2 인 Die 라고 생각하면 된다.

 

Binary 에서 구분이 되는 확률 P 가 있다면 [0, 1) 의 결과가 P 보다 작은지 큰지만 확인하면 앞면 / 뒷면을 알 수 있다는 이야기다.

---

"Simulating a Fair Die with Biased Coins" 는 "Simulating a Fair Die" 를 위 동전으로 계산하는 방법을 말한다.

이전 Die 의 경우 LCM 을 사용했지만 Biased Coin 을 사용하면 LCM 을 사용하지 않을 수 있다.
예를 들어 1/4 확률을 가지는 N 개의 면을 가진 주사위가 있으면 첫번째 면이 나올 확률은 1/4 이다.
그리고 이 확률 p 를 동전의 경우로 치환해서 [0, 1) 의 RandomValue 가 1/4 보다 작으면 Index 값을 "0" 으로 판단한다.
만약 RandomValue 가 1/4 보다 크다면 ( = 3/4 ) Index "0" 은 실패한 것이니 다시 동전을 던져야 한다.
두 번째 던진 동전이 heading ( 성공 ) 했는지 판단하는 확률 p 는 이전과 같이 1/4 가 아니다.
이전에 Index "0" 의 검사가 끝났으므로 두 번째는 나머지 3 개의 면 중에 하나가 나올 확률인 "1/3" 가 되어야 한다.

이런 과정으로 for ( i -> n-1 ) 에서 각 i 가 나올 확률은 1/(n-i) 으로 매 번 동전 던지기를 하면 된다.
* 마지막 면은 항상 p 가 1 이다

이 알고리즘은 소요 시간이 최악의 경우 N 번 수행되어야 하므로 효율이 좋지 않다.

---

동일 확률을 가지는 N 개의 주사위를 봤으니 이번에는 "Simulating a Loaded Die with a Biased Coin" 이다.
각 면이 다른 확률을 가지는 주사위의 예시를 들어보자.

예시는 똑같이 Probability : [1/2, 1/3, 1/12, 1/12] 으로 시작한다.
기본 아이디어는 "Simulating a Fair Die with Biased Coins" 와 비슷한데 각 면이 나올 확률만 다시 계산하면 된다.
* 동전 던지기가 성공하는 확률 p 를 구하는 방법만 다르다

첫 번째 "1/2" 의 경우 당연히 1/2 확률로 RandomValue 가 나오면 성공이고 실패하는 경우는 1 - 1/2 인 나머지 공간 ( 1/2 ) 에서 1/3 의 확률이 나와야 하는데 이런 방식으로 모든 면에서 발생하는 확률을 더하면 전체 확률이 1 이 나와야 하기 때문에 이 때 "1/3" 을 바로 쓰는 것이 아닌 이전 결과가 누적된 확률을 사용한다.

원문에서는 누적되 값을 "mass" 라고 표현했는데 그 내용은 이전 확률의 누적 공간을 뺀 공간에서 p 를 계산한다.
첫 시작은 mass = 1 값부터 시작하고 매 iterator 가 돌면서 p / mass 의 결과인 "p1" 을 동전 던지기 확률로 사용하고 실패한다면 mass 에서 p 를 빼주고 다시 동전을 던진다.

위 방법은 역시나 최악의 경우에 N 번 던지기를 해야하므로 효율이 좋지 않다.

---

다음 "Generalizing Biased Coins: Simulating Loaded Dice" 는 비슷하지만 다른 아이디어를 사용했다.
매번 확률을 구하는 것이 아니라 각 면의 확률을 범위로 먼저 계산하는 아이디어다.

 

예를 들어 probability : [1/4, 1/5, 1/8, 1/8, 1/10, 1/10, 1/10] 인 경우 각 p 의 합을 1로 만들되 LCM 으로 치환하여 누적된 확률 범위 [10/40, 18/40, 23/40, 28/40, 32/40, 36/40, 40/40] 으로 바꾸는 것이다.
* 원문에서는 LCM 을 쓰진 않고 누적만 했다

 

그 다음 RandomValue 를 얻고 위 영역에 도달되는 범위를 Index로 치환하면 된다.

이 경우, RangeArray 를 만드는 과정이 필요한데 간단하게 N 크기의 Array  를 만들고 첫 번째 값은 "0" 다음 값은 이전 값에 p 를 더하면 된다. 이후 RandomValue 값을 Array 내 어디 영역인지 검사해야 하는데 BinarySearch 를 사용하면 된다.
* Array 생성 시 미리 정렬을 했다는 가정이 있다
BinarySearch 는 비용이 logN 이므로 이전 알고리즘의 worst cost (n) 보다는 나은 방법이다.

요 방법을 "Roulette Wheel Selection" 와 유사해서 비슷한 이름으로 부른다고 한다

 

여기서 조금 개선된 아이디어도 있는데 BinarySearch 대신 BinaryTree 를 사용한다.
* 그냥 BinaryTree 가 아닌 well-balanced BT 를 사용한다. 그러면 best (1) 을 얻을 수 있다
* BT 사용 시, Tree 생성 비용인 n^2 가 발생한다

---

드디어.. Alias method 의 아이디어가 되는 "Throwing Darts" 가 설명된다.
위키와 같은 영문 글에서는 다트를 던진다고 하는데 솔직히 무슨 말인지 이해가 안됐다.
이제 천천히 알아보자.

다트를 던지는 방식과 이전 방법의 차이점은 "[0, 1)" 라는 1 차원에서 확률을 표현하기 위해 각 p 를 계산했다면 Darts 는 2 차원 "면" 에서 다트를 던져 맞는 공간을 원하는 결과 ( Index ) 로 사용한다는 점이다.
* 그림이 있다면 이해하기 쉽다, 원본을 꼭 찾아보길 바란다

다트를 던져서 맞는 공간을 Index 로 본다

이번에는 항상 쓰는 예시 [1/2, 1/3, 1/12, 1/12] 를 선으로 긋지 말고 사각형으로 만들어 보자.
사격형의 가로 ( width ) 는 모든 영역이 동일하다고 가정하고 도형들을 줄 세워보면 "1/2" 가 가장 길고 "1/12" 가 가장 짧다.

4 개의 면 도형을 이어붙이면 큰 사각형이 나온다. 이 사각형은 면적이 w * 4 * h ( h = 1/2 ) 로 볼수 있다.
여기서 w 는 그냥 동일하니 제쳐두고 ( w = 1 ) 총 면적은 4 * h ( 1/2 ) 이며 이 도형에 다트를 하나 던지면 어딘가에는 꽂힐 것이다.

여기서 꽂힌 영역의 도형이 가지는 Index 가 우리가 원하는 결과 값이다.

이를 계산하기 위해 p 를 그대로 height 로 쓰면 복잡하니까 기준 height 를 잡을 건데 이건 normalize 과정이라 생각하면 쉽다. Normalize 하면 probability : [1, 2/3, 1/6, 1/6] 으로 변경된다.

 

이 말은 W 값을 [0, w) 의 RandomValue 를 얻어 w 를 결정한다는 의미다. 예를 들어

w 가 0 인 경우는 p 값이 1 이니 결과가 "0" 이다.
w 가 1 인 경우는 p 값이 2/3 이니 동전 던지기를 해서 얻는 [0, 1) 의 RandomValue 가 2/3 ( p ) 보다 작으면 결과는 "1" 이고 크다면 빈 영역 ( empty ) 이니 이 경우 값은 "1" 이 아닌 것이다.
이런 경우는 w 를 다시 계산해야 한다.
이렇게 heading 되는 영역이 나올 때까지 과정을 반복하면 "Throwing Darts" 가 완성된다.

여기서 문제는 빈 영역이 나오면 다시 w 를 계산해야 하는데 이 때 최악의 반복 횟수 가 n 이 나온다.
* 원문에는 수학적인 증명이 나오니 참고하면 된다

위 방법은 다 좋은데 empty 가 발생하는 문제가 있고 이를 개선한 것이 "The Alias Method" 이다.
"The Alias Method" 는 Dart 에서 빈 영역이 나오지 않게 매꾸는 방법이고이제 전체 사각형을 채우는 방법을 이해하면 된다.

사각형을 채우기 위해서는 기존의 비균등 도형들을 적절하게 짤라 붙여야 하는데 여기서 Avg ( 평균 ) 이 나온다. 쉽게 생각해서 각 w ( 면 ) 이 n 개 라면 각 면적의 확률이 1/n 이 나오도록 height 를 Normalize 이 필요하다.

예를 들어, [1/2, 1/3, 1/12, 1/12] 에서 가장 큰 height ( p ) 는 1/2 이고 이 값이 평균인 "1/4" 이 되려면 "4" 를 곱해줘야 한다. 이제 "4" 값을 각 p 에 적용하면 [2, 4/3, 1/3, 1/3] 을 얻을 수 있다.
이제 가장 큰 height ( 2 ) 의 절반 크기의 height ( 1 ) 를 가지는 사각형에 모든 p 들을 짤라서 채울 수 있다.

 

짤라서 채워지는 공간은 다른 p 가 나올 확률과 같으므로 이를 Mapping 할 수 있는 Table 정보 ( Array ) 가 필요한데 이를 Alias 라고 부른다.

위 경우를 예로 들면
Index : [0, 1, 2, 3]
Prob : [2/3, 1, 1/3, 1/3]
Alias : [1, x, 0, 0]

 

Dart 와 같이 w 를 구하고 해당 w 의 prob 을 동전 던지기로 성공하면 그대로 w 를 쓰고 실패하면 alias 의 값 ( alias[w] ) 를 값으로 쓰면 된다.

prob 과 alias 를 만드는 과정은 아래와 같다.

여기서 prob 자체가 "1 보다 작은 경우" 와 "1 보다 큰 경우" 로 나누는데 그 이유는 prob 자체가 "1 보다 크면" 항상 w ( Index ) 가 나옴을 보장할 수 있고 "1 보다 작으면" 부족한 영역을 채운 ( 짤라서 넣은 ) "1 보다 큰" Index 를 쓰기 때문이다.

• N 크기의 alias, prob 이름의 array 를 만들고 N ( w ) 수 만큼 각 prob 에 곱한 다음에 모든 p 를 순회하며 1 보다 크면 large ( array ) 에 1 보다 작으면 small ( array ) 에 넣는다.
• 이후 small 을 순회하며 small 의 첫번째를 빼고 large 에서 첫번째를 빼고 prob 에는 small 에서 뺀 숫자를 넣고 alias 에서는 large 에서 뺀 Index 를 넣는다
• large 에서는 small 에서 부족한 p 만큼을 빼고 1 보다 크면 다시 large 에 1 보다 작으면 small 에 넣는다.
• 위에서 빠진 small 은 p 가 1 이 되었으므로 그냥 제거하고 끝이다
• 모든 과정을 다 끝내고 large 의 p 가 남아 있다면 그냥 그 자체로 p 가 1 임을 충족하므로 large 에서 빼고 prob 에는 1 을 넣는다

NextValue 를 구하는 과정은 [0, n) 으로 RandomValue 를 구하고 w 에 해당되는 prob 으로 동전 던지기를 해서 성공하면 w 를 그대로 쓰고 아니라면 alias[w] 를 쓴다.

위 방법이 "Vose's Alias Method" 이다.

 

그리고 개선하지 않은 버전은 아래의 특징이 있다.

• "Naive Alias Method" 는 데이터 구조 ( stack, queue, array ) 대신 매 순회마다 p 가 1 보다 작은 경우와 큰 경우를 probability array 에서 찾아 사용하는 방법이다. p 를 찾기 위한 iterator 가 한번 더 필요하기 때문에 alias method 를 만드는 비용이 n^2 이다
• "Alias Method" 는 단일 array 대신 BalancedBinaryTree 를 사용해서 비용이 NlogN 발생한다

마지막으로 실용적인 "Vose's Alias Method" 가 나오는데 이는 부동소수점 이슈가 있기 때문에 디자인을 조금 변경했는데 Table 생성이 조금 다르다.

나머지는 과정은 똑같고 small, 과 large 가 모두 빈 것을 기준으로 순회하며 계산된 p 를 구하는 방법이 "small p + larget p - 1" 로 변경되었다. 그리고 모든 순회 과정을 끝내고 small 이 남는 경우에도 prob 을 1 로 두어 사용한다.
* 부동 소수점 계산으로 인해 small 이 남는 경우도 발생할 수 있다

"Vose's Alias Method" 를 사용하면 Table ( 데이터 구조 ) 생성 비용 Time O(n) / 계산 비용 Time O(1) / Memory 비용 O(n) 을 요구하여 Probability 를 기반으로 임의의 Index 값을 구할 수 있다.