예제

2014년 6월에 있었던 제 6대 전국 지방선거에서는 지방 자치 단체장을 선거로 뽑았습니다. 그 당시 JTBC에서는 투표일 직전 현대리서치와 오베이에게 의뢰하여 투표 의향과 지지 후보를 물어보았고 그 결과는 경기도의 경우 다음과 같았습니다.

구역	유효 표본수	남경필 지지 (%)	김진표 지지 (%)	구역	유효 표본수	남경필 지지 (%)	김진표 지지 (%)
과천/군표/의왕	60	53.8	41.7	평택	55	53.8	46.2
성남	129	56.6	43.4	화성	53	45.5	54.5
수원	152	53.4	46.6	고양	138	45.2	54.8
안성/용인	106	43.6	56.4	김포	40	48.2	51.8
안양	81	63.8	36.2	파주	38	24.5	75.5
광명	39	23.4	76.6	동두천/양주/포천/연천	42	68.9	31.1
부천	104	53.6	46.4	의정부	55	36.6	63.4
시흥	41	41.9	58.1	가평/양평/여주/이천	53	66.2	33.8
안산	79	50.5	49.5	광주/구리/하남	71	39.5	60.5
오산	32	41.4	58.6	남양주	71	42.7	57.3

위 자료를 이용한 선거 예측의 정확성을 높혀주기 위해서 2년전에 있었던 대선 투표 결과를 사전정보로 반영한 베이즈 추정을 구현하는것이 하나의 방법일수 있습니다. 어떻게 하면 될까요?

풀이

이 문제는 전형적인 베이즈 추정 문제입니다. 각 구역 $i$에서의 참값을 $Y_i$라고 하면 그에 대한 통계 추정량값은 $\hat{Y}_i$으로표현될수 있습니다. 그리고 2년전 대선에서 그 지역의 참값 (새누리당 대선후보 투표율)을 $X_i$라고 하면 $Y_i$ 와 $X_i$이 어느 정도 상관관계가 있을 것이라고 생각할수 있습니다. 이러한 관계는 일종의 사전모형이 될 것입니다. 즉, 사전모형이란 현재 자료수집 이전의 정보에 대한 모형으로 이 경우에는 \begin{equation} Y_i = \gamma X_i + e_i , \ \ \ e_i \sim N( 0, X_i \sigma^2) \ \ \ (1) \end{equation} 으로 표현할수 있고 이 경우 $\gamma$와 $\sigma^2$은 사전모형에 대한 모수가 됩니다.

다음으로는 현재 자료를 통해서 얻어지는 통계량의 모형을 구해야 합니다. 이 경우에는 \begin{equation} \hat{Y}_i\mid Y_i \sim N( Y_i , V_i) \ \ \ \ \ \ (2) \end{equation} 으로 표현될수 있고 $V_i$는 $\hat{Y}_i$의 분산값으로써 위의 경우에는 $(1- \hat{Y}_i)\hat{Y}_i/ n_i$으로 계산될 수 있습니다.

위의 두 식을 결합하면 베이즈 정리를 이용하여 사후 모형이 다음과 같이 얻어집니다. \begin{equation} Y_i \mid ( X_i , \hat{Y}_i ) \sim N [ \alpha_i \hat{Y}_i + (1-\alpha_i) \gamma X_i, (1-\alpha_i) X_i \sigma^2 ] \ \ \ \ \ \ \ \ (3) \end{equation} 이때 $\alpha_i = \frac{ X_i \sigma^2}{X_i \sigma^2 + V_i }$ 으로 계산됩니다. 이렇게 얻어진 사후 모형의 조건부 기대값 \begin{equation} \hat{Y}_i^* \equiv E ( Y_i \mid X_i , \hat{Y}_i ) = \alpha_i \hat{Y}_i + (1-\alpha_i) \gamma X_i \end{equation} 은 서베이 결과와 과거 대선 결과를 반영하는 예측으로써 사전 모형의 유용성 정도에 따라 예측의 정확성을 높여주는 효과를 가져옵니다. 위의 조건부 기대값 $\hat{Y}_i^*$이 베이즈 추정량이라고도 부르는데 실제로는 이 베이즈 추정량이 사전모형 모수의 함수이므로 이를 추정하여 대입해서 계산해야 합니다.

토론

1. 위의 예에서 (3)의 사후 모형이 정규분포로 얻어진 것은 (1)의 사전 모형도 정규분포이고 (2)의 관측치 모형도 정규 분포를 따르기 때문입니다. (엄밀하게 말하면 (2)는 추정량의 표본분포입니다. ) 만약 (1)의 사전 모형이 정규 분포가 아닌 다른 모형이면 사후 분포는 정규 분포를 따르지 않고 \begin{equation} p(Y_i \mid \hat{Y}_i, X_i) = \frac{ f( \hat{Y}_i \mid Y_i) \pi (Y_i \mid X_i)}{ \int f( \hat{Y}_i \mid Y_i) \pi (Y_i \mid X_i) d Y_i} \end{equation} 를 따르게 됩니다. 이러한 사후 분포가 형태가 알려져 있지 않은 경우에는 몬테카를로 방법을 사용해서 베이즈 추정량을 계산합니다.

2. 모수 추정은 EM 알고리즘이나 베이지안 MCMC 방법을 사용합니다. 만약 EM 알고리즘을 사용하는 경우에는 모수 $\gamma$와 $\sigma^2$에 대한 score equation 을 $Y_i$와 $Y_i^2$의 함수로 표현한 후에 모형 (3)에서 얻어진 사후 기대값을 계산해서 반복적으로 풀어야 합니다. (디테일은 생략하겠습니다) 또는 (1)과 (2)를 결합하여 \begin{equation} \hat{Y}_i \sim N[X_i \gamma , V_i + X_i \sigma^2 ] \end{equation} 을 얻을수 있으므로 $(X_i , \hat{Y}_i)$의 관측치를 바탕으로 최대우도 추정법을 이용하여 모수를 추정할 수 있습니다. 이렇게 모수를 최대우도 추정으로 계산한후 $\hat{Y}_i^*$ 에 대입해서 구하는 것을 Empirical Bayes 추정이라고 합니다.

3. 아래 그림은 위의 방법으로 예측한 결과를 실제 투표 결과와 비교한 오차 절대값의 분포입니다. 왼쪽은 서베이 결과만을 사용한 예측 분석의 오차 분포이고 오른쪽은 서베이 결과와 대선자료를 베이즈 추정으로 결합한 예측 분석의 오차 분포입니다. 서베이 결과는 평균 7% 의 오차를 보여주었는데 베이즈 추정은 평균 2%의 오차를 보여주었습니다. 이 경우 베이즈 추정으로 정확도가 상당히 개선되었다고 이야기 할수 있을 것입니다.

. .

예제

어느 수입농산물 관리업체에서는 수입농산물의 품질을 매달 측정하고 그에 대한 통계를 생산해서 관리하는 시스템을 구축하려고 합니다. 그런데 그에 대한 테스트가 A,B 방법 두가지 종류가 있는데 A 방법은 싸지만 정확성이 떨어지고 B방법은 비용은 비싸지만 매우 정확한 방법이라고 합니다. 이러한 경우 제한된 예산 하에서 어떤 식으로 샘플을 뽑아야 비용 대비 정확성이 높은 조사 방법이 될수 있을까요? 예를들어 A 방법은 1만원이 들고 B 방법은 10만원이 든다고 할때 전체 예산 1천만원하에서 어떻게 효율적인 통계생산 시스템을 만들수 있을까요?

풀이

세가지 방법을 생각할수 있습니다. 하나는 1,000개의 샘플에 A 방법을 적용하는 것입니다. 그 방법은 현실적으로 적절하지 않습니다. 왜냐하면 A 방법 자체의 부정확성 때문에 생산되는 통계가 편향의 위험이 있기 때문입니다. 다른 하나는 그냥 100개의 샘플에 B 방법을 적용하는 것입니다. 이것은 조사에 따른 편향은 없지만 샘플수가 적어서 통계량의 분산이 커지고 통계의 대표성이 떨어질수 있습니다.

마지막 방법은 Two-phase sampling 이라고 불리우는 것입니다. 이 방법은 위의 두 방법의 절충인데 예를 들어 500개의 샘플을 먼저 뽑아서 방법 A 를 적용하고 그 중에서 50개를 랜덤하게 뽑아서 방법 B 를 적용한후 50개의 최종표본으로부터 통계적 모형을 이용하여 A 방법의 편향을 추정한후 첫번째 표본 500 개의 편향을 보정한 예측값의 평균을 사용하는 방법입니다.

좀더 자세히 설명하기 위해서 A방법으로 측정된 값을 $X$ 라고 하고 B 방법으로 측정된 값을 $Y$라고 하고 이 두 변수간에 다음의 모형을 세울수 있다고 합시다. $Y = \beta_0 + \beta_1 X + e$ 이때 $E(e)=0$ 으로 가정합니다. 이러한 경우 $E(Y \mid X)= \beta_0 + \beta_1 X$ 는 $X$의 편향을 제거한 값이 됩니다. 따라서 50개의 최종 샘플로부터 그 추정치 $\hat{\beta}_0$와 $\hat{\beta}_1$를 구하게 되면 $\hat{Y} = \hat{\beta}_0 + \hat{\beta}_1 X$를 사용하여 500개의 전체 표본에서 편향이 제거된 $Y$의 예측치를 구할수 있습니다. 이를 이용해서 500개의 원자료에서 $\hat{Y}_i$값을 계산한후 그값의 평균으로 통계를 생산하는 방법이 Two-phase sampling 을 사용한 추정 방법이 됩니다.

이러한 two-phase sampling 추정법은 $X$와 $Y$의 상관관계가 높을수록 그 효과가 큽니다. 이 상관관계를 $\rho$라고 표현하면 two-phase sampling 추정량의 분산은 \begin{equation} \frac{1}{500} \sigma_y^2 \rho^2 + \frac{1}{50} \sigma_y^2 ( 1- \rho^2) \end{equation} 으로 얻어집니다. 위의 식에서 $\rho^2$이 1에 가까울수록 위의 분산은 작아지게 됩니다. 예를들어 $\rho^2=0.8$라고 하면 위의 분산은 $0.0056 \sigma_y^2$으로서 두번째 방법을 사용했을 경우 표본수 178개를 사용하는 것과 같은 효과를 얻습니다. 조사 비용이 1,780만원이 드는 조사를 1,000만원으로 줄이는 방법이 되는 것입니다.

토론

1. 위의 내용을 좀더 일반화하면 다음과 같습니다. 먼저 $n_1$ 개의 샘플에 A방법을 적용하여 $X$값을 측정하고 그 중에서 일부를 뽑아 $n_2$ 개의 표본을 얻어 그 샘플에서 $Y$값을 측정하는 경우 two-phase sampling 추정량의 분산은 \begin{equation} V= \frac{1}{n_1} \sigma_y^2 \rho^2 + \frac{1}{n_2} \sigma_y^2 ( 1- \rho^2) \end{equation} 으로 표현합니다. 그런데 전체 비용은 $C_1 n_1 + C_2 n_2 =1000$ 만원으로 계산할수 있는데 여기서 $C_1$은 A 방법을 사용할 때 드는 비용 (1만원), $C_2$는 B 방법을 사용할때 드는 비용(10만원)이 됩니다. 이 제약조건 하에서 위의 $V$를 최소화하는 최적화 문제를 생각할수 있습니다. 그 해는 코쉬 부등식을 사용하여 풀면 그 최적해는 \begin{equation} \frac{n_1}{n_2} = \sqrt{\frac{1-\rho^2}{\rho^2} \times \frac{C_1}{C_2}} \end{equation} 으로 구해집니다. 이를 $C_1 n_1 + C_2 n_2=1000$만원에 대입하면 최적 배분이 얻어집니다.

2. 위의 예제에서처럼 $\rho^2=0.8$이라고 하면 그 최적값은 $n_1=390$, $n_2=61$이 되어서 그 분산이 $0.00533\sigma_y^2$으로 되고 따라서 위의 두번째 방법을 사용했을때의 표본수 188개를 사용하는 것과 같은 효과를 갖습니다. 표본수가 188개를 사용하면 1,880만원을 들어야 하는데 1,000만원으로 동일한 효과를 얻게 되는 것입니다. 이렇게 two-phase sampling 은 비용을 줄이고 최대의 효과를 보이는데 유용한 통계학적 방법입니다. 만약 매달 이 통계를 생산한다면 1년에 880*12 = 10,560만원을 절약하는 것입니다.

3. 이러한 two-phase sampling 으로는 $X$값은 다 관측되지만 일부 표본에서 $Y$값이 결측이 되는 자료 구조를 얻습니다. 이렇게 표본 설계 단계에서 고의로 결측자료를 생성하여 비용을 줄이고 효율을 높힐수 있는데 이러한 방법을 계획결측(planned missingness) 이라고도 부릅니다. 보건 통계에서는 질병 유무를 나타내는 변수인 $Y$를 먼저 관측하고 그에 따라 일부 표본에서 $X$를 관측하는 방식인 case control study 가 많이 사용되는데 이 역시 계획결측의 일종입니다.

. .

서론

오일러의 등식이라고 들어보셨나요? \begin{equation} e^{i \pi} + 1 = 0 \end{equation} 으로 표현되는 등식은 종종 수학에서 가장 아름다운 수식이라고도 불리우는데요 여기에는 5가지 아주 중요한 숫자가 나옵니다. 이는 바로 $e, \pi , i , 1, 0$의 다섯개 숫자입니다. 이 다섯개의 숫자는 수학에서 아주 중요한 의미를 갖는 수인데요 이 다섯개가 모두 사용되어서 하나의 식을 만들어 내고 있으니 수학도들에게 큰 매력을 주는 수식입니다.

자, 통계학에서도 위의 오일러의 등식에 해당되는 수식이 있다면 그건 무엇일까? 그렇다면 그 다섯가지 숫자에 해당되는 다섯가지 주요 개념은 무엇일까요?

풀이

하버드 대학의 Meng 교수는 작년 12월 상하이에서 열린 ICSA 학회 기조연설에서 수학의 오일러 공식에 대응되는 공식으로 다음과 같은 수식을 사용할 것을 제안했습니다.
\begin{equation} \hat{\mu} - \mu = \rho \sigma \sqrt{ \frac{ N-n}{n}} \ \ \ \ \ (1) \end{equation} 이 수식 하나에 통계학의 가장 중요한 기호 5가지가 나온다는 것입니다. 첫째는 $\mu$ 입니다. 평균에 해당되는 기호입니다. (그 앞의 $\hat{\mu}$은 $\mu$에 대한 추정량을 의미합니다.) 두번째는 $\sigma$ 입니다. 표준 편차에 해당되는 기호입니다. 다음으로는 $\rho$입니다. 그 기호는 상관계수를 나타내는 기호입니다. 여기서는 y와 관측확률과의 상관계수를 지칭합니다. 네번째로는 $n$입니다. 표본 크기를 나타냅니다. 마지막으로 등장하는 기호가 $N$입니다. 표본이 뽑힌 모집단의 크기를 나타내는 기호입니다. 이렇게 5개의 주요 기호가 합쳐져서 하나의 수식을 나타내니 이것이 오일러의 등식에 해당된다는 이야기입니다.

토론

1. 위의 수식에 대한 좀더 자세한 설명은 다음과 같습니다. 현실 세계에 존재하는 모집단은 유한 모집단이므로 이를 ${x_1, x_2, \cdots, x_N}$이라 표현하고 관심 모수를 $\mu_g = N^{-1} \sum_{i=1}^N g( x_i)$라고 한다면 이 관심모수를 얻어진 자료의 단순 평균으로 $\hat{\mu}_g = \sum_{i=1}^N R_i g(x_i)/(\sum_{i=1}^N R_i )$ 을 사용하여 $\mu_g$를 추정하는 문제를 생각할수 있다. 이때 $R_i$는 모집단 원소 $i$가 샘플에 포함되어 $x_i$값을 관측하게 될수 있을때 1의 값을 갖고 그렇지 않을때 0의 값을 갖는 지시변수를 의미한다. 이러한 경우 Uniform $\{1,\cdots, N \}$의 분포를 갖는 확률변수 $I$라는 것을 생각하면 다음 식이 유도된다. \begin{equation} \hat{\mu}_g - \mu_g = \frac{E_I [ R_I g(x_I) ] }{E_I (R_I)}- E_I [ g (x_I) ] = \frac{ Cov_I [ R_I, g(x_I)]}{E_I (R_I)} = Corr_I [R_I, g(x_I)] \times
   \sqrt{V_I [ g(x_I)]} \times \frac{ \sqrt{V_I (R_I)}}{E_I(R_I)} \end{equation} 여기서 $E_I(R_I)=n/N$이고 $V_I (R_I) = (n/N)(1-n/N)$이므로 이를 간단히 표현하면 $\hat{\mu}_g - \mu_g = \rho_{R,g} \times \sigma_g \times \sqrt{\frac{N-n}{n}}$ 으로 표현되는데 이로써 식 (1)이 얻어진다.

2. 위의 논리를 좀더 발전시키면 다음과 같은 식이 얻어집니다. \begin{equation} E[ \hat{\mu}_g - \mu_g]^2 = D_I \times D_O \times D_U. \end{equation} 여기서 $D_I = E[ \rho_{R,g}^2]$으로 표현되는데 Defect index 라고 불리우는데 Data quality 를 나타내는 지표이고, $D_O =(N-n)/n= (1-f)/f$로써 Dropout Odds 라고 부를수 있는데 Data quantity 를 나타내는 지표이고, $D_u=\sigma_g^2$은 Data uncertainty 라고 부를수 있는데 이는 problem difficulty 를 나타내는 지표라고 이해할수 있습니다. 통계의 핵심은 결국 data quality, data quantity, and data uncertainty 로 요약될수 있습니다.

3. 위의 공식을 바탕으로 빅데이터의 유효표본수를 구할수 있습니다. 즉, 현재 주어진 빅데이터의 표본수와 동일한 MSE (평균오차제곱)을 구현하는 단순임의표본의 표본수는 얼마인가를 계산하는 것입니다. (유효표본수와 관련된 설명은 저의 이전 포스팅 “유효표본수”를 참고하시기 바랍니다.) 공식은 다음과 같습니다. (여기서 $f=n/N$입니다) \begin{equation} n_{eff} = \frac{n }{1+ (1-f) [ (N-1) D_I-1] } \cong \frac{f}{1-f} \frac{1}{D_I} \end{equation} 따라서 $\rho_{R,g} = 0.05$ 라고 한다면 (약 5%의 선택편향이 있다면) $f=1/2$일 경우 $n_{eff}=400$을 얻게 됩니다. 즉, 5% 정도의 약한 선택편향이 있는 경우, 모집단의 반이 응답을 선택하였더라도 이때 얻어지는 추정의 정확성은 제대로된 확률표본 $400$명을 추출한 것과 동일하다는 것입니다. 예를들어서, 이러한 자발적 표본은 서울시 인구의 반인 500만명이 응답을 했더라도 실제 정확도는 400명의 확률표본을 얻어내 조사한 것과 마찬가지입니다. 이 경우 5%의 선택편향이 미치는 영향은 자료의 99.99%를 소용없게 만드는 것입니다.

4. 따라서 우리가 표본수가 크면 클수록 정확할 것 같지만 선택편향이 있는 경우에는 이를 무시한 신뢰구간이 참값을 벗어갈 확률이 1에 가깝게 수렴하기에 표본이 크면 클수록 신뢰구간이 더 참값에서 멀어진다는 것입니다. Meng 교수는 이를 두고 Big data paradox 라고 불렀습니다.

. .

예제

2014년에 콩고 공화국의 어느 지역에는 에볼라 바이러스 증상을 보이는 환자가 나타났습니다. 그래서 콩고 공화국에서는 급히 국제 보건 기구의 도움을 받아 그 주변 지역에서 백신 접종 캠페인을 실시 했습니다. 이 캠페인은 열흘 동안 열리는데 그 캠페인을 통해 바이러스 백신 접종률을 높이는게 아주 중요하다고 합니다. (만약 접종률이 낮으면 전염병이 다른 지역으로 퍼질 확률이 높아지게 된다고 합니다.) 따라서 백신 접종률을 높히기 위해서 다양한 노력을 하는 것도 중요하지만 실제 접종률을 정확히 파악하는 것 역시 중요합니다.

실제 접종률을 파악하기 위해 그 지역의 마을을 다니며 서베이를 했는데 현지 사정상 제 4일차, 5일차, 6일차, 8일차에 서베이를 해서 마을 주민들의 백신 접종률을 아래와 같이 추정했다고 합니다.

날(day)	추정 접종율(%)	표준 오차(%)	접종자수
4	68.8	15.3	240,624
5	75.6	14.1	274,701
6	87.5	10.6	305,807
8	93.8	7.6	341,780
10			366,783

위의 테이블에서 표준오차(Standard error)는 서베이를 통해서 얻어지는 추정치의 신뢰도를 나타내는 지표입니다. 위 테이블의 접종자수는 서베이를 통해서 얻어지는 것이 아니라 캠페인 기간동안 접종을 받은 사람들 숫자를 집계한 것입니다. 이 접종자수에는 해당 마을 주민 뿐만 아니라 인근 마을 주민들도 접종을 했기에 이를 바탕으로 접종률 통계를 산출할수는 없다고 합니다. 위의 정보를 이용해서 결국 마지막날인 제 10일째에 실제 접종률이 얼마였는지를 추정하는 것이 자료분석의 목표입니다. 어떻게 해야 보다 정확한 추정을 할수 있을까요?

풀이

이러한 자료 구조는 실제 문제에서 많이 나타납니다. 우리가 원하는 변수는 실제 접종률인데 (이를 $Y$로 표현합시다) 이를 직접 관측하지는 못하고 $Y$의 추정량인 $\hat{Y}$만을 관측합니다. 다행히 확률 표본 추출로 얻어진 서베이에서는 추정량의 표준오차를 얻을수 있습니다. 위의 자료에서 접종자수는 $Y$와 관련있는 보조정보로써 $X$라고 표현할수 있습니다. 우리는 상식적으로 $Y$와 $X$에 선형관계가 있을 것이라고 모형을 세울수 있습니다. 이러한 자료에서 흔히 가정되는 모형은 ratio model 입니다. \begin{equation} Y_i = X_i \gamma + e_i, \ \ \ e_i \sim N (0, X_i \sigma^2) . \ \ \ \ \ (1) \end{equation} 이러한 모형에서 $\gamma$의 최적 추정량은 $ \hat{\gamma} = (\sum_i Y_i )/( \sum_i X_i )$ 으로 계산됩니다. 이를 바탕으로 $Y_{10}$의 예측을 $\hat{Y}_{10} = X_{10} \hat{\gamma}$ 으로 계산할수 있습니다. 이러한 추정을 비추정(ratio estimation)이라고 부릅니다.

그러나 문제는 우리가 $Y_i$를 관측하는 것이 아니라 $\hat{Y}_{i}$을 관측하므로 위의 비추정을 바로 적용할수 없다는 것입니다. 만약 $u_i = \hat{Y}_i - Y_i$이라고 표현하면 이 분포는 평균이 $0$ 이고 분산이 $v_i$인 정규분포로 가정할수 있으므로 $\hat{Y}_i=Y_i + u_i$의 주변분포는 $\hat{Y}_i \sim N( X_i \gamma , X_i \sigma^2 + v_i )$으로 표현됩니다. 우리 자료에서는 4개의 관측치가 있고 2개의 모수가 있으므로 ($v_i$는 알려져 있음) 모수 추정이 가능해 집니다. 이렇게 해서 모수가 추정되면 $Y_{10}$의 예측치는 $\hat{Y}_{10}=X_{10}\hat{\gamma}$으로 구해집니다. 위의 자료를 가지고 계산하면 $\hat{Y}_{10}=1.024$로 나오게 됩니다.

토론

1. 모수 추정을 위해서 최대우도 추정법을 사용하거나 Iterative reweighted least squares method 를 사용할수 있는데 후자는 $\gamma$와 $\sigma^2$을 번갈아가면서 추정하는 것입니다. 이를 위해서 $W_i (\sigma^2)=1/( X_i \sigma^2 + v_i )$으로 정의하면 $\gamma$의 추정은 $\sum_i W_i ( \hat{\sigma}^2 ) ( \hat{Y}_i - X_i \gamma)^2$ 을 최소화 하는 값인 \begin{equation} \hat{\gamma}= \frac{\sum_i W_i ( \hat{\sigma}^2) X_i \hat{Y}_i }{ \sum_i W_i ( \hat{\sigma}^2) X_i^2 } \label{1} \end{equation} 으로 구하면 됩니다. $\sigma^2$의 추정을 위해서는 $Z_i (\gamma)= (\hat{Y}_i - X_i \gamma)^2$으로 정의하면 $Z_i(\gamma) \sim [ X_i \sigma^2 + v_i , 2( X_i \sigma^2 + v_i )^2 ] $임을 이용하여 $\sum_i \left\{ \frac{ Z_i (\hat{\gamma}) }{ X_i \sigma^2 + v_i } - 1 \right\}^2$ 을 최소화하는 $\sigma^2$을 찾으면 됩니다. 이는 \begin{equation} \hat{\sigma}^{2(t+1)}= \frac{\sum_i { W_i ( \hat{\sigma}^{2(t)} ) }^2 X_i ( \hat{Z}_i - v_i) }{ \sum_i { W_i ( \hat{\sigma}^{2(t)} ) }^2 X_i^2 } \label{2} \end{equation} 으로 계산됩니다. 따라서 위의 두 식을 번갈아가면서 계산하면 그 수렴값이 최적 모수 추정값으로 사용될수 있습니다.

2. 위의 방법대로 추정을 하면 ${Y}_{10}$의 예측값이 1보다 큰 값을 가지게 될수 있습니다. 이는 모형 (1)이 $Y_i \in (0,1)$의 range restriction 을 부여하지 않았기 때문에 발생하는 것으로 이해할수 있습니다. 이를 해결하기 위해서는 다음의 모형을 세울수 있을 것입니다.
   \begin{equation} Y_i = p(\gamma_0+ X_i \gamma_1) + e_i, \ \ \ e_i \sim (0, p_i (1-p_i) \sigma^2) . \end{equation} 이때 $p(x)= \exp(x)/\{1+ \exp(x)\}$ 으로 표현되고 $p_i = p(\gamma_0 + X_i \gamma_1)$입니다. 따라서 $\hat{Y}_i \sim \left[ p(\gamma_0 + X_i \gamma_1) , v_i + p_i (1-p_i) \sigma^2 \right]$ 이 얻어지므로 토론 1에서 소개한 Iterative reweighted least squares method 방법을 사용하여 모수를 추정할수 있습니다. 이렇게 해서 모수 추정을 하게 되면 $Y_{10}$의 예측은 $\hat{Y}_{10} = p(\hat{\gamma}_0 + X_{10} \hat{\gamma}_1)$으로 계산됩니다. 위의 자료를 적용하면 $\hat{Y}_{10}$=95.2%이 얻어집니다.

3. 또한 이러한 모수 추정을 바탕으로 $\hat{Y}_i$가 관측된 날에 대하여 $Y_i$에 대한 보다 개선된 예측이 가능해질 것입니다. $Y_i$에 대한 추정량으로써 $\hat{Y}_i$와 $\hat{p}_i = p(\hat{\gamma}_0 + X_{i} \hat{\gamma}_1)$의 두가지를 생각할수 있는데 전자는 직접 조사에 의한 추정량으로 직접추정량(direct estimator)이라고 부르고 후자를 모형에 의하여 얻어진 추정량으로 모형기반 추정량이라고 부를수 있습니다. 이 두개를 최적 결합한 예측값은 \begin{equation} \hat{Y}_i^* = \hat{\alpha}_i \hat{Y}_i + (1- \hat{\alpha}_i) \hat{p}_i \label{5} \end{equation} 으로 계산되고 이때 최적계수 $\hat{\alpha}_i$는 \begin{equation} \hat{\alpha}_i = \frac{ \hat{p}_i (1-\hat{p}_i) \hat{\sigma}^2 }{v_i + \hat{p}_i (1-\hat{p}_i) \hat{\sigma}^2} \end{equation} 으로 계산됩니다.

4. 위의 예제는 제가 했던 컨설팅에서 얻어진 실제 데이터를 사용한 것입니다. 위의 방법론을 바탕으로 각 날짜별로 백신 접종률의 추정값을 구하면 다음과 같습니다. 괄호 안의 값은 해당 추정량의 표준오차에 해당하는 값이 됩니다. 모든 경우에 대해서 최적추정의 표준 오차가 제일 작은 것을 확인할수 있을 것입니다. 이는 모형을 반영하여 일종의 shrinkage estimation 을 구현했기 때문입니다.

날(day)	직접 추정(%)	모형기반추정(%)	최적추정(%)
4	68.8 (15.3)	66.6 (15.0)	68.0 (10.7)
5	75.6 (14.1)	78.9 (12.9)	77.2 (9.6)
6	87.5 (10.6)	87.0 (10.7)	87.2 (7.6)
8	93.8 (7.6)	93.0 (8.1)	93.1 (5.7)

. .

. .

예제

김재광 교수는 지난 학기에 카이스트에서 통계학 개론 과목을 강의했습니다. 그런데 홍길동(가명)이라는 학생이 기말고사 기간에 미국에서 열리는 학회에 참석해야 해서 시험을 치룰수 없게 되었다고 하면서 어떻게 해야 할지 알려달라는 문의 메일을 보내온 것이었습니다. 김재광 교수는 중간고사 성적만으로 그 학생의 기말 고사 점수를 예측해서 그 예측점수의 90%에 해당하는 값을 그 학생의 기말고사 점수로 대체해서 최종 성적을 내겠다고 공지했습니다. 이럴 경우 어떻게 해야 합리적인 방법으로 홍길동 학생의 기말 고사 성적 예측을 할 수 있을까요?

풀이

다른 분야에서도 그렇겠지만 특히 통계학에서는 지금까지 관측된 자료를 바탕으로 일어나지 않은 사건에 대한 예측을 하는 것이 주된 관심사 중의 하나입니다. 위의 예제처럼 중간고사 성적을 $X$ 라고 하고 기말고사 성적을 $Y$라고 한다면 다른 학생들의 $X$와 $Y$의 관계를 이용해서 통계적 모형을 세운후 이 모형을 바탕으로 홍길동 학생의 특정 $X$값으로부터 $Y$ 값을 예측할수 있는데 이러한 통계 분석 기법을 회귀 분석이라고 부릅니다.

실제로 이 과목의 중간고사 평균은 $61.38$였고 기말고사 평균은 $64.11$이었습니다. 그리고 이 홍길동 학생의 중간고사 점수는 35점 이었습니다. 그렇다면 위의 문제는 $X=35$일때 $Y$의 조건부 기대값이 얼마이겠냐는 문제로 귀결됩니다. 이러한 문제를 어떻게 풀수 있을까요?

한가지 접근법은 다른 학생들 중에서 홍길동 학생과 가장 비슷한 중간고사 성적을 가진 학생의 기말고사 점수를 사용하는 것일 겁니다. 만약 김철수라는 학생이 중간고사에서 똑같은 $35$점을 받았다면 이 학생이 중간고사 시점 당시에는 홍길동 학생과 실력이 비슷하다고 볼수 있으므로 김철수 학생의 기말고사 성적으로 홍길동 학생의 기말고사 성적을 예측해 볼수 있을 것입니다. 이러한 방법은 $X$ 를 기준으로 가장 가까운 이웃을 찾아서 그 값을 사용한 것이므로 nearest neighbor 방법이라고도 부릅니다. 그런데 이러한 방법을 사용하면 한개의 이웃값만을 사용하므로 조금 곤란한 일이 발생할수 있습니다. 만약 김철수 학생이 기말에서 무단으로 결시를 해서 0점을 받았다고 합시다. 그러면 홍길동도 0점 처리 해야 할까요? (그건 합리적인 방법이 아니겠지요. )

그러한 위험성을 없애기 위해서는 이웃의 범위를 조금 넓히고 싶을 것입니다. 성적이 가장 가까운 학생 1명보다는 점수가 가까운 학생 $K$명(예: $K=5$)을 사용해서 그들의 평균을 사용하는 것이 조금 더 안전할 것입니다. 이를 KNN (K-th nearest neighbor) 방법이라고 합니다. 그렇지만 $K$를 너무 크게 잡으면 그것 또한 문제가 생기겠지요. K를 작게 하면 예측값의 분산이 커지고 K 를 크게 하면 예측값의 편향 제곱값이 커지게 됩니다. 최적의 $K$를 찾아내는 것은 통계학 연구의 주된 주제입니다.

만약 중간고사와 기말고사의 관계식이 평균적으로 일차식의 관계가 있다고 가정하면 회귀 모형을

$$Y = \beta_0 + \beta_1 X + e$$

으로 가정할수 있고 여기서 $e$는 $X$로 설명되지 않는 부분을 나타내는 것으로써 흔히 오차라고 부릅니다. 이러한 방법은 모형을 $\beta_0$, $\beta_1$ 같은 유한개의 모수(parameter)로 결정되도록 표현하고 문제를 모수를 전체 자료를 사용하여 결정하는 문제로 바꾸어 줌으로써 (이를 추정이라고 합니다) 예측치의 효율을 높혀주고자 하는 방법입니다. 이렇게 유한개의 모수로 모형을 표현하여 추론하는 방법을 모수적 접근법 (parametric approach)라고 하고 앞의 KNN 처럼 직접적인 모형 가정을 하지 않는 방법을 비모수적 접근법 (nonparametric approach) 이라고 부릅니다. KNN을 사용한 비모수적 방법에는 K 개의 관측치가 사용되지만 모수적 방법에는 N 개의 관측치가 사용됩니다. 따라서 모수적 방법이 전체 자료를 다 사용하는 방법이므로 관측된 정보를 최대한 활용한다고 할수 있지만 모수적 가정을 해야 하니 위험이 따릅니다. 각기 장단점이 있는 것이지요.

토론

1. 비모수적 접근법에서 모형에 대한 가정은 $Y=f(x)+ e$로 표현하되 $f(x)$가 $x$에 대한 연속함수라는 정도의 약한 가정입니다. 이 가정의 특수한 경우가 $f(x)=\beta_0+ \beta_1 x$가 될 것입니다. 따라서 비모수적 접근법이 보다 약한 가정을 사용한 방법론이므로 더 넓게 적용될 수 있다고 볼수도 있습니다. 모수적 접근법은 가정된 모수적 모형이 어느 정도는 정확하게 자료의 분포를 기술해 주고 있다는 가정 하에서 그 예측의 정당성이 뒷받침 됩니다. 모수적 접근법은 더 많은 가정을 바탕으로 하였기에 그 가정이 맞을 경우에는 비모수적 접근법으로 얻어진 예측보다 더 정확한 예측을 해 줍니다.

2. 만약 중간고사를 한번만 본게 아니라 여러번 보았다면 비모수적 접근법은 보다 어려워집니다. 이러한 경우에는 $X$가 벡터가 되므로 그것에 대한 이웃(neighbor)를 정의하기가 어려울 뿐만 아니라 $X$의 차원이 커질수록 (중간고사 횟수가 많아질수록) 모든 $X$ 값에 대해 일정 범위안에 있는 가능한 이웃의 숫자가 급격히 줄어듭니다. 예를 들어 중간고사 성적이 5점 차이 이내에 있는 경우를 이웃으로 정의한다면 중간고사 횟수가 많아지면서 모든 시험내에서 성적이 5점 차이를 같는 이웃의 숫자는 급격하게 줄어들게 됩니다. 이러한 문제를 차원의 저주 (curse of dimensionality) 라고 부릅니다. 이 차원의 저주 때문에 비모수적 접근법은 $X$의 차원이 커지면 (표본수 N 이 아주 크지 않는한) 현실적으로 불가능한 방법론이 됩니다.

3. 위의 사례는 작년에 실제 일어났었던 일입니다. 자료분석을 하여 얻어진 회귀 분석 결과는 다음과 같습니다. $\hat{Y} = 48.95 + 0.247 * X$ 그리고 이때 $\hat{\beta}_1$의 표본 오차(standard error) 추정값은 0.00673 값이 나왔습니다. 즉, 중간고사 성적이 높은 사람이 기말고사 성적이 높은 경향이 있다는 것을 알수 있습니다. (아래 plot 을 참조하세요.) 그래서 홍길동의 중간고사 점수 $x=35$점을 위의 식에 대입하면 $\hat{Y}=57.59$ 가 나옵니다. 아래 그림에서 살펴보면 알수 있듯이 $x=35$점은 분명 아주 낮은 점수 중에 속하는데 위의 회귀분석 예측을 통해 계산된 57.59 라는 점수는 기말고사의 낮은 점수대에 속한다기 보다는 그보다 좀더 평균값에 가까운 예측을 하게 됨을 보여줍니다. 이러한 현상은 종종 평균으로의 회귀(regression to the mean) 이라고 불립니다. 따라서 위와 같은 단순 회귀 분석으로 예측을 해주면 중간고사에서 꼴찌는 유리해지고 1등은 불리해지는 결과를 갖게 됩니다. 그래서 저는 이 방법을 사용하지 않고 순위(rank)를 고려하는 다른 방법을 사용했습니다. 이렇게 통계 모형을 사용할 때에는 특정 도구에 대해 맹목적인 확신을 갖기 보다는 보다는 조심스럽게 접근하는 자세를 갖는것이 더 바람직합니다.

. .

. .

예제

어느 벤쳐회사가 온라인에서 영화를 추천해주는 시스템을 개발하는 경우를 생각해 봅시다. 그 시스템 안에는 K 개의 영화가 있고 고객은 한편당 얼마씩 지불하며 영화를 스트리밍으로 볼수 있는데 영화 감상이 끝나면 영화에 대한 간단한 추천 여부를 클릭해야 한다고 합시다. (만약 중간에 영화 상영을 중단하면 비추천인 것으로 간주한다고 합시다.)

이 회사에서는 고객들의 영화 감상 취향을 파악하여 로그인한 고객들에게 기존 데이터를 바탕으로 그들의 영화 취향을 파악한후 맞춤형으로 영화를 추천하는 것을 목적으로 할때 어떤 식으로 추천 시스템을 구축해야 할까요?

풀이

이러한 문제는 전형적인 결측 자료 (missing data) 예측 문제입니다. $Y_k$를 영화 $k$를 추천하는지를 나타내는 여부를 나타내는 지시변수 (추천이면 1의 값을 갖고, 비추이면 0의 값을 갖는 변수)라고 정의하면 각각의 고객에게서 우리는 $Y_k$을 (1 또는 0로) 관측하거나 또는 $Y_k$를 관측하지 못합니다. 각각의 고객에게서 $Y_1$, $Y_2$, $….,$ $Y_K$ 에 대한 정보를 일부 수집하게 되는데 이중 많은 값들은 결측(missing)이 됩니다.

예를들어 간단히 $K=5$이라고 합시다. 즉, 이 회사에는 5개의 영화만 상영합니다. 그런데 홍길동이라는 사람은 이중에서 3개의 영화를 보았고 (영화 1,2,3 을 보았다고 합시다) 그 사람의 추천값이 $Y_1=1$, $Y_2=0$, $Y_3=1$ 라고 합시다. 그렇다면 그 사람에게 영화 4를 추천하는게 좋을까요 아니면 영화 5를 추천하는게 좋을까요?

이에 대답하기 위해서는 현재의 고객 데이터로부터 전체 영화들의 선호도 관계를 추정한후 $P(Y_4=1 \mid Y_1=1, Y_2=0, Y_3 =1 )$와 $P(Y_5=1 \mid Y_1=1, Y_2=0, Y_3 =1 )$의 조건부 확률을 계산하여 이 두개의 조건부 확률 중에서 더 큰 값을 갖는 항목에 해당하는 영화를 추천해 주면 되는 것입니다. 이러한 조건부 확률을 추정하기 위해서는 EM 알고리즘이라는 방법론을 사용하는 것이 정석입니다.

이러한 EM 알고리즘을 사용하기 위해 먼저 전체 변수의 결합 확률 $P(Y_1=a, Y_2=b, Y_3=c, Y_4=d, Y_5=e):=\pi_{abcde}$을 정의합니다. 이 결합확률로부터 조건부 확률을 구할수 있습니다. 예를 들어 $Y_4$와 $Y_5$가 결측인 경우에 대해 \begin{equation} P( Y_4=d , Y_5=e \mid Y_1=a, Y_2=b, Y_3=c) = \frac{ \pi_{abcde}}{\sum_{d=0}^1 \sum_{e=0}^1 \pi_{abcde} } \ \ \ \ \ \ (1) \end{equation} 를 구할수 있을 것입니다. 만약 이 조건부 확률값이

$$P( Y_4=d , Y_5=e \mid Y_1=a, Y_2=b, Y_3=c) = \left\{ \begin{array}{ll} 0.6 & \mbox{ if } (Y_4, Y_5)=(1,1) \\ 0.2 & \mbox{ if } (Y_4, Y_5)=(1,0) \\ 0.1 & \mbox{ if } (Y_4, Y_5)=(0,1) \\ 0.1 & \mbox{ if } (Y_4, Y_5)=(0,0) \\ \end{array} \right.$$

으로 얻어진다면 이 값을 바탕으로

고객	영화1	영화2	영화3	영화4	영화5
홍길동	추천	비추천	추천	?	?

의 고객 데이터를

고객	영화1	영화2	영화3	영화4	영화5	가중치
홍길동	추천	비추천	추천	추천	추천	0.6
홍길동	추천	비추천	추천	추천	비추천	0.2
홍길동	추천	비추천	추천	비추천	추천	0.1
홍길동	추천	비추천	추천	비추천	비추천	0.1

으로 가상적으로 만들어낼수 있을 것입니다. 여기서 가중치란 결국 조건부 확률을 지칭합니다. 만약 모든 항목에 대해 관측이 된 케이스는 가중치가 1이 되는 것입니다. 만약 김재동 이라는 고객이 영화 2와 영화 4에 대해서만 평가를 했다고 하면 이에 대한 조건부 확률은

$$P( Y_1=a , Y_3=c, Y_5=e \mid Y_2=b, Y_4=d) = \frac{ \pi_{abcde}}{\sum_{a=0}^1\sum_{c=0}^1 \sum_{e=0}^1 \pi_{abcde} } \ \ \ \ \ \ \ (2)$$

으로 얻어지게 됩니다. 이 경우 김재동의 가상데이터는 $(Y_1, Y_3, Y_5)$의 모든 가능한 값에 대한 8개의 조합을 바탕으로 구현하는데 이에 대한 가중치는 (2)에서 얻어진 조건부 확률을 사용합니다. 이런 식으로하면 모든 고객에 대해 가상 데이터를 만들수 있는데 이는 각 고객에게 모든 가능한 경우의 조합으로 자료를 확장하고 그에 대한 가중치를 구현해 줌으로써 전체 자료에 결측이 없는 것처럼 만들어 주는 것이 됩니다. 이를 E-step 이라고 합니다.

이를 이용하여 만들어진 전체 가상 자료를 이용하여 $\pi_{abcde}$의 추정치를 간단하게 업데이트 할수 있습니다. 이 추정에는 가중치를 반드시 사용해야 합니다. 즉, $\pi_{abcde}$의 추정치는 ($Y_1$=a, $Y_2$=b, $Y_3$=c, $Y_4$=d, $Y_5$=e)에 해당되는 케이스의 가중치를 모두 더해서 전체 가중치의 합으로 나누어주면 됩니다. 이렇게 모수 추정값을 다시 계산해 주는 것은 M-step 이라고 하는데 이렇게 업데이트된 모수 추정값은 식 (1)과 (2) 등에서 조건부 확률을 업데이트하는데 사용될수 있습니다. 이렇게 E-step 과 M-step 을 반복적으로 업데이트 하다보면 결국 수렴하는데 그 최종값으로 확률을 계산하면 그것을 EM 알고리즘이라고 합니다.

토론

1. 위의 사례는 EM 알고리즘의 가장 간단한 형태입니다. 다변량 범주형 자료에서 무응답이 있는 경우에는 위의 EM 알고리즘이 매우 간편하게 사용될수 있습니다. 그 외에도 EM 알고리즘은 엄청나게 많은 통계 모형 추정 및 예측 문제에 적용될수 있습니다. 앞으로 이와 관련된 내용을 좀더 다루려고 합니다.

2. 위의 예제에서는 원자료가 다항분포를 따르는 것으로 가정한 것입니다. 만약 0/1 의 binary 반응변수가 아닌 5점 척도의 ordinal data 라고 한다면 잠재변수에 대한 모형은 다항분포를 사용하는 것보다는 ordinal data 에 적합한 모형을 사용하는 것이 좋습니다. 예를 들면 cumulative logit model 같은 것이 그런 종류의 데이터에 더 적합한 모형입니다.

3. 위의 사례에서 가상 데이터를 만드는 방법론은 일종의 fractional imputation 이라는 방법으로 이해할 수도 있습니다. 제가 2011년 Biometrika 라는 저널에 발표한 논문 “Parametric fractional imputation for missing data analysis”에도 소개된 방법론입니다. SAS 14.1 의 Proc Surveyimpute 라는 프로시져의 FEFI 옵션이 바로 이 방법론을 구현하는 것입니다. 제가 쓴 책에도 자세히 나와 있습니다. (링크 참조: https://www.crcpress.com/Statistical-Methods-for-Handling-Incomplete-Data/Kim-Shao/p/book/9781439849637)

. .

. .

이번주 블로그 포스팅은 아래의 기고글로 대신합니다.

http://www.dongascience.com/news/view/14553

예제

1936년 미국에서는 민주당의 루즈벨트 대통령과 공화당의 랜든 후보가 대선 후보로 지명되어 열띤 선거전을 했었습니다. 그 당시 미국은 대공황을 겪었고 루즈벨트 대통령은 뉴딜 정책으로 경제를 살려보려고 했었고 경기가 조금씩 회복하는 분위기여서 대부분의 신문사들은 루즈벨트 대통령의 재선을 낙관했었다고 합니다.

그러나 그 당시 Literacy Digest 라는 잡지사에서는 루즈벨트 대통령이 전체 유권자의 43% 의 득표밖에 받지 못해서 재선에 실패할 것 이라고 예측치를 내어 놓았는데 이 예측치는 대략 천만명으로 구성된 샘플에게 우편 설문지를 보내 그 중에서 응답한 2백 40만명을 바탕으로 해서 얻어진 값이었습니다. (이 천만명의 주소는 전화번호부와 여러 클럽 회원 명부에서 얻었다고 합니다) 그 당시 Literacy Digest 지는 1916년부터 대선 예측을 다 맞추었고 240만명이라는 엄청난 크기의 자료를 통해 (아마 최초의 빅데이터?) 예측을 했기에 이 예측 결과는 정치계에 큰 충격을 주었다고 합니다.

이때 조지 갤럽이라는 젊은 여론 조사 학자는 그당시 2백 40만명의 응답을 취합하는데 몇달이 걸리므로 자신은 이중 일부만 샘플링을 통해서 미리 추정을 하고 (44%가 될것이라고 예측했음) 그와는 별도로 5만명의 랜덤 샘플을 통해 루즈벨트가 56% 을 득표할 것이라고 예측하여 루즈벨트 대통령이 재선에 성공할 것을 맞추었습니다. (실제 루즈벨트 대통령의 대선 특표율은 62%였다고 합니다.) 240만명을 바탕으로 얻어진 Literacy Digest 예측이 단지 5만명을 바탕으로 얻어진 갤럽의 예측보다 훨씬 부정확했던 것입니다. 왜 그런 일이 생겼을까요?

풀이

위의 사건은 여론조사 분야에서 매우 유명한 사건입니다. 모집단이 ${y_1, \cdots, y_N}$으로 구성되었다고 하고 이중 일부를 표본으로 뽑아서 이 표본 평균으로 추정한다고 할때 모수는 $\theta=N^{-1} \sum_{i=1}^N y_i=\bar{Y}$으로 표현될수 있고 표본 평균은 $\hat{\theta}= (\sum_{i=1}^N \delta_i y_i )/ (\sum_{i=1}^N \delta_i )$으로 표현될수 있는데 이때 $\delta_i$는 원소 $i$가 표본에 뽑히면 1 아니면 0 을 갖는 지시변수입니다.

이때 추정량 $\hat{\theta}$의 모수 $\theta$에 대한 편향은 $Bias(\hat{\theta})= E( \hat{\theta})- \theta$으로 정의될수 있는데 $\delta_i \sim Bernoulli(\pi_i)$ 이라 가정하면 \begin{eqnarray} Bias( \hat{\theta}) &\doteq & \frac{\sum_{i=1}^N \pi_i (y_i -\theta) }{\sum_{i=1}^N \pi_i}
&= & Corr( \pi, Y) \times CV(\pi) \times CV(Y) \times \theta \end{eqnarray} 으로 계산되는데 여기서 $CV(\pi)= \frac{\sqrt{V(\pi)}}{ E(\pi)} = \frac{ \sqrt{ \sum_{i=1}^N ( \pi_i - \bar{\pi})^2/N} }{ \bar{\pi} }$ 이고 (여기서 $\bar{\pi}=\sum_{i=1}^N \pi_i/N$) $CV(Y)= \frac{\sqrt{V(Y)}}{ E(Y)} = \frac{ \sqrt{ \sum_{i=1}^N ( y_i - \bar{Y})^2/N} }{ \bar{Y} }$ 이다. 따라서 $\pi_i$들의 상대 변동이 작거나 (표본 포함 확률이 거의 비슷하거나), Y 들의 상대 변동이 작거나 (Y들이 대체적으로 동질적이거나) 하면 표본 평균을 사용한 추정이 편향이 작아지게 됩니다. 그러나 이보다 더 중요한 것은 $Corr( \pi, Y) = \frac{\sum_{i=1}^N (\pi_i - \bar{\pi}) (y_i - \bar{Y})}{\sqrt{ \sum_{i=1}^N ( \pi_i - \bar{\pi})^2} \sqrt{ \sum_{i=1}^N ( y_i - \bar{Y})^2}}$ 으로 정의되는 $\pi$와 $Y$ 의 상관 계수로써 이 상관계수가 0 이 아니면 편향이 발생하게 되는 것이다. 즉, 표본의 대표성 문제 (표본 평균의 편향 문제)는 결국 추출 확률이 관심변수와 상관관계가 있다는 데에서 연유합니다.

위의 예제에서는 표본이 그 당시 집에 전화가 있는 사람들 위주로 추출되었기에 부자들이 체계적으로 더 많이 뽑혔고 전통적으로 부자들은 공화당을 지지하므로 위의 correlation 값이 음수를 갖게 됩니다. 따라서 Bias 역시 음의 값을 가지는 결과를 가져온 것입니다.

토론

1. 모집단의 개체들이 표본으로 뽑힐 확률이 체계적으로 다른 경우 발생되는 추정량의 편향 문제를 선택 편향 (selection bias)라고 합니다. 확률 표본 설계를 통해 얻어지는 자료가 아닌 자발적 참여를 통해 얻어지는 자료의 경우에는 이러한 selection bias 가 큰 문제가 될수 있습니다. 예를 들어 어떤 교육 프로그램에 참여하는 사람과 그렇지 않은 사람들 간의 평균 성취도의 차이가 교육 프로그램 자체의 효과로 인한 것인지 아니면 프로그램 참여라는 선택에 미치는 잠재 변수의 차이로 인한 것인지 알아낼 수가 없기 때문입니다. 그래서 이런 selection bias 를 줄여주기 위해 선택에 대한 모형을 이용하여 보정해주는 연구가 계량경제학을 위시로 하여 (Heckman 모형) 아직도 활발히 연구되는 분야입니다.

2. 위의 Bias 공식에서 중요한 포인트는 Bias 가 표본수에 의존하지 않는다는 것입니다. 즉, 아무리 샘플수를 늘려도 Bias 는 줄어들지 않는다는 것입니다. 이 Bias 때문에 선택편향이 있는 자료에서는 대수의 법칙(law of large numbers)이 성립하지 않게 되는 것입니다. 그래서 Literacy Digest 지의 추정은 240만명의 샘플을 사용했지만 단지 5만명의 샘플을 사용한 갤럽의 조사보다도 더 부정확한 결과를 가져온 것입니다. 일반적으로 추정량의 정확성은 오차제곱 평균 (Mean squared error; MSE)이 라는 것으로 평가될수 있는데 $MSE(\hat{\theta}) = Var( \hat{\theta}) + \{ Bias( \hat{\theta}) \}^2$ 으로 표현됩니다. 여기서 분산 $Var( \hat{\theta})$은 표본의 크기가 클수록 그 값이 줄어들게 되지만 ${ Bias( \hat{\theta}) }^2$은 표본수와는 관계없고 $r=Corr( \pi, Y)$에 관계 있으므로 Literacy Digest 의 표본처럼 선택 편향이 있는 자료에서는 결국 $Bias$의 제곱값이 줄어들지 않고 남아있어서 $r=0$이 되도록 샘플링을 한 갤럽의 예측치보다 (훨씬 많은 표본수에도 불구하고) MSE 가 더 커지게 되는 결과를 가져온 것입니다.

3. 위의 표본 포함 확률 $\pi_i$를 아는 확률 표본 추출의 경우에는 표본 평균을 사용하지 않고 $\pi_i$의 역수를 가중치로 하는 추정량을 사용하면 추정량의 편향이 사라지게 됩니다. 이를 Horvitz-Thompson 추정량이라고 하는데 표본 조사론 분야에서는 많이 사용되는 방법입니다. 자발적 선택으로 표본이 얻어진 경우에는 $\pi_i$를 알지는 못하고 이에 대한 모형을 써서 그 추정치를 가중치로 사용할수 있습니다. 이를 Propensity score weighting 이라고도 합니다. . .

. .

예제

매년 각국 학생들의 수학 및 과학, 독해력 등의 성취도를 측정하는 PISA (Programme for International Student Assessment) 프로그램에서는 각 나라별로 수천명 이상의 학생을 랜덤하게 뽑아서 표본으로 뽑힌 학생들에게 표준화된 시험을 보게 하여 학생들의 평균 학업 능력을 측정하여 발표합니다. 이런 PISA 프로그램을 한국에서 실시하고자 할때 학생 표본을 10,000 명을 뽑기로 하였다고 합니다. 이때 20개의 학교를 랜덤하게 뽑고 그 학교에서 학생을 500 명씩 뽑는 방법 (A 방법) 과 10 개의 학교를 랜덤하게 뽑고 그 학교에서 학생을 1,000 명씩 뽑는 방법 (B 방법)을 고려하고 있다고 합시다. 두 방법 모두 10,000 명을 뽑는 것인데 통계학 적으로 무슨 차이가 있을까요? 만약 차이가 있다면 어떤 방법이 더 좋은 것일까요?

풀이

표본수는 추정의 정확성에 영향을 미칩니다. 추정의 정확성을 신뢰구간의 길이로 표현한다고 하면 정확성은 표본수의 제곱근에 비례한다고 이야기 할수 있습니다. 그러나 이때 말하는 표본수라는 것은 엄밀하게 말하면 simple random sampling (SRS)으로 얻어진 자료로 가정했을때의 표본수, 즉 유효표본수(effective sample size)를 지칭하는 것입니다. 위의 예에서는 일단 표본 학교를 먼저 뽑고 그 학교에서부터 학생을 뽑는 집락 추출 (cluster sampling)을 사용했으므로 표본 자료는 집락내에서 양의 상관관계를 가집니다. 즉, $y_{ij}$를 학교 $i$ 내의 $j$-번째 학생의 성적이라고 하면 $y_{ij} = \mu+ a_i + e_{ij}$으로 표현되어서 $\mu$는 전체 평균, $a_i\sim (0, \sigma_a^2)$ 는 학교로 인한 차이, $e_{ij}\sim (0, \sigma_e^2)$는 개인별 차이로 설명될수 있습니다. 따라서 $\rho \equiv Corr (y_{ij}, y_{ij'}) = \frac{ \sigma_a^2}{ \sigma_e^2 + \sigma_a^2} , \ \ \ \mbox{ for } j \neq j'$ 으로써 양의 상관관계를 갖게 됩니다. 만약 학교별 평균 성적 차이가 크면 $\sigma_a^2$ 값은 커지게 되고 그 경우 $\rho$의 값도 커지게 됩니다.

문제는 이러한 경우 추정치의 정확성을 저해한다는 것입니다. 왜냐하면 위의 모형 하에서 표본평균 $\bar{y}$ 의 분산은 대략 $V( \bar{y})= \left\{ 1+ ( m-1) \rho \right\} \sigma_y^2 /N$ 으로 표현되는데 여기서 $N$은 총 표본수이고 $m$는 학교 내의 표본 학생수가 됩니다. 따라서 SRS 를 했을때의 분산값인 $\sigma_y^2/N$ 보다 $\{ 1+ ( m-1) \rho \}$만큼 분산이 증가하게 됩니다. 이때 이 증가분 ${ 1+ ( m-1) \rho }$은 집락 추출을 사용함으로써 발생되는 분산의 증가분을 나타내는 것으로써 표본 설계 효과 (design effect) 라고도 불립니다. 따라서 학교내에 표본 학생수가 적을수록 design effect 는 작아지게 되므로 B 방법 보다는 A 방법이 추정량의 더 분산이 작게 나오는 방법으로 표본수 자체는 같지만 더 효율적이 되는 것입니다.

토론

1. 위의 예에서 $\rho$는 집락내 상관계수 (intracluster correlation)이라고 불리우는데 같은 집락 내의 동질성의 척도입니다. 집락내 상관 계수가 낮더라도 집락 크기($m$)가 크면 표본 설계 효과는 크게 되고 따라서 유효표본수는 작아집니다. 예를 들어 $\rho = 0.1$ 이라고 한다면 $m=1,000$인 경우에는 표본 설계 효과는 $1+ (1000-1)*0.1 \doteq 100$ 이 됩니다. 실제 표본수가 10,000 이라도 유효 표본수는 실제 표본수를 표본 설계 효과로 나누어 주기에 100 밖에 되지 않는 것입니다. 이러한 자료로부터 $t$-test 를 하고자 할때 유효 표본수를 사용하지 않고 그냥 실제 표본수를 사용하여 분석을 하면 분산을 너무 작게 추정하는 오류를 범하게 되는 것입니다.

2. 빅데이터를 마치 표본 자료처럼 분석하고자 할때에도 위와 같은 오류를 범하기 쉽습니다. 자료들간의 상관관계를 무시하고 마치 독립된 자료인 것으로 간주하고 분석을 하면 분산을 과소추정하여 치명적인 오류인 일종오류를 범하게 됩니다. 자료의 구조가 통계 모형에 제대로 반영되어야 유효표본수를 계산해 낼수 있고 그로부터 통계학적으로 근거를 가지는 분석을 해낼수 있는 것입니다. 사이즈가 큰 자료를 가지고 통계 분석하고자 할때에는 반드시 체크해야 하는 사항입니다.

. .

. .

예제

어느 국가에서 갑자기 발생한 전염병의 감염율을 조사하기 위하셔 표본을 뽑고자 한다. 표본 추출은 단순 임의 추출로 한다고 할때 표본으로부터 얻어진 감염율이 참값과의 차이가 95% 의 확률 내에서 d=0.01 이내가 되기 위해서는 대략 몇개의 표본을 뽑아야 할 것인가?

풀이

어떤 특정 모집단에서 얻어진 $n$ 개의 자료를 $X_1, \cdots, X_n$ 이라고 할때 이 자료로부터 얻어지는 표본 평균 $\bar{X}_n$은 근사적으로 정규분포를 따릅니다. 이것이 그 유명한 중심극한 정리(central limit theorem)입니다. 본 예제에서는 $X_i$ 는 감염률이 $p$ 인 베르누이 분포에서 $n$ 개를 관측한 것으로 볼수 있으므로 평균이 $E(\bar{X}_n)=p$ 이고 분산이 $V(\bar{X}_n)=p(1-p)/n$ 인 정규분포를 따르게 될 것입니다.

따라서 $P \left( \left| \frac{ \bar{X}_n -p }{\sqrt{p(1-p)/n} } \right| <1.96 \right)= 0.95$ 을 이용하여서 $n$ 값을 구하면 $\sqrt{p(1-p)/n} \times 1.96 < d =0.01$ 으로써 즉 $n > \left( \frac{ 1.96}{d} \right)^2 p (1-p)$ 이 되는데 우항의 최대값이 $p=1/2$에서 위치하므로 대략적으로 $n > d^{-2}$으로 계산되어 $d=0.01$인 경우에는 $n=10,000$ 명의 표본이 필요하게 되는 것입니다.

토론

1. 어떤 특정 모집단에서 얻어진 $n$ 개의 자료를 $X_1, \cdots, X_n$ 이라고 할때 이 자료로부터 얻어지는 특정 통계량 $Y=u(X_1, \cdots, X_n)$과 관련된 통계적 성질에 더 관심이 많을 경우가 있습니다. 이때 얻어지는 통계량 $Y$의 확률 분포를 표본 분포(sampling distribution)이라고 부릅니다. 즉, 통계량 $Y$의 표본 분포는 자료 $(X_1, \cdots, X_n)$의 확률분포로 부터 유도되는 것입니다.

2. 중심 극한 정리는 자료의 확률 분포에 대한 법칙이 아니라 통계량의 표본 분포에 대한 법칙입니다. 표본 수가 증가한다고 모집단의 확률 분포는 변하지 않습니다. 적절한 조건 하에서 자료에서 얻어지는 통계량의 표본 분포가 정규 분포로 근사될수 있다는 것이 표본 분포(sampling distribution)의 정규 근사(normal approximation)의 핵심 내용입니다.

3. 자료의 확률 분포를 아는 경우 통계량의 표본 분포를 찾는 방법은 몬테 카를로 시뮬레이션을 통해서 이해할수 있습니다. 주어진 확률 분포로부터 자료를 발생시킨후 통계량을 발생시키는 것을 반복하면 통계량에 대한 히스토그램을 그릴수 있는데 이 히스토그램이 그 통계량의 표본 분포에 대한 근사(approximation)가 되는 것입니다. 문제는 자료의 확률 분포를 정확히 모른다는 것이지요. 위의 예에서 베르누이 분포의 모수값 p 를 모르기에 몬테 카를로 시뮬레이션으로 자료를 발생시키지 못합니다. 하지만 표본 평균으로 p 를 추정하여 자료의 확률 분포가 Bernoulli $(\hat{p})$ 인 것으로 간주하여 시뮬레이션을 하는 것을 생각할 수 있는데 이것이 바로 parametric bootstrap 의 개념입니다. 즉, 중심극한 정리를 이용하지 않고 통계량의 표본 분포를 계산하는 방법으로 parametric bootstrap 을 이용할 수 있습니다.

4. 위의 예에서 d 는 보통 신뢰구간의 길이로써 오차 한계(margin of errors) 라고도 부릅니다. 여론 조사에서 오차 한계가 3% 이다 라고 이야기 하는 것은 대부분 표본수가 $(0.03)^{-2} \doteq 1,000$ 이라는 이야기를 좀더 포장하여 표현하는 것입니다. 실제로는 여론 조사 자체에서 발생하는 여러가지 편향 때문에 실제 오차는 종종 이보다 더 크게 나타납니다.

. .

예제

어느 시골 마을 입구 큰길 앞에서는 택시가 1시간당 10대가 지나간다고 합니다. 그 10대중 2대만 빈차라고 하고 합승을 해주지 않는다고 합니다. 당신이 택시 타는 줄의 두번째에 있다면 당신이 택시를 잡기 위해 1시간 이상 기다리게 될 확률은 얼마일까요?

풀이

택시가 오는 시각이 서로 독립적인 경우 한 시간동안 지나가는 택시의 갯수는 포아송 분포를 따릅니다. 위의 예에서는 한 시간 동안 오는 빈 택시의 갯수는 평균이 2인 포아송 분포, 즉 Poisson(2),가 되는 것입니다. 그런데 이를 t 시간으로 확장하면 t 시간동안 오는 빈택시의 빈도수 $N_t$ 는 t 에 따라 증가하는 확률 변수로써 그 분포가 Poisson (2t) 를 따르고 이를 Poisson process 라고 합니다.

포아송 프로세스는 사건 발생 횟수에 대한 확률 변수라고 한다면 거꾸로 사건이 발생하는 시간에 대한 확률 변수를 생각할수 있습니다. 즉, $T_1$을 첫번째 택시가 오는데 걸리는 시간, $T_2$를 두번째 택시가 오는데 걸리는 시간, 이런 식으로 정의하면 이 문제는 $P( T_2 >1)$을 계산하고자 하는 문제입니다. 그런데 포아송 프로세스의 정의상 ${ T_2 >t }$ 이라는 사건은 ${N_t \le 1}$이라는 사건과 동등하므로 (왜냐하면 두번째 택시가 t 시간 이후에 온다는 이야기는 t 시간 까지 빈택시 횟수가 1개 이하라는 이야기이므로) 우리가 원하는 확률은 $P( T_2 >1) = P ( N_1 \le 1) = P(N_1 =0) + P(N_1=1) = e^{-2} (1+2)=3e^{-2}=0.41$ 으로 계산됩니다. 참고로 $X$ 가 Poisson(m) 를 따른다는 것은 $P(X=k)= \frac{ m^k}{k!} e^{-m}$ 임을 의미합니다.

토론

1. 포아송 프로세스 $N_t$는 처음에는 0부터 시작해서 $t$가 늘어남에 따라 하나씩 증가합니다. (단위 시간당 발생 비율은 $\lambda$ 값을 가집니다.) 이 포아송 프로세스가 처음으로 1 의 값을 갖게 되는 사건이 발생하는 시간 (제 1회 사건 발생시간)의 확률분포가 바로 지수분포가 됩니다. 왜냐하면 ${N_t=0}$ 이라는 사건은 시간 $t$까지는 관심 사건이 아직 일어나지 않았다는 이야기이고 이는 처음 발생은 t 이후에 일어난 다는 의미이므로 ${T_1>t}$과 동등합니다. 따라서 $P( T_1>t)=P(N_t=0) = \frac{(\lambda t)^0}{0 !}\exp (-\lambda t)$ 으로 되므로 $T_1$ 의 누적 분포 함수 (cumulative distribution function)은 $1-\exp (-\lambda t)$이 되어서 평균이 $1/\lambda$ 인 지수 분포를 따르는 것을 확인할수 있습니다.

2. 만약 첫번째 사건이 발생하는 시간에 관심이 있는 것이 아니라 첫 $r(>1)$ 번째 사건이 발생하는 시간에 관심이 있는 경우에는 서로 독립인 지수분포가 $r$개 합해진 형태인 감마 분포를 따르게 됩니다. 앞의 정의를 이용하면 $T_r$은 첫 $r$-번째 사건이 발생하기까지 걸리는 시간이 되는데 ${T_r > t}$라는 것은 그 소요 시간이 $t$보다 크다는 것이고 따라서 시간 $t$에는 아직 $r$-번째 사건이 발생하지 않았다는 의미이므로 ${N_t \le r-1}$과 동일해 집니다. 따라서 $P (Y_r \le t ) = 1- P ( Y_r > t) = 1 - P ( N_t \le r-1 ) = 1 - \sum_{k=0}^{r-1} \frac{ (\lambda t)^k}{ k!} e^{-\lambda t } = \sum_{k=r}^{\infty} \frac{ (\lambda t)^k}{ k!} e^{-\lambda t }$
   이 얻어지고 이를 $t$에 대해 미분하면 $f_Y( t ) = \frac{\lambda^k }{(k-1)! } t^{k-1}e^{- \lambda t }$
   이 되어서 감마 분포의 확률 밀도 함수와 동일해 집니다. 위의 택시 예에서 두번째 줄에 있으므로 택시를 잡는 사건이 두번째 일어나는 시간에 관심이 있는 것입니다. 그 시간이 감마 분포를 따르게 되는 것입니다.

3. 지수 분포의 특징은 memoryless property 라는 것인데 이를 수식으로 표현하면 다음과 같습니다. $P(T_1 > t+s \mid T_1 >t) = P(T_1>s).$ 즉, 그 전에 얼마 동안 기다렸다는 것이 아무런 상관이 없다는 것입니다. 이는 결국 포아송 프로세스의 정의에서 연유하는 것입니다. 어떤 사건이 발생하는 사건(arrival time)이 서로 독립적으로 분포하기 때문에 이런 성질을 갖는 것입니다.

. .

예제

몬테 카를로 시뮬레이션이란 자료가 생성되는 자료 생성 메카니즘을 고정시켜 놓고 난수 발생을 통해 반복적으로 자료를 발생시켜서 원하는 분석을 반복적으로 실시함으로써 그 분석 방법에 대한 과학적 결론을 얻어내고자 하는 계산 방법론입니다. 특히 통계학에서는 통계량을 계산하는데 이러한 시뮬레이션을 통해 통계량의 분포를 얻어내는 것이 이론적으로 구하는 것보다 쉬운 경우가 종종 있어서 많이 사용됩니다. 홍길동씨는 특정 통계량을 이런 몬테 카를로 실험을 통해서 (1,000회 반복했다고 합시다) 얻어낸 후에 그것들의 평균과 분산을 구했더니 각각 1.021646 이라는 값과 0.00564345 이라는 값을 얻었다고 합시다. 이럴때 최종 결과를 1.025646 이라고 발표하는 것이 어떤 문제가 있을까요?

풀이

몬테 카를로 시뮬레이션은 통계학 뿐만 아니라 물리학이나 생물학과 같은 자연과학과 공학 및 사회과학에서도 많이 사용되는 방법입니다. 한번의 시뮬레이션을 돌리면 얻어지는 통계량은 시그널(signal) 부분과 노이즈(noise) 부분이 합쳐져서 얻이지지만 대수의 법칙에 의해서 이를 무한히 반복하여 평균을 내면 노이즈 부분은 사라지고 시그널 부분만 남습니다.

그러나 위의 예처럼 1,000번의 시뮬레이션을 하더라도 노이즈 부분은 완전히 사라지지 않고 조금 남습니다. 위의 예에서 1.021646 이라는 값은 시그널과 노이즈가 결합된 형태에서의 1000 개의 통계량의 평균이므로 시그널 파트는 동일하지만 노이즈 파트는 시뮬레이션 횟수가 증가함에 따라 달라질수 있기 때문입니다. (다만 시뮬레이션 횟수가 많으면 노이즈의 평균은 0 으로 수렴합니다. 이게 대수의 법칙입니다)

따라서 시뮬레이션의 결과를 report 할 경우에는 시그널 부분만을 보여주는 것이 좋습니다. 왜냐하면 그렇게 하는 것이 재현성(reproducibility)의 측면에서 볼때 바람직 하기 때문입니다. 동일한 시나리오에서 얻어지는 시뮬레이션 결과가 사람마다 달라질텐데 노이즈 까지 포함해서 발표를 하는 것은 독자를 혼란하게 하고 발표자의 과학자적 자질에 의심을 품게 만드는 것입니다.

그러면 어떻게 시그널을 찾을수 있을까요? 가장 쉬운 방법은 이 시뮬레이션 결과로부터 모수의 신뢰구간을 구하면 됩니다. 즉, 시뮬레이션 자료가 표본수가 1,000 인 표본 자료로 보고 표본 평균이 1.021646 이고 표준 편차는 $\sqrt{0.00564345}=0.0751229$ 이므로 95\% 신뢰구간은 $(1.021646 - 2*0.0751229/\sqrt{1000}, 1.021646 + 2*0.0751229/\sqrt{1000})=(1.016895, 1.026397)$ 이 됩니다. 즉, 소수점 두자리 숫자 자체가 이미 정확한 값이 아니라 노이즈를 반영하는 값이라는 것입니다. 따라서 이 시뮬레이션 결과를 발표할때에는 유효수자 소수점 첫째자리까지인 1.0 으로 발표하는 것이 맞는 방법입니다.

토론

1. 통계학 초기 역사에서 가장 유명한 몬테 카를로 시뮬레이션은 T-분포의 분포 테이블을 만든 것입니다. 기네스 맥주 회사에서 근무하는 고셋 경이 T-통계량의 유한 표본에서의 sampling distribution 을 구하고자 몬테 카를로 시뮬레이션을 했는데 그 당시에는 컴퓨터가 없어서 난수표를 이용해서 했다고 합니다. 이와 관련하여 제가 몇년 전에 작성한 블로그 링크를 추가합니다. (

http://blog.naver.com/kim00020/60058915894 )

1. 몬테 카를로 시뮬레이션의 이론적 근거는 대수의 법칙입니다. 우리가 원하는 모수가 어떤 사건이 일어날 확률이라고 한다면 그것은 적분으로 표현되는 것이지만 그 적분을 정확하기 계산하기가 어려우므로 컴퓨터 시뮬레이션으로 자료를 랜덤하게 발생시켜 그 중에서 특정 사건이 일어나는 빈도를 계산해서 근사적으로 계산하는 것입니다.

2. 몬테 카를로 시뮬레이션을 설계할 때 주의할 점으로는 두 추정량을 비교할때 동일한 샘플을 사용해야 한다는 것입니다. 즉, A 라는 방법론을 사용한 추정량 $\hat{\theta}_A$ 와 B 라는 방법론을 사용한 추정량 $\hat{\theta}_B$ 의 sampling distribution 을 몬테 카를로로 구현하고자 할때 두 방법론이 동일한 자료로 부터 계산되어야 더 효율적입니다. 즉, $V( \hat{\theta}_A - \hat{\theta}_B) = V( \hat{\theta}_A) + V( \hat{\theta}_B)- 2 Cov ( \hat{\theta}_A, \hat{\theta}_B)$ 에서 두 통계량이 서로 별도의 독립적인 자료로부터 계산되면 Covariance 가 0 이 되지만 같은 자료로부터 계산되면 대부분의 경우 covariance 가 양수가 나오므로 두 통계량의 비교에 더 효율적인 실험 설계법이 되는 것입니다.

예제

우생학의 창시자로 알려진 프랜시스 골턴(Francsis Galton)의 에피소드는 집단지성의 위력을 이야기할 때 종종 인용됩니다. 골턴은 1906년 우연히 우시장에서 황소 무게를 눈대중으로 알아맞히는 대회를 구경하는데 이 대회에는 무려 800여 명이 참가했는데 그 누구도 황소 무게를 정확히 맞히지 못했다고 합니다.

하지만 놀라운 것은 대회 참가자들이 적어낸 무게의 평균을 계산해본 결과 황소의 실제무게와 거의 비슷했다는 것이다. 추정 평균치는 1197파운드였고 황소의 실제 무게는 1198파운드였다고 합니다. 황소 전문가들이 써낸 추정치보다 전체 참가자 추정치의 평균값이 더 정확했던 것입니다.

골턴은 이 이해할 수 없는 놀라운 이야기를 1907년 과학저널 ‘네이처’에 발표합니다. 개별적인 한 사람 한 사람의 지식이나 지혜는 미미해보일지라도 그것이 모이면 누구도 예상 못한 놀라운 힘을 발휘할 수 있음을 시사해주는 것입니다.

풀이

각 개인의 황소 몸무게 예측치를 $Y_i$ 라고 표현하고 황소의 몸무게 참값을 $\mu$ 라고 한다면 각 개인 추정치의 오차는 $e_i = y_i - \mu$ 로 표현될 것입니다. 어떤 사람은 황소 몸무게를 지나치게 크게 추정하여 $e_i$ 값이 양의 값을 갖게 되고 또 어떤 사람은 황소 몸무게를 작게 추정하여 $e_i$ 값이 음의 값을 갖게 되겠지만 평균적으로는 $E(e_i)=0$ 이라고 볼수 있을 것입니다. (즉, 사람들의 예측에 편향이 없다고 볼수 있을 것입니다. ) 그러한 경우에는 표본 평균이 참값에 확률적으로 수렴합니다. 즉, $\bar{Y}n=n^{-1} \sum{i=1}^n Y_i$과 참값 $\mu$ 과의 차이를 표본 평균의 오차라고 한다면 이 오차는 표본수가 커질수록 확률적으로 0 에 수렴합니다. 이를 대수의 법칙이라고 합니다. 이러한 대수의 법칙을 수식으로 표현하면 다음과 같습니다. $\lim_{n \rightarrow \infty} P \left\{ \left| \bar{Y}_n - \mu \right| > \epsilon \right\} = 0 , \ \ \mbox{for any }\epsilon >0.$

이를 좀더 알아보기 위하여 0.5 의 확률로 1 또는 0 의 값을 갖는 베르누이 확률변수 Y 를 n=1000 개 독립적으로 발생시킨후 각 $t=1,2, \cdots, n$ 에 대해서 $\bar{Y}t=\sum{i=1}^t y_i / t$ 를 계산한후 이를 $Y$ 축으로 놓고 $t$ 값을 $X$ 축에 놓아 그래프를 그려보면 다음과 같이 나옵니다.

처음 표본수가 낮은 경우에는 어느 정도 변동이 있지만 표본수 $n$이 증가하면서 표본평균 값이 0.5 를 중심으로 크게 변하지 않음을 알수 있습니다.

토론

1. 대수의 법칙은 체비세프 부등식을 이용해서 쉽게 구할수 있습니다. 체비세프 부등식을 사용하면 $P \left\{ \left| \bar{Y}_n - \mu \right| > \epsilon \right\} \le \frac{E\{ \left| \bar{Y}_n - \mu \right|^2 \} }{\epsilon^2}$ 이 얻어지는데 $E(Y_i)=\mu$ 이 성립하면 위 식의 오른쪽 항은 $Var ( \bar{Y}_n)/ \epsilon^2$ 이 되고 이 분산은 $Y_i$ 들이 서로 독립인 경우에 분산이 $\sigma^2/n$이 되어 따라서 $n$이 증가함에 따라 0으로 수렴하게 되므로 위 식의 오른쪽도 0 으로 수렴하게 되는 것입니다.

2. 대수의 법칙을 처음 발견한 사람은 야곱 베르누이 입니다. 그는 체비세프 보다 100년 넘게 앞서 살았던 사람이었으므로 체비세프 부등식을 이용한 증명을 하지 못하고 $Y_i$ 가 지시변수인 경우 $\sum_{i=1}^n Y_i$ 의 분포가 이항분포를 따른다는 것을 이용하여 이를 바탕으로 복잡하게 증명을 완성했습니다. 그래서 그의 이름을 따서 $Y$ 가 1 또는 0의 값을 갖는 확률 분포를 베르누이 분포라고 부릅니다.

3. 대수의 법칙이 가장 많이 사용되는 분야는 아마도 몬테 카를로 계산일 것입니다. 몬테 카를로 계산법은 원하는 모수(parameter)를 $\mu = E(X)$ 의 형태로 표현한후 시뮬레이션을 통해 $X_i$ 들을 아주 많이 발생시킨 후 이들의 표본 평균으로 근사하는 방법입니다. 조건부 확률같은 확률 역시 일종의 평균 이므로 몬테 카를로 시뮬레이션을 통해서 근사적으로 구현할수 있습니다. 이 몬테 카를로 계산은 컴퓨터의 발달과 함께 통계학의 적용 범위를 급격히 넓혀준 유용한 기법이며 그것의 이론적 근거는 대수의 법칙입니다.

4. 대수의 법칙과 관련하여 흥미로운 것은 표본수가 2배 증가한다고 해도 정확도가 2배 증가하는 것이 아니라 2의 제곱근인 1.414 배만큼 증가한다는 것입니다. 표본수가 $n$에서 $2n$ 이 되면 표본평균의 표준 오차가 $\sigma/\sqrt{n}$ 에서 $\sigma/\sqrt{2n}$ 으로 $\sqrt{2}$ 만큼 줄어듭니다. 따라서 비용은 표본수에 비례해서 증가하지만 그에 따른 효용은 그의 제곱근에 비례해서 증가하게 되므로 적절한 수준에서 비용과 효용을 고려한 표본수를 결정할 필요가 있습니다.

5. 아래 그래프는 두개의 베르누이 변수로부터 각각 표본 평균을 구한후 그것들의 오차가 얼마나 빨리 0으로 수렴하는가를 보여주는 그래프입니다. $Y$는 Bernoulli (0.5) 분포로 부터 발생했지만 $Z$는 Bernoulli (0.1) 분포로부터 발생했습니다. 이 그래프로부터 확인할수 있듯이 $Y$의 평균보다는 $Z$의 평균이 더 빨리 모평균에 가깝게 수렴합니다. 그 이유는 $V(Y)=0.5 0.5=0.25$ 이지만 $V(Z)=0.90.1=0.09$ 이므로 $\bar{Y}_n$ 의 분산이 $\bar{Z}_n$의 분산보다 크기 때문입니다.

예제

컴퓨터용 팬(fan)을 만드는 회사에서 제품이 고장나면 무상으로 교체해주는 워런티 기간을 3년으로 잡았다고 합시다. 컴퓨터용 팬과 같은 기계적 제품의 생존 시간을 X 라고 할때 그 확률 값은 $P( X < x ) = 1- \exp \left( - \frac{x}{\theta}\right), \ \ x>0$ 을 따르는 것이 알려져 있고 여기서 모수 $\theta=E(X)$값은 회사 제품별로 다른데 그 회사의 과거 자료를 분석한 결과 그 제품의 평균 수명은 5만 시간이라는 것이 밝혀졌다고 합니다. 그렇다면이 3년 워런티에 따라 발생하는 추가 비용은 제품값의 몇 \% 으로 계산해야 할까요?

풀이

3년은 최대 365 day $\times$ 24 hours $\times$ 3 years = 26,280 hours 로 계산될수 있으므로 3년 내에 고장이 날 확률은 $P(X<26,280)$을 계산하는 문제로 귀결됩니다. 따라서 위의 식에 $\theta=50,000$을 대입하면 $P(X<26,280)= 1- \exp (-\frac{26280}{50000}) =0.4088$ 으로 계산되어 40.88\% 의 확률로 3년내에 고장이 발생할 것입니다. 따라서 원 제품값에 추가로 40.88\% 에 해당되는 비용을 워런티의 값으로 계산해 주어야 손해가 발생하지 않습니다. 즉, 제품 한대값이 원래 1만원이라면 워런티를 포함한 판매가는 1만 4천 88원보다 크게 잡아야 합니다. 만약 1년 워런티로 한다면 추가비용은 16\%로 줄어들 것입니다.

토론

1. 우리가 어떤 확률 변수 $X$에 대해 특정 범위에서의 확률값을 계산하기 위해서는 확률 분포(probability distribution)를 알아야 합니다. 위의 예제에서는 평균이 $\theta$ 인 지수 분포(exponential distribution)를 사용했는데 생존시간에 대한 확률 분포로써 지수 분포가 아닌 다른 분포 (예: Weibull 분포)도 사용될수 있습니다.

2. 확률 분포란 확률 변수의 실현치들의 분포가 어떤 규칙을 따른다고 보고 그 법칙을 수리적으로 기술한 것입니다. 현실 세계에서 이를 정확하게 파악하는 것은 매우 어려운 일이고 대신 이를 좀더 단순한 방식으로 근사적으로 표현하고자 하는 것이 통계적 모형의 핵심입니다.

3. 즉 통계 모형이란 어떤 확률 변수에 대한 확률 분포를 좀더 이해하기 쉬운 형태로 근사적으로 표현하여 문제를 단순화하여 이해하고자 하는 것입니다. 통계모형은 확률분포를 범주화하여 몇개의 모수(parameter)로 표현해 줌으로써 확률 분포를 결정하는 문제를 모형 선택과 모수 추정의 문제로 단순화 시키게 됩니다. 모형 선택은 해당 분야의 특정 도메인 지식으로 알려져 있는 경우가 종종 있고 모수 추정은 수집된 데이터를 통해서 얻는 경우가 많이 있습니다.

4. 위의 예에서 흥미로운 것은 평균 수명이 5만시간 이지만 그의 반 정도에 해당되는 2만 6천 시간이내에 고장날 확률이 40%나 된다는 것입니다. 평균 수명의 반도 안되서 고장나는 사건이 발생하는 빈도를 결코 무시할수 없다는 것입니다. 그런 면에서는 평균이 주는 착시 효과라고도 부를수 있을 것입니다. 이를 제대로 계산하지 못해서 워런티 값을 낮게 책정하면 많이 팔고도 손해를 볼수 있게 됩니다.

예제

유방암을 진단하는 검사 방법중의 하나인 유방촬영검사(Mammography) 검진법은 의료 전문가들 사이에서도 효용성 논란이 일고 있다고 합니다. 왜 그런지 좀더 알아보기 위하여 다음과 같은 배경지식을 가정하도록 합시다. (이 값은 실제와는 조금 다를수 있습니다.)

1. 성인 여성중 1%가 유방암에 걸린다.
2. 유방촬영검사는 유방암이 있는 경우 80%를 발견한다. (따라서 유방암이 있는 여성의 20%를 놓친다.)
3. 유방촬영검사는 유방암이 없는 경우에도 10%를 유방암이 있는 것으로 판독한다 (따라서 유방암이 없는 여성의 90%를 정상으로 밝혀낸다.)

이러한 경우 어떤 여성이 이 검사에서 양성 결과가 나왔다면 이 환자가 진짜로 유방암에 걸렸을 확률은 얼마인가요?

풀이

이 문제는 베이즈 정리를 이용한 사후 확률을 계산하는 가장 대표적인 예 중의 하나입니다. 현실 세계에서는 검사결과와 질병유무가 완전히 일치하지 않기 때문에 우리는 검사결과로 부터 질병 유무를 알아내고 싶을 때 사후 확률을 계산하고자 합니다. 검사 결과와 질병 유무에 따른 각 경우의 분류는 다음과 같습니다.

검사결과	질병 유	질병 무
양성(positive)	진양성(true positive)	위양성(false positive)
음성(negative)	위음성 (false negative)	진음성 (true negative)

위의 유방암 예제에서는 1%만이 유방암 환자이므로 10,000 명 중에 100명이 유방암 환자이고 9,900 명이 정상입니다. 따라서 100명의 환자 중에서 검사에서 양성을 갖는 진양성 환자는 80명이 될 것이고 그렇지 않은 위음성 환자는 20명이 될 것입니다. 또한 9,900명의 정상인 중에서 10%인 990명은 검사가 양성을 가지는 위양성 환자가 되고 나머지 90%인 8,910명은 검사가 음성을 가지는 진음성 환자가 됩니다. 따라서 이를 도표로 정리하면 다음과 같습니다.

검사결과	질병 유	질병 무	합계
양성(positive)	80	990	1,070
음성(negative)	20	8,910	8,930

따라서 검사 결과가 양성이 나왔더라도 대부분은 1만명 중의 990명에 해당되는 위양성 그룹에서 나온 것입니다. 검사 결과가 양성인 사람 중에서 실제 질병이 있는 사람의 비율은 $80/1,070 \doteq 0.0748$, 즉, 7.48\% 에 불과한 것입니다.

토론

1. 질병 진단법의 정확성은 두가지 측면에서 계산됩니다. 하나는 민감도(sensitivity)라고 불리는 것으로써 질병이 있는 사람 중에서 얼마나 높은 비율로 양성을 찾아내느냐 하는 것입니다. 다른 하나는 특이도(specificity)라고 불리는 것으로써 질병이 없는 사람 중에서 얼마나 많이 음성 결과를 얻어 내느냐를 나타냅니다. 위의 유방촬영검사는 80%의 민감도와 90%의 특이도를 가진다고 할수 있을 것입니다.

2. 베이즈 정리는 증거를 바탕으로 사건이 일어날 확률을 사후적으로 계산하는 방법에 관한 정리입니다. 위의 예에서 검사 결과가 양성이라는 것은 하나의 증거입니다. 이 증거를 통해서 해당 여성의 유병 확률은 1%에서 7.48% 로 높아지게 되는 것입니다. 이를 수식으로 나타내면 다음과 같습니다. \begin{eqnarray} P(\mbox{Cancer } \mid \mbox{ Positve} ) &=& \frac{ P( \mbox{Cancer } \& \mbox{ Positve} ) }{ P( \mbox{ Positive} )}
   &=& \frac{ P( \mbox{Cancer} ) P( \mbox{ Positive} \mid \mbox{Cancer} )}{P( \mbox{Cancer}) P( \mbox{ Positive} \mid \mbox{Cancer}) + P( \mbox{Normal} ) P( \mbox{ Positive} \mid \mbox{ Normal}) }
   &=& \frac{ 0.01 \times 0.8}{ 0.01\times 0.8 + 0.99 \times 0.1 }=0.07476636 . \end{eqnarray}

3. 이렇게 상태가 두개인 경우 (질병 유무)에는 사후 확률을 Odd 를 이용해서 계산하는게 더 쉽습니다. 위의 조건부 확률 공식을 사용하면 질병 유무에 대한 사후 확률의 비는 다음과 같이 표현될수 있습니다. $\frac{ P(\mbox{ Cancer} \mid \mbox{ Positve} ) }{ P(\mbox{ Normal} \mid \mbox{ Positve} ) } = \frac{P( \mbox{ Positive} \mid \mbox{ Cancer} ) }{P( \mbox{ Positive} \mid \mbox{ Normal} ) } \times \frac{P( \mbox{Cancer} ) }{P(\mbox{Normal})}.$ 위의 식에서 등식 오른쪽의 두번째 항은 초기 확률의 비입니다. 거기에 “Test positive”라는 증거를 반영한 보정 승수(evidence adjustment factor) $\frac{P( \mbox{ Positive} \mid \mbox{ Cancer} ) }{P( \mbox{ Positive} \mid \mbox{ Normal} ) }$ 을 곱해서 사후 확률의 비를 계산하게 되는 것입니다. 위의 유방암 검진 사례를 이용하면 초기 유방암 발생 비율은 암:정상의 비율이 1:99 입니다. 그런데 테스트 결과가 양성인 경우의 증거 보정 승수는 (0.8/0.1)=8 로 계산됩니다. 따라서 테스트 양성을 관측한 후의 암:정상의 비율은 8:99 가 되는 것입니다. 따라서 이 경우 사후 확률은 8/(8+99)=0.07476636 으로 위의 결과와 동일해 지는 것입니다.

4. 만약 검사를 한번 더 했다면 어떻게 될까요? 즉, 검사를 두번 독립적으로 실시해서 둘다 양성이 나왔다면 어떻게 될까요? 그러면 증거 보정 승수는 8*8=64 으로 높아지고 따라서 사후 확률은 64/(64+99)=0.392638 으로 계산됩니다. 즉, 두번 다 양성 결과가 나왔다고 하더라도 그 사람이 실제 유방암 환자일 확률은 39\% 밖에 되지 않는다는 것입니다.

5. 이처럼 확률은 증거를 통해서 계속 업데이트 됩니다. 유방암 관련 초기 확률이 1\% 에서 일회의 양성 결과를 통해 7.48\% 라는 값으로 업데이트 되었고 그 확률은 다시 한번의 양성 결과를 통과하면 39\% 로 높아집니다. 즉, 1회 관측으로 얻어진 사후 확률이 제 2회 관측을 반영하는 사후 확률 계산에서는 사전 확률로 사용되는 것입니다.

예제

몬티홀(Monty Hall) 문제를 들어보셨나요? 몬티홀 문제는 미국의 어느 방송국에서 했던 “Let’s Make a Deal”이라는 쇼 프로그램에서 진행자 Monty Hall 씨의 이름을 따서 붙혀진 흥미로운 게임입니다. 이 게임은 다음과 같이 진행됩니다.

1. 방송국 스튜디오에는 3개의 문이 닫힌채로 보여집니다. 그 문 뒤에는 각각 2개의 염소와 1개의 멋진 스포츠카가 임의의 순서대로 배치되어 있습니다. 주인공이 하나의 문을 선택하게 되는데 이 스포츠카가 있는 문을 선택하면 그 차를 갖게 되는 쇼입니다.

2. 먼저 이 프로그램에서는 주인공이 나와서 문을 하나 고릅니다. (예를 들어 첫번째 문을 골랐다고 합시다)

3. 그러면 그 프로그램의 진행자인 몬티 홀씨는 어느 문 뒤에 스포츠카가 있는지를 알기에 나머지 두 문 중에서 염소가 있는 문을 열어 보여 줍니다. (만약 문 1에 자동차가 있으면 문2와 문3중에서 하나를 임의로 선택해서 보여줍니다.) 그리고는 주인공에게 묻습니다. “당신은 처음의 선택을 바꾸겠습니까?” 하고 말입니다. 주인공은 처음의 선택한 문을 고르거나 아니면 나머지 남은 문을 고르거나를 결정해야 합니다.

주인공은 처음의 선택을 바꾸는게 나을까요? 아니면 둘다 마찬가지 인가요?

풀이

이 문제는 조건부 확률의 개념을 설명하는데 흔히 사용되는 사례입니다. 아무런 정보가 없을때에는 자동차가 특정 문 뒤에 있을 확률은 각각 1/3입니다. 즉, C1 을 첫번째 문에 자동차를 놓는 사건, C2 를 두번째 문에 자동차를 놓는 사건, C3 를 세번째 문에 자동차를 놓는 사건으로 한다면 P(C1)=P(C2)=P(C3)=1/3 입니다.

그런데 여기에 새로운 정보가 들어갑니다. 몬티홀씨의 행동은 어느 문에 자동차가 있는지에 대한 정보를 제공합니다. 주인공이 처음에 문 1을 선택했을때, D1 을 몬티홀씨가 첫번째 문을 보여주는 사건, D2 을 몬티홀씨가 두번째 문을 보여주는 사건, D3을 몬티홀씨가 세번째 문을 보여주는 사건이라고 정의한다면 몬티홀씨의 행동은 참값에 따라 다르게 결정되므로 다음과 같은 조건부 확률을 갖게 됩니다. $P(D3 \mid C1)= 1/2, \ \ P(D3 \mid C2)= 1, \ \ P(D3 \mid C3)=0.$ 즉, 자동차가 첫번째 문에 있다면 문2 와 문3에 둘다 염소가 있으므로 몬티홀씨는 두 문중에 하나를 같은 확률로 뽑아 보여줍니다. 하지만 자동차가 두번째 문에 있다면 나머지 두 문중에 염소가 있는 문은 문3 뿐입니다. 그래서 $P(D3 | C2)= 1$이 되는 것입니다. 마찬가지로 생각하면 자동차가 세번째 문에 있으면 $P(D3 | C3)= 0$이 됩니다.

따라서 만약 몬티홀씨가 세번째 문을 열었다면 조건부 확률의 정의를 이용하면 $P(C1 | D3 )= \frac{P(D3 |C1) P(C1)}{ P(D3 |C1) P(C1) + P(D3 |C2) P(C2) + P(D3 |C3) P(C3) }=1/3$ 으로 계산됩니다. 하지만 $P(C2 | D3 )= \frac{P(D3 |C2) P(C2)}{ P(D3 |C1) P(C1) + P(D3 |C2) P(C2) + P(D3 |C3) P(C3) }=2/3$ 이므로 문2에 자동차가 있을 확률이 더 높아집니다. 만약 몬티홀씨가 두번째 문을 열었다면 마찬가지 방법으로 $P(C1 \mid D2)=1/3$이고 $P(C3 \mid D2)=2/3$임을 계산할수 있습니다. 따라서 주인공은 자신의 선택을 바꾸는 것이 항상 유리합니다.

토론

1. 문 1에 자동차가 있을 확률은 몬티홀씨의 행위에 상관없이 1/3 로 일정합니다. 그런데 나머지 문에 자동차가 있을 확률은 몬티홀씨의 행위에 영향을 받습니다. 왜냐하면 몬티홀씨는 나머지 두개의 문 중에서 어느 문에 자동차가 없는지를 보여줌으로써 확률을 한쪽으로 몰아주었기 때문입니다. 즉, 처음 확률 1/3 의 확률로 문2 와 문 3에 각각 자동차가 있을수 있었는데 몬티홀씨가 문 3을 열어서 자동차가 없음을 보여주게 되면 그 사건을 통해 문3에 자동차가 있을 확률은 0 이 되고 따라서 문2에 자동차가 있을 확률은 2/3가 되는 것입니다.

2. 만약 문이 10개 있다고 하고 주인공이 문 1을 선택했을때 몬티홀씨가 자동차가 없는 나머지 8개 문을 열어서 확인시켜 준다면 이 원리는 더 명확해 집니다.

3. 관찰을 통해 확률이 업데이트 되는 것은 학습(learning)의 주요 개념입니다. 위의 예제는 사전확률과 사후확률의 관계로도 표현됩니다. 특정 실험으로부터 어떤 결과가 얻어진다고 할때 그 결과를 얻기 이전에서의 확률은 사전확률이라 불리우고 그 결과를 바탕으로한 조건부 확률은 사후 확률이라고 부릅니다. 사후 확률을 구하는 방법을 베이즈 정리라고 부릅니다. 다음 시간에 좀더 자세히 다루도록 하겠습니다.

예제

어느 마을에 4개의 농장이 있고 그 농장에서 생산하는 곡물의 올해 생산량을 알고 싶다고 합시다. 편의상 다음과 같은 값을 가진다고 합시다.

ID	농지 면적	생산량 (톤)
1	4	1
2	6	3
3	6	5
4	20	15

위의 테이블에서 볼수 있듯이 농지 전체 면적은 36 이고 생산량은 24톤 될 것입니다. (농지 면적은 알고 있지만 생산량은 방문 조사하기 전까지는 알지 못한다고 합시다.) 그런데 예산상의 이유로 이 농장 4군데를 다 방문해서 생산량을 조사하지 못하고 2개의 농장 만을 조사해야 한다고 합시다. 이를 위해서 이 마을의 이장님은 농가 면적이 중간에 해당되는 농장 2와 농장 3을 방문해서 조사하려고 하였습니다. 그러나 그 마을에 사는 어느 통계학자는 그래서는 안되고 {1,2,3,4} 중에서 2개를 같은 확률로 램덤하게 뽑아서 그것에 해당되는 농장을 방문해야 한다고 주장했습니다. 왜 그래야 할까요?

풀이

위의 문제는 전형적인 샘플링 문제입니다. 통계학에서 검증된 샘플링 방법은 확률 표집법(probability sampling)입니다. 오늘은 이것을 공부해 보려고 합니다.

먼저 4개에서 2개를 뽑는 것은 다음과 같은 6가지 경우의 수가 있습니다.

Case	표본 ID	표본 평균	오차
1	1,2	2	-4
2	1,3	3	-3
3	1,4	8	2
4	2,3	4	-2
5	2,4	9	3
6	3,4	10	4

여기서 오차라는 것은 표본 평균이 모집단 평균(여기서는 24/4=6임)을 추정하고자 할때 발생하는 오차로써 표본 평균과 모평균의 차이를 나타냅니다. 위의 표에서 볼수 있듯이 어느 경우에도 오차는 발생하게 됩니다.

먼저 이장님이 생각한 표본 추출법은 위의 6가지 경우에서 한개를 주관적으로 결정하는 방법으로써 비확률 표집법(non-probability sampling)이라고 할수 있습니다.

반면 확률 추출법은 가능한 표본의 집합에 확률을 부여한후 그 확률에 따라 랜덤하게 뽑는 방법입니다. 위의 통계학자가 제안한 방법은 단순임의추출(simple random sampling)이라고 불리는 방법으로써 가장 단순한 형태의 확률추출법입니다. 이 방법은 모든 가능한 표본의 집합에 다음처럼 동일한 확률을 부여하는 것입니다.

Case	표본 ID	표본 평균	추출 확률
1	1,2	2	1/6
2	1,3	3	1/6
3	1,4	8	1/6
4	2,3	4	1/6
5	2,4	9	1/6
6	3,4	10	1/6

이렇게 추출확률을 부여하고 나서 그에 따라 표본을 뽑는 방법이 확률표본입니다. 예를 들어, 위의 단순임의추출 방법에서는 주사위를 굴려 1이 나오면 농가 1과 농가 2를 표본으로 뽑고 주사위 눈이 2가 나오면 농가 {1,3}이 표본 농가가 되는 것입니다. 이러한 표본의 확률 분포로부터 통계량의 확률분포가 유도될수 있는데 이를 표본 분포(sampling distribution)라고 합니다. 즉, 표본 평균이라는 통계량 $\hat{\theta}$ 의 확률 구조가 위의 표본 추출법으로부터 결정됩니다.

위의 테이블에서 볼수 있듯이 표본평균은 {2,3,4,8,9,10} 에서 하나를 각각 1/6의 확률로 그 값을 갖는 이산형 확률변수가 됩니다. 따라서 $E(\hat{\theta}) = \frac{1}{6} (2+3+4+8+9+10)=6$ 으로 되어서 모평균과 동일해 짐을 확인할수 있습니다. 여기서 기대값의 의미는 위의 sampling 을 여러번 반복했을때 얻어지게 되는 통계량 값들의 평균의 극한값으로 이해할수 있습니다.

또한 분산도 다음과 계산할수 있습니다. $V(\hat{\theta}) = \frac{1}{6}\{ (2-6)^2 + \cdots + (10-6)^2\} = 9.67 .$ 이렇게 확률 표본법은 가능한 표본에 추출확률을 부여하여 통계량의 표본 분포를 얻어냄으로써 기대값이나 분산과 같은 통계학적 계산을 가능하게 한다는 장점이 있습니다.

토론

1. 위의 사례는 확률 표본의 가장 간단한 예입니다. 확률 표본은 위의 예처럼 각 경우에 확률이 부여되어야 하고 각 원소가 뽑힐 확률이 0 보다 커야 합니다. 확률 표본은 비편향 추정 (unbiased estimation)을 구현할 수 있다는 대표적 장점이 있습니다. 확률 표본은 비편향 추정의 충분 조건입니다.

2. 모수 $\theta$의 추정량 $\hat{\theta}$ 의 오차(error)는 아래와 같이 나눌수 있습니다. \begin{eqnarray} \hat{\theta}-\theta &=& { \hat{\theta}-E( \hat{\theta} ) } + { E ( \hat{\theta} )- \theta } \ &=& \mbox{variation} + \mbox{bias } \end{eqnarray} 그러면 확률 표본의 경우 bias 가 0 이 되도록 추정량을 결정할수 있고 따라서 이 경우 variation 만 낮추어 주면 추정량의 정확성이 보장되는 것입니다. 한편 variation 을 낮추어 주는 가장 확실한 방법은 표본수를 증가시켜 주는 것입니다.

3. 게다가 확률 표본의 추가적 장점은 대수의 법칙이나 중심 극한 정리와 같은 중요한 통계학적 성질이 잘 규명되어 있다는 것입니다. 즉, 통계학 분야에서는 확률 표본에 대해서 많은 연구가 진행되었기 때문에 신뢰구간 건설이나 가설검정과 같은 통계적 추론들이 확률표본으로부터는 가능해 집니다.

4. 위의 단순 임의 추출은 확률 표본의 하나의 특수한 경우입니다. 아래와 같은 sampling design 도 또다른 형태의 확률 표본입니다.

표본 ID	표본값	추정량	추출 확률
1,4	1,15	4.5	1/3
2,4	3,15	6	1/3
3,4	5,15	7.5	1/3

이 경우에는 평균의 추정량이 {4,5, 6, 7.5}의 값을 각각 1/3 의 확률로 갖는 확률변수가 됩니다. 이러한 경우 그 기대값은 6 으로 모집단 평균과 동일해지며 (unbiased 추정) 분산은 $(1/3)*{ (1.5)^2 + 0^2 + (1.5)^2 } = 1.5$ 가 되어 simple random sampling 에서의 분산값인 9.67 보다 현저하게 작아집니다. 따라서 이 sampling design 이 더 효율적인 표본 설계가 되는 것입니다. (이 표본 설계는 층화 추출의 일종입니다. 이에 대해서는 다음 기회에 설명하도록 하겠습니다.)

.

예제

1975년 사이언스지에서는 ``Is there a sex bias in graduate admissions?’‘라는 제목의 흥미로운 논문이 발표되었습니다. 이 논문에서는 버클리 소재 캘리포니아 주립대학의 대학원 입시 자료를 바탕으로 남성과 여성 지원자간에 합격률의 차이가 있는지를 분석하였습니다. 먼저 이 자료에 의하면 특정기간 동안 현황을 살펴보면 8,422명의 남학생이 지원하여 44% 가 합격하였고 4,321명의 여학생이 지원하여 35% 가 합격했다고 합니다. 그러니 이 자료로부터 대학원 입시 과정에서 성차별이 존재하는 것이 아니냐고 합리적인 의심을 할수 밖에 없었을 것입니다.

그래서 전공 학과별로 자료를 분석하여 어떤 학과에서 문제가 있었던 것인지를 밝혀내기로 하고 다음과 같은 결과를 얻었다고 합니다.

성별	전공	지원자수	합격률 (%)
남성	A	825	62
남성	B	560	63
남성	C	325	37
남성	D	417	33
남성	E	191	28
남성	F	373	6

한편 여학생 자료는 다음과 같습니다.

성별	전공	지원자수	합격률 (%)
여성	A	108	82
여성	B	25	68
여성	C	593	34
여성	D	375	35
여성	E	393	24
여성	F	341	7

위의 6개 전공학과의 자료를 비교해 보면 남학생과 여학생간의 합격률 차이가 별로 없습니다. 오히려 여학생이 더 높은 것처럼 나타납니다. 뭐가 잘못된 것일까요?

풀이

위의 사례처럼 세부 집단으로 들어가서 보았을때 자료로 부터 얻어진 결론이 바뀌어서 나타나는 현상을 Simpson’s paradox라고 부릅니다. 이러한 파라독스가 발생하는 원인은 지원자들의 전공 선택과 성별이 교락(confounding)되었기 때문인데요 이를 좀더 자세히 설명하면 다음과 같습니다. 전공 A와 B는 합격률이 높은 전공이고 나머지 전공은 합격률이 낮은 전공입니다. 즉, 전공 A와 B는 들어가기 쉬운 학과인 것이지요. 위의 테이블에 의하면 남학생들은 이 전공에 많이 지원했지만 여학생들은 그렇지 않았습니다. 그래서 들어가기 쉬운 전공에 남학생들이 많이 지원해서 합격률이 높게 나타난것처럼 보이는 것입니다. 극단적으로 학과가 2개밖에 없다고 할때 A 학과에는 모두 남자가 지원하고 B 학과에는 모두 여자가 지원했다면 합격률의 차이가 남녀의 합격률의 차이에 의한 것인지 학과 자체의 합격률 차이에 의한 것인지 구분해 낼수 없습니다.

다시 말하자면 남학생들은 들어가기 쉬운 학과에 많이 지원했기 때문에 전체 남자 지원자 중에서 보면 합격률이 높이 나오는 것처럼 착시 현상이 나타난것 뿐입니다. 좀더 공평한 비교를 위해서는 위의 테이블처럼 전공별로 쪼개어 보던가 아니면 남녀별 합격률 계산을 제대로 다시해야 합니다.

위의 테이블에서 보다 공평한 남여별 합격률을 계산하기 위해서는 가중치를 남녀별 지원자수가 아닌 전체 지원자를 바탕으로 계산해야 합니다. 즉, 아래 테이블에서처럼 전체 지원자수를 계산하고 이를 바탕으로 합격률을 구하면 다음과 같습니다.

전공	전체 지원자수	남자 합격률 (%)	여자 합격률
A	933	62	82
B	585	63	68
C	918	37	34
D	792	33	35
E	584	28	24
F	714	6	7
Total	4,526	39	43

즉, 전체 지원자수를 바탕으로 합격률을 가중 평균하면 남자는 39%, 여자는 43% 으로 계산되어 오히려 여자의 합격률이 높은 것으로 나타납니다. 성별 지원자수가 아닌 전체 지원자수를 바탕으로 가중평균을 구하는 방법은 교락요인을 통제해 주는 효과를 가져옵니다.

토론

1. 위의 데이터는 관측 연구(observational study) 에서 얻어진 것입니다. 실험(experiment)가 아닌 관측 연구 자료의 분석은 교락요인이 통제가 되지 않아서 잘못된 결론을 얻기가 매우 싶습니다. 일반적으로 말해서 인과관계는 실험을 통해서 입증을 할수 있고 관측연구에서는 인과(causality) 관계가 아닌 상관(association) 관계 만을 밝혀낼 뿐입니다.

2. 관측 연구에서 가장 어려운 점은 모든 교락요인을 다 통제하지 못할수 있다는 것입니다. 실험이 아니므로 randomization 을 이용한 통제가 불가능하고 사후적으로 분석 단계에서 통제를 해야 하는데 문제는 관측되지 않는 교락 요인이 있을수 있다는 것입니다.

3. 위의 예에서는 전공 선택이 명백한 교략요인이었기에 이를 통제하기가 쉬웠는데 다른 경우에는 그게 명백하지 않을수 있습니다. 회귀분석(regression analysis)을 통해서 관측연구 자료를 분석할때 관심 요인 외에도 다른 요인을 설명변수에 넣어주는 것도 그러한 요인을 통제하고자 하는 것입니다. 즉, 다른 (관측된) 요인이 동일하다고 했을때 관심 요인이 반응 변수에 얼마나 영향을 미치는지 보고자 하는 것입니다. 그런데 만약 중요한 교락 요인이 이 회귀 분석에서 누락되었다면 인과관계 입증은 실패하는 것입니다. 유명한 통계학자 R.A. Fisher 경도 그러한 이유로 담배가 폐암의 직접적인 원인이 된다는 것을 받아들이지 않았습니다.

예제

소아마비 백신을 개발한 조너스 소크 박사 이야기를 들어보셨나요? (아래 링크 참조 https://ko.wikipedia.org/wiki/조너스_소크)

조너스 소크 박사와 그의 연구팀이 개발한 소아마비 백신은 공식적으로 사용되기 이전에 그 백신이 안전하고 효과적임을 증명하는 임상실험을 통과해야 했습니다.

소크가 백신을 개발한 날이 1952년 3월 26일이었으며, 소크는 자신에게 이 백신을 처음 사용해 성공을 거둔 후, 공개적으로 지원자를 모집해 효과를 검증하기 시작했습니다. 1955년 4월 12일, 미국 44개주에서 소크 백신은 부작용이 없는 아주 효과가 좋은 소아마비 예방을 위한 백신이라는 연구결과가 최종 발표되기까지 보건당국은 많은 노력을 기울여 실험을 해야 했습니다.

처음의 임상 실험 계획은 자원자 모두에게 백신을 투여하여 백신을 받은 어린이들의 소아마비 발병률이 예전보다 더 낮아지는 지를 비교하는 것이었다고 합니다. 하지만 보건통계학자들은 그렇게 해서는 안되고 지원자 집단을 랜덤하게 두 그룹으로 나누어서 하나에는 백신을 주고 다른 하나에는 가짜약(placebo)을 주어서 그 두 집단간의 소아마비 발병률의 차이를 비교해야 한다고 주장했습니다. 왜 그래야 하는 것일까요?

풀이

통계학에서 가장 중요한 개념 중의 하나는 통제(control)입니다. 비교를 제대로 하려면 동일한 조건에서 비교를 해야 한다는 것이지요. 즉, 실험군(treatment group)과 대조군(control group)과의 차이는 treatment 여부 외에는 없어야 한다는 것입니다. 그래야 인과관계에 대한 과학적인 결론을 얻어낼수 있다는 것입니다.

만약 처음 계획처럼 자원자 모두에게 백신을 투여한다면 자원자들의 속성이 비자원자들의 속성과 체계적으로 다를수 있습니다. (실제로 자원자 집단의 가계 소득이 평균보다 높았다고 합니다.) 그러면 실험군의 소아마비 발병률이 대조군보다 낮게 나온다 하더라도 이게 백신 때문인지 아니면 소득에 따른 영양상태의 차이와 같은 다른 속성 때문인지 구분해 낼수가 없습니다. 이렇게 두 요인이 엉켜 있어서 자료를 통해서 관심 요인에 의한 효과를 구분해 낼수 없는 경우를 교락(confounding)되었다고 표현합니다. 이런 confounding effect 를 통제하기 위해서 동일한 조건의 대조군을 설정해야 하는데 randomization 이라는 마술이 이를 가능하게 한다는 것입니다. 두 집단을 나눌때 동전을 던저서 앞면이 나오면 실험군에 속하게 하고 뒤면이 나오면 대조군에 속하게 하면 우연히 두 집단의 속성이 다를 가능성이 표본수가 커지면서 급격하게 낮아지게 되는 것입니다.

토론

1. 통계학은 과학적인 비교를 하고자 하는 학문입니다. 비교의 핵심은 비교 대상을 잘 찾는 것입니다. 만약 비교 대상이 모든 조건에서 잘 통제되었다면 실험후 두 집단의 차이는 treatment 에 의한 효과일 확률이 큽니다. 통제된 실험을 하기 위해서는 random 하게 두 그룹을 나누어 주는 방법이 사용될수 있습니다. 이런 실험을 randomized controlled experiment 라고 합니다.

2. Randomized controlled experiment 는 실험 이전의 두 집단간의 차이를 줄여줍니다. 그리고는 실험 과정의 차이를 없애기 위해서 treatment 집단과 control 집단에게 동일한 방식으로 주사를 하는 것이 좋습니다. 또한 의사 역시 이 환자가 진짜 백신을 주사받는 것인지 가짜 백신을 주사받는 것인지 모르게 하는 것이 좋습니다. 즉, 환자 뿐만 아니라 의사 역시 이 환자가 treatment 집단인지 control 집단인지 모르도록 하는 것이 실험 과정에서 발생할수 있는 다른 요인을 통제하는데에 좋습니다. 이런 방법을 double-blinded experiment 라고 합니다.

3. 이렇게 해서 얻어진 자료를 분석하기 위해서 T-test 라는 것을 씁니다. 두 집단의 표본 평균의 차이가 얼마나 통계적으로 유의한 차이가 있는지를 검정하는 방법입니다. 이런 T-test 가 가장 큰 검정효과를 가져오기 위해서는 두 집단의 표본수를 같게 해 주는게 좋습니다.

4. 현실 세계에서는 위의 randomized control experiment 가 항상 가능한 것은 아닙니다. 예를 들어, 담배의 유해성을 알아보기 위해 사람들을 두 집단으로 나누어서 한 집단에는 담배를 피게 하고 다른 집단에는 담배를 못피우게 할수 없는 것입니다. 실험이 아닌 자료에서는 다른 요인을 통제하기가 어렵습니다. 이에 대해서는 다음 시간에 이야기 하도록 하겠습니다.

---

layout: post title: “기대값의 의미” author: “JKKim” date: “August 26, 2016”

comments: true share: true

---

예제

회사원 홍길동 씨는 두달후에 있을 여행을 위해 비행기표를 구입했습니다. 예를들어 백만원짜리 비행기표라고 합시다. 그런데 이 비행기 표가 non-refundable 이라서 갑자기 부득이한 사정이 생겨서 비행기를 타지 못하는 상황이 되어도 환불을 받지 못하고 그냥 돈 백만원을 날리게 된다는 것입니다. 이러한 상황을 대비해서 보험을 들고자 해서 알아봤더니 5만원을 내야 한다고 합니다. 이때 홍길동씨는 5만원을 내고 보험을 드는 것이 더 나을까요? 여기서 보험금 5만원은 어떤 의미일까요?

풀이

보험금 5만원은 기대 손실(expected loss)에 대한 가격에 보험회사의 수수료 및 이윤이 추가된 금액입니다. 홍길동씨가 두달후에 있을 비행기 출발 시점에 비행기를 타지 못할 확률을 $p$ 라고 할때 다음과 같은 손실 함수를 가지게 됩니다.

$$L= \left\{ \begin{array}{ll} 1,000,000 & \mbox{ with probability } p \\ 0 & \mbox{ with probability } 1-p \end{array} \right.$$

즉, 홍길동씨의 손실 함수 $L$ 은 두 개의 값을 가지는 확률 변수 (random variable) 이 되는 것입니다. 이러한 손실 함수에 대한 기대값은 다음과 같이 계산될 것입니다.

$$E(L) = 1,000,000 \times p + 0 \times (1-p) = 1,000,000 \cdot p .$$

따라서 홍길동씨 입장에서는 자신이 비행기를 못타게 될 확률이 5\% 보다 더 크다고 생각이 든다면 이 보험을 드는 것이 더 합리적인 행동인 것입니다.

토론

1. 위의 예제에서 $L$ 은 두개의 값을 가지는 확률변수입니다. 이렇게 셀수 있는 갯수의 (countable한) 값을 가지는 확률 변수를 범주형 (discrete) 확률변수라고 하고 몸무게나 소득처럼 연속적인 값을 가지는 확률 변수를 연속형 (continuous) 확률 변수라고 합니다. 범주형 확률 변수의 기대값은 위의 예제에서처럼 $E(X) = \sum x_i P(X=x_i)$ 의 형태로 나타나고 연속형 확률 변수의 기대값은 적분을 사용해서 계산합니다.

2. 보험회사 입장에서는 확률 $p$ 를 좀더 정교하게 계산하는 것이 중요할 것입니다. 이 확률을 모든 사람에게 동일하다고 가정하는 것보다는 연령대와 성별 등에 따라 다르게 계산하는 경우를 생각해 볼수 있는데 이때의 확률은 조건부 확률이라고 부르고 $p(x)$ 로 표시할수 있을 것입니다. 여기서 $x$ 는 성별이나 연령대 같은 조건값이 되겠습니다. 조건부 확률을 사용하여 계산된 기대값을 조건부 기대값(conditional expectation)이라고 합니다.

3. 보험회사 입장에서는 이러한 조건부 확률을 정확하게 추정하는 것이 매우 중요합니다. 조건부 확률을 지나치게 크게 추정하면 보험 상품이 비싸져서 보험이 안팔리게 되니 매출이 낮아지고 조건부 확률을 지나치게 낮게 추정하면 보험 상품은 싸지겠지만 보험금 지불액이 예상보다 높아지므로 손실이 커지게 됩니다. 따라서 보험회사에서는 이러한 조건부 확률을 좀더 정확하게 추정하기 위해서 양질의 데이터를 구입하여 분석하고자 하는 욕구가 생기는 것입니다.

4. 조건부 확률을 자료로부터 얻어내는 작업을 통계학에서는 추정 (estimation) 이라고 하고 전산과학 쪽에서는 학습 (learning)이라고도 부릅니다. 결국 특정 사건이 일어날 확률에 대한 지식을 데이터를 통해서 알아내는 것이 통계학의 기본 작업입니다. 이런 사례는 무궁무진 합니다. 은행에서 대출신청을 받아서 대출 여부 및 이자율을 결정할때에도 조건부 파산 확률을 계산해서 결정해야 하고 정부에서 곡물 생산량을 미리 추정해서 곡물값이 폭등하거나 폭등하지 않도록 수출/수입을 통제하는 것도 결국은 이러한 추정 작업을 통해서 결정하게 됩니다. 즉, 통계학은 정확한 조건부 확률을 계산하여 최선의 결정을 내리고자 하는 학문입니다.

5. 만약 위의 손실함수에서 특정 손실값이 무한대에 가까우면 그 확률이 아무리 낮더라도 기대값 역시 무한대 입니다. 무한대 곱하기 $p(>0)$ 는 무한대이기 때문입니다. 원자력 발전소 건설 같은 것이 그러한 예가 될 것입니다. 원자력 발전소가 고장나거나 폭파 되는 경우가 그 확률이 매우 낮지만 그 확률이 0 이 아닌 이상 기대 손실에 반영됩니다. 원자력 발전소가 폭파되어서 수백만의 생명이 죽는다는 것은 무한대의 값을 가지는 손실이라고 볼수 있습니다. 그래서 이러한 원자력 발전소 건설에 대한 기대손실은 무한대 입니다. 하지만 기대 이익은 유한합니다. 따라서 통계학적 상식으로는 원자력 발전을 하지 않는것이 합리적인 결정입니다.