베이지안통계분석 전종준 1 1 University of Seoul, Korea Spring 2017 1/49

Outline 들어가며베이지안추론베이지안분석모형 Monte-Carlo Markov Chain (MCMC) 2/49

들어가며 확률변수 동전던지기실험 실험결과는앞면또는뒷면확률변수 X 는앞면일때 1, 뒷면일때 0 의값을갖는다. 주가의로그수익률실험 ( 관찰 ) 결과는주가의로그수익률값확률변수 X 는주가의로그수익률값을갖는다. 3/49

들어가며 확률분포확률분포 : 확률변수의확률이흩어져패턴, 확률 ( 밀도 ) 함수로표현됨. X Bernoulli(θ): 1 과 0 의값에확률이각각 θ, 1 θ 씩흩어져있음. Pr(X = x) = θ x (1 θ) x X N(µ, σ 2 ): Pr(a < X b) = b a ) 1 (x µ)2 exp ( 2πσ 2σ 2 dx 4/49

들어가며 확률분포의요약 확률분포는확률이흩어진패턴즉, 함수로표현되므로숫자로요약할수있다면유용할것이다. 분포의요약치평균, 분산, 왜도, 첨도, 중앙값, 분위수, 최빈값... 베르누이분포는평균혹은 1 이나올확률 (θ) 로완전히요약된다. 정규분포는평균 (µ) 과분산 (σ 2 ) 으로완전히요약된다. 5/49

들어가며 데이터가잘알려진모수분포를따르지않을때분포의성질을어떻게확인할수있을까? 해당분포를따르는랜덤샘플을생성하여그샘플에기반하여평균을구해보거나분산을구해보거나... 히스토그램을통해분포의개형을확인한다. 6/49

들어가며 확률분포의요약 감마분포 : x <- rgamma(n = 1000, shape = 3, scale = 2); summary(x); hist(x) 베타분포 : x <- rgamma(n = 1000, shape1 = 0.5, shape2 = 2); summary(x); hist(x) 7/49

베이지안추론 통계적추론의한방법 추론하는대상을확률변수 ( 어떤확률분포를따르는랜덤한대상 ) 로다룸. 추론해야하는대상의사전확률분포 (prior distribution) 가있음. 데이터가주어진경우주어진사전분포와우도함수를이용하여추론해야할대상에대한사후분포 (posterior) 를구함 구해진사후분포를이용하여의사결정을함. 8/49

베이지안추론 빈도론방법과베이지안방법의비교 빈도론 데이터가특정한분포를따른확률변수라고가정. 그특정분포의모수는알려지지않은고정된상수라고가정 통계적추론은그고정된모수를추정하는것임. 추론에추정량의분포를이용함. 베이지안 데이터가특정한분포를따른확률변수지만, 이미관찰한 ( 주어져있는 ) 상태라고가정. 그특정분포의모수가상수가아니라랜덤 ( 확률변수 ) 하다고가정. 모수에대한분포에대해사전정보가있음. 통계적추론은랜덤한모수의분포에대한추정하는것임. 9/49

베이지안추론 빈도론방법과베이지안방법의비교 : 동전던지기실험 빈도론 X i Bernoulli(θ) (i = 1,, n) θ (0, 1) (constant) ˆθ = X (ˆθ nθ)/ nˆθ(1 ˆθ) 는표준정규분포로근사할수있다. 베이지안 X i θ Bernoulli(θ) θ π( 사전분포 ) θ X 1,, X n? ( 사후분포 ) 사후분포를통해 θ 의추론을한다. 평균 mode 분위수등.. 10/49

베이지안추론 사후분포의유도 : 베이즈정리 데이터 : Y (realization y) 모수 : θ = (θ 1,, θ p ) 우도함수 : L(y θ) ( 모수가주어진경우데이터가생성되는확률모형으로표현됨 ) 사전분포 : π 0 (θ) 사후분포 π(θ y) L(y θ)π 0 (θ) 11/49

베이지안추론 빈도론방법과베이지안방법의비교 : 동전던지기실험동전의앞면이나올확률 θ를추정하는상황을생각하자. 20번의실험을진행하였다. 동전의앞면이 9번나왔다. 12/49

베이지안추론 빈도론자의분석 θ (0, 1) 인알지못하는고정된상수 X i iid Bernoulli(theta) X = 9/20 = 0.45 이므로 ˆθ = 0.45 θ 를약 95% 의확률로포함하는신뢰구간은 (ˆθ 1.96SE(ˆθ), ˆθ + 1.96SE(ˆθ)) = (.232,.668) 이다. ( 단, SE= 0.45 0.55/20 =.111 ) 13/49

베이지안추론 베이지안의분석 θ 는 (0, 1) 위에서값을같는확률변수 θ π 0 = Beta(10,10) (Beta(10,10) 는모수가 10, 10 인베타분포 ) X i θ iid Bernoulli(θ) π(θ x) = 상수 L(x θ) π 0 = Beta(19, 21) 14/49

베이지안추론 베이지안의분석 Beta(19,21) 분포의개형 : 평균 0.487, 2.5% quantile 0.333, 97.5% quantile 0.642 15/49

사전분포와사후분포의비교 사전분포 : 평균 0.5, 2.5% quantile 0.289,.97.5% quantile 7113 사후분포 : 평균 0.487, 2.5% quantile 0.333, 97.5% quantile 0.642 16/49

베이지안추론 동전던지기실험결과비교 빈도론자 ˆθ = 0.45 95% 신뢰구간 (.232,.668) 베이지안 π(θ x) Beta(19, 21) 17/49

유용한모형들 포아송 - 감마모형 우도함수포아송모형 : 0 을포함한자연수의값에서확률을갖는모형예 : 사고건수, 특허건수... 양수위에서정의된평균모수 1 개있음 (µ). 감마 ( 사전 ) 분포양수위에정의된분포함수형태모수 (α), 규모모수 (β) 가있음. 18/49

유용한모형들 포아송-감마모형데이터 : x i 는 X i µ iid Poisson(µ) (i = 1,, n) 의실현값우도함수 : 포아송확률밀도함수의곱사전분포 : µ π 0 = Gamma(α, β) 사후분포 : n µ x Gamma(α + x i, (n + β 1 ) 1 ) i=1 19/49

유용한모형들 포아송 - 감마모형 1 년동안발생하는특허분쟁건수가 20 건으로관측되었다. 데이터는포아송모형으로가정하고모형모수의사전분포를 Gamma(4, 4) 라고가정하자 n = 1 이고 x 1 = 20 따라서, π(µ x) Gamma(4 + 20, (1 + 1/4) 1 ) 20/49

유용한모형들 포아송 - 감마모형 21/49

유용한모형들 지수-감마모형데이터 : x i 는 X i µ iid Exp(λ) (i = 1,, n) 의실현값우도함수 : 지수분포의확률밀도함수의곱사전분포 : λ π 0 = Gamma(α, β) 사후분포 : n θ x Gamma(α + n, ( x i + β 1 ) 1 ) i=1 22/49

유용한모형들 정규 - 정규모형 데이터 : x i 는 X i µ iid N(µ, σ 2 ) (i = 1,, n) 의실현값 우도함수 : 정규분포의확률밀도함수의곱 사전분포 : µ π 0 = N(m, s 2 ) 사후분포 : µ x N(m, s ) m = w x + (1 w)m ( 단 w = ns 2 /(ns 2 + σ 2 )), 분산이 (s ) 2 = ( n σ 2 + 1 s 2 ) 1 23/49

사후분포의분석 예제 : 정규 - 코시모형 Let Y 1,, Y n i.i.d. N(θ, 1) π 0 (θ) = 1/π(1 + θ 2 ) Posterior: π(θ y) n i=1 exp{ i θ) 2 } 1 2 1 + θ 2 n(θ ȳ)2 exp{ } 1 2 1 + θ 2 사후분포의개형이표준형태가아니기때문에분포의특성을파악하기힘들다. 24/49

사후분포의분석 사후분포의특성값사후분포의평균 : E(θ y) 사후분포의분산 : Var(θ y) θ 의 credible set {a(y), b(y)}: Pr{a(y) < θ < b(y) y} = 0.95 25/49

Monte Carlo method 일반적인문제적분의계산 E π [h(x)] = h(x)π(x)dx : 어려운문제만약아래와같은성질을가지는독립표본을추출할수있다고하자. X (1), X (2),, X (N) π(x) 그렇다면다음과같은방법으로원하는적분값을근사할수있을것이다. E π [h(x)] h N = 1 N Monte Carlo 적분이라부른다. N h(x (t) ) t=1 26/49

Monte Carlo method Remark 만약앞서말한독립표본을추출할수있다면대수의법칙에의해다음이성립하여이론적인근거를확보하게된다. h N E π [h(x)] (1) as N 하지만 π 분포를따르는독립표본을얻는것은일반적으로어렵다. 독립표본을얻는대신표본평균이참평균으로수렴하는 (ergodic) 랜덤표본을추출하는방법을생각한다. 27/49

Monte Carlo method 상태전이모형인 Markov chain 에서극한분포의성질을활용하고자한다. Markov chain 은상태전이에대한확률모형이고, 전이확률이현재의상태에만의존한다는성질이있다. 여러가지상태를옮겨다니는어떤확률적규칙 ( 전이확률 ) 이주어져있을때, 그규칙들이좋은성질을만족시키면, 전이상태를매우많이한이후부터는각상태의상대방문빈도가좋은성질을가지는유일한확률분포로수렴하게된다. 즉, 극한확률이우리가원하는분포 π 를따르게전이확률구조를잘만들면, Markov chain 을통해랜덤샘플을생성하는것이가능하다. 이러한방법을 MCMC(Monte-Caro Markov Chain) 이라고한다. 28/49

Monte Carlo method 예를들어어떤확률변수가값을 1,2,3,4 만가지고각각의확률이 3/16, 1/4, 14/48, 13/48 이라고하자. 이산형확률변수이므로독립표본추출이매우쉽지만다음같은방법으로표본을추출해보자 다음과같은전이행렬을준비한다. P = 1/4 1/4 1/2 0 0 1/4 1/2 1/4 1/4 1/4 1/4 1/4 1/4 1/4 0 1/2 전이행렬의 i 행, j 열원소를 p ij 라고하자. p ij 의미는상태 i 에서상태 j 로옮겨갈 ( 전이 :transition) 확률이다. 29/49

Monte Carlo method 먼저출발점을 1 이라고하자 (X (0) = 1) 상태 1 에서 (p 11, p 12, p 13, p 14 ) 의확률로상태를옮겨간다. p 14 = 0 이므로상태 1 에서상태 4 로옮겨갈수는없다. 상태 3 에왔다고하자 (X (1) = 3). 다시여기서 (p 11, p 12, p 13, p 14 ) 의확률로상태를옮긴다. 이작업을반복하면아주많은시행 N 번후에 X (N+1), X (N+2), 는정상극한분포로수렴하며그확률분포가 1, 2, 3, 4 에서각각의확률 3/16, 1/4, 14/48, 13/48 가된다. 즉, 전이확률구조를이용하여원하는분포를가지는표본을추출할수있다. 원하는분포로수렴하는극한분포를가진전이확률을찾아내는것이중요한문제임. 30/49

깁스표본수집 이변량분포에서 MCMC 를알아본다. (X, Y ) π(x, y) π(x, y) 로부터랜덤표본을추출하는것이어렵다고하자. 하지만다음같이조건부확률이잘알려져있어 X Y = y π(x y) 과 Y X = x π(y x) 의표본수집이쉬운상황을가정하자 이때깁스표본수집법은극한에서정상분포 π 를가지는 Markov chain 을생성할수있다. 31/49

깁스표본수집 알고리즘 1 Initialization : Set X (0) = x (0) and Y (0) = y (0) 2 For i = 1 to N 1 Generate X (i) π(x y (i 1) ) 2 Generate Y (i) π(y x (i) ) 32/49

깁스표본수집 Remark (X (1), Y (1) ), (X (2 ), Y (2 )), 은정상인 π(x, y) 극한분포를가지는 Markov chain 이다. 각점들이상태를의미하고선이상태의변화를의미한다. 깁스표본수집방법은상태전이확률이조건부확률로주어진 Markov Chain 을생성해 π 를따르는랜덤샘플을추출한다. Figure: 깁스표본수집방법에서상태의전이과정 33/49

깁스표본수집 예제 Y i i.i.d N(µ, σ 2 ), π(µ, σ 2 ) 1 σ 2. 사후분포 : π(µ, σ 2 y) ( 1 σ 2 )n/2+1 exp{ τ = 1/σ 2 이라하면다음을보일수있다. π(µ σ 2, y) = N(ȳ, σ 2 /n) π(τ µ, y) = Gamma(n/2, (y i µ) 2 /2) (yi µ) 2 2σ 2 } 두조건부확률이모두잘알려진표준형분포이기때문에깁스표본수집으로 MCMC 를수행하면된다. 34/49

깁스표본수집 일반적인케이스에서깁스표본수집 1 Initialization: Set X 1 (0) = x 1 (0),, X p (0) = x p (0) 2 For i = 1 to N 1 Generate X (i) 1 π(x 1 x (i 1) 2,, x (i 1) p ) 2 Generate X (i) 1 π(x 2 x (i) 1, x (i 1) 3,, x (i 1) p ) 3 Generate X (i) 1 π(x 3 x (i) 1, x (i) 2, x (i 1) 4,, x (i 1) 4 5 Generate X (i) p π(x p x (i) 1,, x i) p 1 ) p ) 35/49

MH 알고리즘 Metropolis-Hastings algorithm 사후분포의형태는 π(θ y) L(y θ)π 0 (θ) 이며, 많은경우 π(θ y) = ConstantL(y θ)π 0 (θ) 의상수항 (constant) 를계산하기힘들다. 이때조건부확률을구하기힘드므로깁스표본수집을사용하기어렵다. 36/49

MH 알고리즘 예제 : 정규 - 코시모형 Y 1,, Y n i.i.d. N(θ, 1) π 0 (θ) = 1/π(1 + θ 2 ) 사후분포 : π(θ y) n i=1 exp{ i θ) 2 } 1 2 1 + θ 2 n(θ ȳ)2 exp{ } 1 2 1 + θ 2 θ π(θ y) 를생성하고자한다. 37/49

MH 알고리즘 알고리즘 1 Choose a transition function q(y x) of a certain Markov chain 2 Initialize x (0). 3 For i = 1 to N 1 Generate x q(x x (i 1) ) 2 With probability α(x (i 1), x) = min{1, π( x)q(x (i 1) x) π(x (i 1) )q(x x (i 1) ) set x (i) = x (acceptance) else set x (i) = x (i 1) (rejection) 38/49

MH 알고리즘 Remark 사용된전이확률모형은 α(x, y)q(y x) 다. π(x) 의비만사용하므로 π 를구성하는상수항을몰라도알고리즘을적용할수있다. q(y x) 는실제로샘플을생성해야하는함수므로사용하기쉬운확률밀도함수를사용한다. 이론적으로는 q( x) 과 π( ) 이같은구간에서확률을가지기만하면 M-H 알고리즘이잘작동하지만, 실제로는 q 를선택하는것은매우중요한문제다. 39/49

MH 알고리즘 모형추정시고려사항수렴성 : Markov chain 의정상분포로수렴여부가중요함 MCMC sample chain 에서앞부분을버린다 : Burn-in MCMC sample 을몇개씩건너뛰면서표집한다 : thinning 서로다른초기값을가진다수의 MCMC chain 을비교한다 40/49

베이지안회귀분석 선형회귀분석 Y i = x i β + ɛ i for i = 1,, n ( 단, x i, β R k 벡터, ɛ 1,, ɛ n iid N(0, σ 2 )) 전통적으로는 β 의추정량으로 LSE 를사용했지만베이지안방법은 β 의분포를추정하고자한다. 기본적으로베이지안방법에서는알지못하는대상은모두확률변수로놓는다. 41/49

베이지안회귀분석 β에대한사전분포 : 다변량정규분포 σ 2 에대한사전분포 : 역감마분포데이터 : (y i, x i ) (i = 1, n) (β, σ) {(y i, x i ) : i = 1, n} 의분포에관심이있음 42/49

베이지안회귀분석 회귀모형계수의사후분포의추정 43/49

베이지안회귀분석 일반화선형모형에대한회귀분석 R-package MCMCpack 에서제공하는모형 : 모수에대한사전분포의설정이제한적임 로지스틱회귀분석 프로빗회귀분석 포아송회귀분석 44/49

베이지안회귀분석 변화점분석 데이터 : 정규분포를따름데이터의평균모수가변화에대한확률적추론을하고싶음변화점의개수어떤시점에변화점이있을확률 45/49

베이지안회귀분석 변화점분석 46/49

실습 : 직접해보기 예제 : Boston Housing 자료 Boston Housing 자료는보스턴의 506 개교외지역의주거와관련된수치들을담고있는자료이다. 자료에대한간략한설명은 Help 메뉴의 Boston 을통하여확인할수있다. 이데이터는 MASS 패키지에포함되어있다. 47/49

실습 : 직접해보기 R 의 MASS 패키지로부터 Boston 자료를불러와이중 chas 와 rad 변수를제외하고 Boston.data 라는이름으로저장하자. 48/49

실습 : 직접해보기 crim 변수를종속변수로하고나머지를설명변수로하는베이지안회귀모형을적합하시오. ( 단, 이때 MCMCpack 패키지를이용하고적절한 burn-in period 와 iteration 횟수를사용하시오.) summary() 를이용하여결과를확인하고결과를해석해보시오. plot() 을이용하여각변수의마코프체인과사후분포를확인해보시오. 다른난수를이용하여베이지안회귀모형을새로적합하고, gelman.diag() 명령어를이용하여각변수에따른 MCMC 의수렴여부를확인하시오. 49/49