읽기일기

Think Bayes (15) 차원 다루기


파이썬을 활용한 베이지안 통계, 앨런 B. 다우니 지음, 권정민 옮김/한빛미디어

15.1 배꼽 박테리아

배꼽 박테리아 관련 프로젝트의 선행 연구 단계에서 연구자들은 60명의 지원자의 배꼽에서 표본을 추출했다. 여기에서 DNA조각을 추출하여 염기 서열은 분석하면 DNA조각의 종이나 속을 파악할 수 있다. 각 조각을 리드(read)라고 부른다.

이 데이터를 이용하여 아래 질문들에 답을 할 수 있다.

  • 관측된 종 수를 기반으로 하여 생태계에 있는 종 수를 추정한다.
  • 각 종의 확산 정도, 즉 각 종의 비율을 알 수 있다.
  • 샘플을 더 추가 수집했을 때 얼마나 많은 새로운 종이 발견될 것인가.
  • 관측 종의 비율이 주어진 상한선까지 증가하려면 얼마나 더 많은 리드가 필요한가.

이 문제를 보이지 않는 종 문제(Unseen Species Problem)라고 한다.

15.2 사자와 호랑이와 곰

문제를 단순화 해보자. 사자, 호랑이 곰의 3개의 종이 있다고 하였을 때, 야생 동물 보호소에서 사자 3마리, 호랑이 2마리, 곰 1마리를 봤다고 하자. 만약 보호소에서 어떤 동물이든 볼 확률이 동일하다면 우리가 볼 수 있는 각 종의 수는 다항 분포를 생성할 것이다. 각 동물의 확산도를 p_lion, p_tiger, p_bear라고 하면, 사자 3마리와 호랑이 2마리와 곰 1마리를 보는 것에 대한 우도는 아래와 같다.

앞에서 나온 베타 분포를 사용하여보자. 이 방법은 각 종의 개체 분포를 따로따로 나타낸다. 예를 들어, 3마리의 사자와 3마리의 사자가 아닌 동물을 봤다고 표현하면 p_lion의 사후 분포는 다음과 같다.

p_lion의 최대 우도 추정 값은 관측률로 50%이다. 비슷하게 호랑이는 33%, 곰은 17%정도이다.

하지만 여기에는 두 가지 문제가 있다.

  1. 여기에서는 모든 종이 균등 분포를 따른다고 가정하였지만 세 가지 종이 있다는 사실을 알게 되는 순간 사전 분포는 1/3이 된다.
  2. 확산도를 다 더하면 1이므로, 각 종의 분포는 독립적이지 않다.

이를 해결하기 위해 디리클레 분포를 사용하자. 디리클레 분포는 베타 분포를 다차원적으로 일반화 한 것이다. 즉 두 개의 가능한 결과값이 아니라 다양한 수의 결과값을 나타낼 수 있다.

다음은 thinkbayes.py에서의 디리클레 분포 클래스다.

디리클레 분포가 주어졌을 때 각 개체의 확산도에 대한 주변분포는 베타 분포로 다음과 같이 계산할 수 있다.

여기서 i는 구하고자 하는 주변 분포의 인덱스다.

이 예제에서 각 종의 사전 분포는 Beta(1, 2)이다.

예상한대로 각 종의 사전 평균 확산도는 1/3이다. 이제 갱신을 위한 코드를 추가하여 보자. 갱신을 위해서는 변수에 관측값을 더해야한다.

dataparams와 동일한 순서로 기록된 개수에 대한 수열이다.

사후 평균 관측 확산도는 44%, 33%, 22%이다.

15.3 계층 버전

원래 문제로 돌아가서, 전체 종의 숫자를 추정해보도록 하자. 이 문제를 풀려면 다른 스윗을 가지는 메타 스윗을 사용하여야 한다. 최상위 스윗은 전체 종의 숫자에 대한 가설을 가지고 있고, 하위 스윗은 확산도에 대한 가정을 가지고 있다.

가능한 종의 수를 3~30으로 정의하였다. 이 범위 내에서 한 종이 보일 확률은 동일하다고 가정하였다.

계층 모델을 갱신하는 것은 보통 하위 레벨부터 갱신한 후 다음 단계로 넘어가지만 여기서는 최상위 레벨부터 갱신하도록 하자.

data는 관측된 숫자의 수열이고 hypo는 디리클레 객체다. LikelihoodDirichlet.Likelihood를 1000번 호출하고 전체 값을 돌려준다. 실은 Dirichlet.Likelihood는 샘플 하나를 골라내서 우도를 계산하기 때문에 1000번을 반복 수행하여 누적하였다.

15.4 랜덤 샘플링

디리클레 분포에서 임의의 샘플을 만드는 방법은 두 가지가 있는데, 하나는 주변 베타 분포를 사용하는 것이다. 이는 한 번에 한가지만 선택하고 나머지는 더하여 사용하여야 한다.

또 다른 방법은 n개의 감마 분포에서 값을 선택해서 정규화하는 것이다.

결과는 3에서 7까지는 근사값이 유의미 하지만 그 이후로는 확률이 급강하하기 시작한다.

15.5 최적화

속도가 느린 문제를 해결하기 위한 내용을 진행해보고자 한다. 이후 내용은 다음의 것들을 고려한다.

  • 먼저 동일한 데이터로 디리클레 분포를 갱신할 때 처음 m개의 변수는 전체에 대해 동일하다. 유일한 차이점은 가설에서의 보이지 않는 종의 숫자다. 따라서 n개의 디리클레 객체가 필요한 것이 아니고, 계층의 최상위 객체에 이 변수를 저장할 수 있다.
  • 가설 전체에 대해 임의의 값에 대한 동일한 집합을 사용하면 값을 생성하는 시간도 단축디ㅗ기도 하지만, 전체 가설을 샘플 공간에서 동일하게 선택하므로 결과로 수렴하는데 더 적은 횟수의 반복만 하면 된다.
  • 종의 수가 증가할 수록 임의의 확산도에 대한 배열이 커지고 그 중 하나가 선택될 확률은 작아진다. 따라서 대부분의 반복은 전체 합에 도움이 되지 않은 작은 우도를 만들어버릴 것이다.
  • 이를 해결하기 위해 한번에 한 종씩 갱신하고자 한다.
  • numpy 배열 연산자를 사용하여 속도를 높인다.

자세한 내용은 생략한다.

15.9 배꼽 박테리아 데이터

이제 다시 배꼽 박테리아 문제로 돌아오도록 하자. 데이터는 다음과 같은 형식이다.

61개의 종에 대한 400개의 리드 샘플이 있는 B1242이다.

92, 53, 47, 38, 15, 14, 12, 10, 8, 7, 7, 5, 5,
4, 4, 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1

전체에서 큰 비율을 차지하는 몇개의 주요 종이 있지만, 많은 종은 하나의 리드에서만 나타난다. 이러한 단일종(singleton)의 수는 일부 보이지 않는 종이 있을 가능성이 잇다는 것도 알려준다. 따라서 사자와 호랑이 예제에서 처럼 각 박테리아가 리드에 나타날 확률은 동일한다고 가정하였다. 실제로도 수집 프로세스에서 편향성이 나타날 수 있다.

이제 몇가지 배꼽 박테리아를 다뤄보자. 각 세부 내용에 대한 정보를 담고 있는 Subject라는 클래스를 정의했다.

각 Subject는 위에서 본 것과 같은 문자열 코드가 있고, (개수, 종 이름) 쌍으로 된 리스트가 있으며, 이는 개수가 작은 순서부터 큰 순서로 정렬되어있다. 이러한 개수와 종 이름에 쉽게 접근하도록 메서드가 준비되어 있다.

자세한 코드는 http://thinkbayes.com/species.py 에서 참고하자.

Subject에는 n의 분포와 확산도를 나타내는 Species5 스윗을 만들고 갱신하는 Proces메서드가 있다.

또한 Suite2n의 사푸 분포를 구하는 DistN이라는 메서드를 제공한다.

다음으로는 각 종의 확산도의 사후 분포를 구한다

index는 사용할 종을 가리킨다. 각 n에 대하여 가능한 값과 그 확률을 계산한다. 그것을 메타 Pmf를 통하여 결합한다.

15.10 예측 분포

4개의 문제 중 두 문제는 앞에서 n의 사푸 분포와 각 조의 확산도를 계산하여 풀었다.

나머지 두 문제를 풀려면, 사후 분포를 사용하여 가능한 미래 사건을 시뮬레이션 하고 종의 수에 대한 예측 분포를 계산한 후 우리가 보고자 하는 전체 종의 비율을 구하면 된다.

그 과정은 다음과 같다.

  1. 사후 분포에서 n을 고른다.
  2. 디리클레 분포를 통해 보지 못했을 수 있는 종을 포함한 각 종의 확산도를 구한다.
  3. 미래 관측 결과에 대해 임의의 수열을 생성한다.
  4. 추가 리드 수의 비율 k로 새로운 종의 수 num_new를 계산한다.
  5. 앞의 단계를 반복하면서 num_new와 k의 결합 분포를 누적한다.

이에 대한 코드는 다음과 같다. RunSimulation은 단일 시뮬레이션을 실행한다.

num_reads는 시뮬레이션에 사용할 추가 리드의 수다. m은 확인한 종의 숫자고 seen은 각 종의 고유 이름에 대한 문자열 집합이다. n은 사후 분포에서 임의로 선택한 숫자고, observations는 종 이름에 대한 임의 수열이다.

매 반복마다 seen에 새로운 관측 결과를 추가하고 리드와 새로 발견된 종의 수를 기록한다.

RunSimulation의 결과는 리드의 수와 새로 발견되는 종의 수의 쌍으로 이루어진 리스트로 표현되는 희박화 곡선(rarefaction curve)이다.

GetNames에서 종 이름의 리스트를 반환하고 그 이름이 주어지면 SpeciesGenerator에서는 이름 뒤에 숫자를 붙여 종을 만든다. 이름을 다 사용하면 unseen으로 이름이 붙여진다.

cdf는 보이지 않는 종을 포함한 CDF이다.

SamplePosterior는 n값에 따른 확산도 샘플을 생성하는 메서드이다.

15.11 결합 사후 분포

앞의 시뮬레이션을 사용하여 num_newk에 대한 결합 분포를 추정할 수 있고 이를 통해 k에 각 값에 따른 num_new 분포를 구할 수 있다.

curvesRunSimulation으로 생성된 희박화 곡선의 리스트다. 각 곡선은 knum_new의 쌍 리스트를 가진다. 결과 결합 분포에는 각 값과 그에 대한 확률이 연결되어있다.

다음은 ks의 리스트를 받아서 knum_new의 조건 분포를 구한다.

15.12 범위

마지막 질문에 대답하기 위해서는 새로운 종의 수 대신 관측된 종의 비율을 계산해주도록 RunSimulation을 고쳐야 한다.

이제 각 곡선을 생성하는 것을 반복하면서 추가 리드 k의 수를 fracs의 리스트에 연결해주는 딕셔너리 d를 만든다. fracsk개의 리드 확보 후의 범위 값의 리스트이다.

이제 k의 값에 대해 fracs의 CDF를 구한다. 이 CDF는 k리드 후의 범위의 분포를 나타낸다. CDF는 주어진 상한선 이하로 확률이 떨어질 확률을 알려주므로 CDF의 여함수(complementary CDF)는 이를 초과할 확률을 알려준다.


Add a Comment Trackback