Think Bayes (10) 근사 베이지안 계산
 파이썬을 활용한 베이지안 통계, 앨런 B. 다우니 지음, 권정민 옮김/한빛미디어 |
10.1 변이 가설
여기서는 변이 가설(Variability Hypothesis)을 다루고자 한다. 최근 수업에서 CDC의 행동 위험 요인 감시 시스템에서 미국의 남녀 성인이 키를 직접 입력한 데이터를 살펴보았다. 데이터셋에는 15507명의 남성과 254722명의 여성의 응답 내용이 들어있었다.
- 남성의 평균키는 178cm이고 여성의 평균 키는 163cm이다.
- 남성의 표준편차는 7.7cm이고 여성은 7.3cm이다.
- 표준편차를 평균으로 나눈 변동 계수(coefficient of variation, CV)를 계산해보면 남성은 0.0433, 여성은 0.0444이다.
이 데이터셋이 변이 가설을 증명하기에는 빈약하다고 결론 내릴 수 없으나 베이지안 방법을 사용하여 결론을 더 명확하게 만들 수 있다.
- 우선 간단한 구현 내용에서 시작하자. 하지만 이 구현 내용은 1000개 이하의 데이터셋에서만 동작한다.
- 로그 변환하에서 확률을 계산하려면 전체 데이터셋을 스케일 변환을 해야하지만, 계산속도가 느려질 것이다.
- ABC로 알려져있는 근사 베이지안 계산을 통해 속도를 높일 것이다.
코드는 http://thinkbayes.com/variability.py 를 참고한다.
10.2 평균과 표준편차
가우시안 분포의 변수인 평균과 표준편차를 9장과 동일한 방식으로 추정할 것이다.
먼저 각 평균 mu와 표준편차 sigma의 쌍과 이 쌍이 나타날 확률의 연결인 Height스윗을 정의하도록 하자.
|
1 2 3 4 5 6 7 8 |
class Height(thinkbayes.Suite, thinkbayes.Joint): def __init__(self, mus, sigmas): thinkbayes.Suite.__init__(self) pairs = [(mu, sigma) for mu in mus for sigma in sigmas] thinkbayes.Suite.__init__(self, pairs) |
mus는 가능한 mu의 연속들이고, sigmas도 sigma값의 연속이다. 이들은 모두 균등분포이다.
이에 대한 우도 함수는 쉽다. mu와 sigma가 주어졌을 때 특정 값 x의 우도를 계산할 수 있다.
|
1 2 3 4 5 6 7 8 |
class Height: def Likelihood(self, data, hypo): x = data mu, sigma = hypo like = thinkbayes.EvalGaussianPdf(x, mu, sigma) return like |
이 문제의 가장 어려운 부분은 mus와 sigmas의 범위를 고르는 일이다. 범위가 너무 작으면 중요한 확률 값들을 지나칠 수도 있고, 범위가 너무 크면 연산 능력을 낭비하게 될 것이다.
따라서 고전적인 방법을 사용하였다. 이를 위해 표준편차를 사용할 것이다.
만약 n개의 샘플을 이용한 평균과 표준편차의 추정치가 $\mu, \sigma$라면, 이들의 표준오차는 $ s / \sqrt{n}, s / \sqrt{2(n-1)} $이다.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 |
def FindPriorRanges(xs, num_points, num_stderrs=3.0): n = len(xs) m = numpy.mean(xs) s = numpy.std(xs) stderr_m = s / math.sqrt(n) mus = MakeRange(m, stderr_m) stderr_s = s / math.sqrt(2 * (n-1)) sigmas = MakeRange(s, stderr_s) return mus, sigmas |
이들을 이용하여 범위를 만드는 함수는 다음과 같다.
|
1 2 3 4 5 |
def MakeRange(estimate, stderr): spread = stderr * num_stderrs array = numpy.linspace(estimate - spread, estimated + spread, num_points) return array |
10.3 갱신
이제 스윗을 갱신하여보자.
|
1 2 3 4 5 |
mus, sigmas = FindPriorRanges(xs, num_points) suite = Height(mus, sigmas) suite.UpdateSet(xs) print suite.MaximumLikelhood() |
보통 같은 데이터를 두 번 사용하면 허위이긴 하나, 이 경우에는 괜찮다. 사전 분포의 범위에서 매우 작은 값들에 대해서 계산하지 않기 위해 데이터를 선택하엿다. 실제로는 사전 분포는 균등 분포이지만 필요 없는 값들은 모두 무시하자.
10.4 CV의 사후 분포
CV의 분포 계산을 위해 mu와 sigma의 쌍을 만들자.
|
1 2 3 4 5 6 |
def CoefVariation(suite): pmf = thinkbayes.Pmf() for (mu, sigma), p in suite.Items() pmf.Incr(sigma/mu, p) return pmf |
나중에 thinkbayes.PmfProbGreater를 사용해서 남자가 여자보다 변이성이 높을 확률을 계산할 수 있다.
하지만 고려해야 할 문제가 있다.
- 데이터셋의 크기가 커질 수록 소수점의 제한으로 계산 문제가 생긴다.
- 데이터셋에 높은 극단 값이 포함되면 아웃라이어로 작용하여 계산이 올바르게 되지 않을 것이다.
10.5 언더플로
앞의 데이터셋에서 100개의 값을 선택하여 분석하였다면 오류 없이 분석이 될 것이지만, 1000개의 값을 이용해서 프로그램을 돌리면 오류 메시지를 보게 될 것이다.
문제는 우도 계산에서 너무 작은 값들이 계속적으로 곱해져서 0에 수렴해버리게 된다. 이런 경우를 언더플로(underflow)라고 한다.
한가지 해결 방법은 매번 갱신하거나 100개를 다시 돌려서 정규화하는 것이다. 하지만 속도가 매우 느릴 것이다. 더 나은 대안으로는 로그로 변환하여 우도를 계산하는 것이다. 로그로 변환하였기 때문에 곱하는 대신 더하는 것으로 변경된다.
Pmf는 Log를 제공한다.
|
1 2 3 4 5 6 7 8 9 10 |
class Pmf def Log(self): m = self.MaxLike() for x, p in self.d.iterItems(): if p: self.Set(x, math.log(p/m)) else: self.Remove(x) |
가장 높은 확률 m을 찾은 뒤 모든 확률을 이것으로 나누어 가장 높은 확률을 1이 되도록 정규화한다. 따라서 로그화된 확률의 최대값은 0, 나머지는 음수를 갖게 된다. 만약 확률이 0인 데이터가 있으면 이는 제거해버린다.
로그화된 Pmf는 Update, UpdateSet, Normalize를 사용하면 결과가 엉뚱하게 나올 것이므로 사용하면 안된다. 대신 LogUpdate와 LogUpdateSet을 사용하자.
|
1 2 3 4 5 6 7 8 9 10 11 |
class Suite: def LogUpdateSet(self, dataset): for data in dataset: self.LogUpdate(data) def LogUpdate(self, data): for hypo in self.Values(): like = self.LogLikelihood(data, hypo) self.Incr(hypo, like) |
로그 우도를 사용한 뒤에는 다시 Exp를 사용해서 원래대로 되돌리자.
|
1 2 3 4 5 6 7 |
class Pmf: def Exp(self): m = self.MaxLike() for x, p in self.d.iterIetms(): self.Set(x, math.exp(p-m)) |
10.6 로그 우도
|
1 2 3 4 5 6 7 8 |
class Height: def LogLikelihood(self, data, hypo): x = data mu, sigma = hypo loglike = scipy.stats.norm.logpdf(x, mu, sigma) return loglike |
다음은 전체 갱신 프로세스를 실행하는 코드이다.
|
1 2 3 4 5 |
suite.Log() suite.LogUpdateSet(xs) suite.Exp() suite.Normalize() |
하지만 여전히 느리다.
10.7 약간의 최적화
- 가우시안의 로그를 씌울 때 각 데이터에 의존하지 않는 값들을 따로 분리하여 미리 계산하고 의존성이 있는 최소한의 데이터만 반복하여 계산하도록 하였다.
- 반복적으로 걔산되는 값을 캐시에 저장하여 다시 사용하도록 하였다.
10.8 근사 베이지안 계산 (ABC)
하지만 이렇게 일일이 계산하지 않고 근사 베이지안(Approximate Bayesian Compuation) 방법을 사용할 수 있다. ABC는 특정 데이터셋의 우도가 다음과 같을 때 사용된다.
- 큰 데이터셋의 경우로 매우 작아서 로그 변환이 어려운 경우
- 계산 비용이 많이 들어서 최적화를 해야할 경우
사실 이 문제에서는 특정 값의 집합을 볼 확률에 대해서는 관심이 없고, 관측 평균과 분산에 관심이 있다. 만약 원래 분포가 가우시안이라고 가정하면 이미 샘플 통계에서 분포를 정확히 찾아낼 수 있기 때문에 보다 효과적으로 답을 알아낼 수 있다.
사실 앞에서 설명한 것이 바로 그것이다.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 |
def LogUpdateSetABC(self, data): xs = data n = len(xs) m = numpy.mean(xs) s = numpy.std(xs) for hypo in sorted(self.Values()): mu, sigma = hypo stderr_m = sigma / math.sqrt(n) loglike = EvalGaussianLogPdf(m, mu, stderr_m) stderr_s = sigma / math.sqrt(2 * (n-1)) loglike += EvalGaussianLogPdf(s, sigma, stderr_s) self.Incr(hypo, loglike) |
10.9 로버스트 추정
이제 아웃라이어를 처리하는 방법을 알아보자.
4분위 범위를 사용해볼 수 있다. 좀더 일반화 시켜 중간값을 중심으로 범위 p내의 범위만을 사용해볼 수 있다.
|
1 2 3 4 5 6 7 8 9 |
def MedianIPR(xs, p): cdf = thinkbayes.MakeCdfFromList(xs) median = cdf.Percentile(50) alpha = (1-p)/2 ipr = cdf.Value(1-alpha) - cdf.Value(alpha) return median, ipr |
이제 가우시안 CDF를 사용하여 ipr을 sigma의 추정값으로 변화하여 범위 내의 분포를 계산할 수 있다.
def MedianS(xs, num_sigmas):
half_p = thinkbayes.StandardGaussianCdf(num_sigmas) - 0.5
|
1 2 3 4 5 |
median, ipr = MedianIPR(xs, half_p * 2) s = ipr / 2 / num_sigmas return median, s |
이를
10.10 누가 변이성이 높은가?
num_sigmas = 1일 때 CV의 사후 분포를 보았다. 남성의 평균은 0.0410, 여성은 0.0429이다. 따라서 여성이 변이성이 높다고 볼 수 있다. 하지만 num_sigmas = 2를 사용하면 남성이 더 높다는 결론을 볼 수 있다. 남성이 작은 쪽에 치우친 경우가 많은데 num_sigmas를 높일 수록 극단 값의 비중이 높아지기 때문이다.
10.11 토의
ABC에 대한 두가지 해석 방법
- 원 값보다 빠르게 계산하기 위해 사용하는 추정 값
- 주어진 가설하에서 데이터의 우도를 계산하기 위한 대안 모델