Think Bayes (6) 의사 결정 분석
파이썬을 활용한 베이지안 통계, 앨런 B. 다우니 지음, 권정민 옮김/한빛미디어 |
6.1 ‘그 가격이 적당해요’ 문제
미국에는 ‘그 가격이 적당해요’라는 게임 쇼가 있다. 이 쇼는 진열대 위에 물건들이 있고, 그 물건들의 실제 가격의 합보다 높지 않은 가격을 부르면서도, 그 가격 합과 가장 가까운 금액을 부르는 참석자가 이기는 게임이다.
베이지안 사고방식을 가진 사람이 볼 때 이 시나리오는 몇 가지 질문을 던지게끔 한다.
- 결과를 보기 전에, 진열대의 가격에 대해 참가자들은 어떤 사전 믿음이 있을까?
- 결과를 본 후, 참가자들의 이런 믿음은 어떻게 바뀌었을까?
- 사후 분포 기반으로 참가자들은 예상가를 얼마나 적어야 할까?
세번째 질문이 바로 베이지안 분석의 의사 결정 분석을 나타낸다. 사후 분포가 주어지면, 참가자들의 기대값을 최대화하는 금액을 결정할 수 있다.
여기서 사용한 코드는 http://thinkbayes.com/price.py에서, 데이터 파일은 http://thinkbayes.com/showcases.2011.csv, http://thinkbayes.com/showcases.2012.csv에서 받을 수 있다.
6.2 사전 분포
사전 분포에 대해서는 ‘그 가격이 적당해요’의 예전 데이터들을 사용해보도록 하자. 여기서는 2011, 2012년 데이터를 사용하도록 하겠다.
6.3 확률 밀도 함수
지금까지는 확률 누적 함수나 확률 질량 함수(PMF)를 다루엇다. PMF는 가능한 각 값과 그것이 나올 확률을 연결한 도표와 같은 것이었다. 이제 확률 밀도 함수(PDF)라는 것을 사용하도록 하자.
수학적으로 PDF를 표기하면 함수 형식으로 나타낼 수 있다. 잘 알려진 것으로 평균이 0, 표준 편차가 1인 가우시안 분포인 PDF를 예로 들 수 있다.
$$f(x) = \frac{1}{\sqrt{2 \pi}} \exp (-x^2 / 2)$$
이 함수는 주어진 x에 대한 확률 밀도를 계산한다. 밀도 자체는 확률이 아니지만(?) 밀도를 연속적 범위로 적분한다면 그 결과값으로 확률이 나올 것이다.
6.4 PDF 나타내기
thinkbayes.py
는 Pdf
라는 클래스를 제공한다. 이 클리스 역시 추상형으로 되어있다. Pdf
는 Density
와 MakePmf
두가지 메서드를 포함한다.
1 2 3 4 5 6 7 8 9 10 |
class Pdf(object): def Density(self, x): raise UnimplementedMethodException() def MakePmf(self, x): pmf = Pmf() for x in xs: pmf.Set(x, self.Density(x)) pmf.Normalize() return pmf |
Density
는 값 x를 받아서 이에 해당하는 확률 밀도를 반환한다. MakePmf
는 PDF의 이산 추정 값을 계산한다.
PDF를 상속받아 만들어진 GaussianPdf
가 있다.
1 2 3 4 5 6 7 |
class GaussianPdf(Pdf): def __init__(self, mu, sigma): self.mu = mu self.sigma = sigma def Density(self, x): return scipy.stats.norm.pdf(x, self.mu, self.sigma) |
여기서는 scipy.stats
의 함수를 사용하도록 되어있다.
실생활에서 가우시안에 근접한 분포를 갖는 수치들이 많아 유용하게 사용할 수 있다. 하지만 실제 데이터의 분포를 가우시안이나 다른 간단한 수학 함수로 나타낼 수 잇다는 보장이 없다. 이런 경우에는 전체 개체의 PDF를 추정하기 위해 샘플을 골라서 사용한다.
예를 들어, ‘그 가격이 적당해요’의 데이터는 1회차에 대해 313개의 금액이 있다. 한가지로 가능한 모든 회차의 가격에서 샘플을 추출하는 것을 고려해볼 수 있다. 다음은 샘플의 예이다.
1 2 |
28800, 28868, 28941, 28957, 28958 |
샘플을 보면 28801부터 28867까지의 값은 없지만 실제로 여기에서 값이 존재하지 않는다는 근거는 없다. 따라서 PDF 곡선이 부드럽게 이어지도록 유추하여야 할 것이다.
커널 밀도 추정(KDE)는 샘플을 이용하여 실제 데이터에 적합한 추정 평활 PDF를 찾아내는 방법이다. thinkbayes
에서는 scipy의 KDE를 이용하여 이를 구현해 놓았다.
1 2 3 4 5 6 |
class EstimatedPdf(Pdf): def __init__(self, sample): self.kde = scipy.stats.gaussian_kde(sample) def Desity(self, x): return self.kde.evaluate(x) |
__init__
에서 샘플을 입력받으면 이에 대하여 커널 밀도 추정값을 계산한다. Density
에서는 만들어진 kde
객체에서 추정된 밀도 값을 반환한다. 다음은 예제코드이다.
1 2 3 4 5 6 7 |
prices = ReadData() pdf = thinkbayes.EstimatedPdf(prices) low, high = 0, 75000 n = 101 xs = numpy.linspace(low, high, n) pmf = pdf.MakePmf(xs) |
6.5 참가자 모델링
만약 이제 쇼에 참가할 사람이라면, 각 회차의 가격에 대한 사전 믿음의 값의 분포를 살펴볼 수 있을 것이다. 사전 값을 보기 전에 다음 질문을 생각해보자.
- 어떤 데이터를 고려해야 하고 어떤 데이터를 수집해야할까?
- 우도 함수를 계산할 수 있을까? 즉 각 price에 대한 가설 값을 정하면서 데이터의 상태 우도를 계산할 수 있을까?
이 질문을 위해 참가자 모델을 만들기로 했다. 이는 틀린 정보가 주어졌을 때 이를 기반으로 금액을 추측하는 용도로 사용하고자 한다. 참가자가 상품을 보고, 각 상품의 금액을 추정하고 이렇게 구해진 금액을 더한다. 이를 총 guess
라 하자.
그렇다면 실제 가격이 price
일때 참가자가 guess
라고 추정했으면 이 때 우도는 어떻게 될까? 이를 다시 표현해보자.
error = price - guess
이라고 하고, 참가자들의 배당 금액 추정 값을 error
로 구한 우도는 어떻게 될까?
데이터를 이용하여 참가자들의 배당 금액과 회차의 실제 가격 간의 차이를 나타내는 diff = price - bid
에 대한 누적 분포가 나타나 있다. diff
가 음수일 경우, 예상가가 실제 값보다 높아지는 것이다.
여기서는 error
의 분포가 평균 0, 분산은 diff
와 동일한 가우시안을 따른다고 가정해보자.
Player
클래스에서 이 모델을 구현했다.
1 2 3 4 5 6 7 8 9 10 11 12 |
class Player(object): def __init__(self, prices, bids, diffs): self.pdf_price = thinkbayes.EstimatedPdf(prices) self.cdf_diff = thinkbayes.MakeCdfFromList(diffs) mu = 0 sigma = numpy.std(diffs) self.pdf_error = thinkbayes.GaussianPdf(mu, sigma) def ErrorDensity(self, error): return self.pdf_error.Density(error) |
prices
, bid
, diffs
는 각각 각 회차에 대한 값과 입찰 금액, 그 차이 값들을 순서대로 가지고 있는 값이다.
ErrorDensity
는 주어진 error
값에 대하여 PDF를 제공한다.
6.6 우도
이제 우도를 작성해보자.
1 2 3 4 5 6 7 8 9 10 11 12 13 |
class Price(thinkbayes.Suite): def __init__(self, pmf, player): thinkbayes.Suite.__init__(self, pmf) self.player = player def Likelihood(self, data, hypo): price = hypo guess = data error = price - guess like = self.player.ErrorDensity(error) return like |
hypo
는 예상 가격이다. data
는 참가자들이 가격을 잘 추정했을 때의 값이다. ErrorDensity
를 통하여 error
에 대한 우도를 계산한다.
6.7 갱신
Player
에 참가자들이 추측을 통해 사후 분포를 계산하는 메서드를 추가하자.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 |
class Player: # 초기화 부분에 추가 n = 101 self.prices_xs = numpy.linspace(0, 75000, n) def MakeBeliefs(self, guess): pmf = self.PmfPrice() self.prior = Price(pmf, self) self.posterior = self.prior.Copy() self.posterior.Update(guess) def PmfPrice(self): return self.pdf_price.MakePmf(self.price_xs) |
PmfPrice
는 사전 분포 값을 위하여 가격의 PDF 이산 추정 값을 생성한다.
사후 확률을 계산하는 방법은 다음과 같다. 우선 각 가설에 대한 Likelihood
(우도)를 계산하고, 여기에 사전 분포를 곱한다. 다음 정규화를 위해 Update
를 호출한다.
6.8 최적 입찰
최적 예상가를 계산하여 보자
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 |
class GainCalculator(object): def __init__(self, player, opponent): self.player = player self.opponent = opponent def ExpectedGains(self, low=0, high=75000, n=101): bids = numpy.linspace(low, high, n) gains = [self.ExpectedGain(bid) for bid in bids] return gains def ExpectedGain(self, bid): suite = self.player.posterior total = 0 for price, prob in sorted(suite.Items()): gain = self.Gain(bid, price) total += prob * gain return total def Gain(self, bid, price): if bid > price: return 0 diff = price - bid prob = self.ProbWin(diff) if diff <= 250: return 2 * price * prob else: return proce * prob def ProbWin(self, diff): prob = (self.opponent.ProbOverBid() + self.opponent.ProbWorseThan(diff)) return prob |
1 2 3 4 5 6 |
class Player: def ProbOverbid(self): return self.cdf_diff.Prob(-1) def ProbWorseThan(self, diff): return 1 - self.cdf_diff.Prob(diff) |
player
와 opponent
는 Player
객체다.
ExpectedGains
는 각 입찰금과 입찰금에 대한 예상 이익의 배열을 생성한다.
ExptecedGain
은 예상가의 예상 이익을 계산한다. 이 메서드는 사후 값을 반복해서 생성한면서, 주어진 회차의 실제 가격에 대한 이익과 그 확률에 대한 가중치를 곱한 합을 구하여 예상 이익을 계산한다.
Gain
은 은 입찰 금액과 실제 가격을 인수로 받는데, 너무 높게 입찰하였으면 0을, 아니라면 입찰금과 실제 가격의 차이를 계산하여 예상 이익을 구한다.
망냑 상대방이 입찰을 높게했다면 이기지만, 그게 아니라면 diff
보다 차이가 더 많이 나기를 바래야 한다.
마지막으로 최적의 예상가를 구하는 코드는 다음과 같다.
1 2 3 4 5 6 7 8 |
class Player: def OptimalBid(self, guess, opponent): self.MakeBeliefs(guess) calc = GainCalculator(self, oppnent) bids, gains = cals.ExpectedGains() gain, bid = max(zip(gains, bids)) return bid, gain |
6.9 토의
베이지안의 추정은 결과가 사후 분포로 나온다는 것이다. 고전 추정법은 단일 점 추정이나 신뢰구간 추정으로 이루어진다는 점에서 차이가 있다. 이는 추정 값을 연관된 다른 분석의 입력 값으로 사용하려 할 때는 큰 도움이 되질 않는다. 베이지안 사고를 처음 하는 사람은 종종 사후 확률을 평균과 최대 우도 추정값만 계산해서 요약하려는 경향이 있으나, 이 값만이 필요하다면 초반부터 베이지안 방법을 쓸 필요가 없을 것이다.