표본 범위와 표본 표준편차 비교

1. 개요

이것저것을 하다보니 통계적 공정 관리(Statistical Process Control; SPC)를 접할 일이 있었다. SPC란 통계적인 방법론을 사용하여 품질을 관리하는 시스템이다. 1920년대 초반에 아이디어가 나왔고 지금까지 많은 공장에서 사용되고 있다. SPC의 구체적인 내용은 이번 포스팅의 주제가 아니므로 다루지 않고 링크로 대체한다. SPC에서 나의 관심을 끈 것은 아래와 같은 그림이었다.^[각주:1] 

 

Sample size($n$)가 2 이상, 10 미만인 경우 Xbar-R chart를, 10 이상일 때에는 Xbar-S chart를 사용하라고 한다.^[각주:2] 즉, $n$이 작을 때에는 표본의 범위($r$)로 모표준편차($\sigma$)를 추정하고 $n$이 클 때에는 표본의 표준편차($s$)로 $\sigma$를 추정하자는 것이다.

$$r=max(x_i)-min(x_i)$$ $$s=\sqrt{s^2}=\sqrt{\frac{1}{n-1}\sum _{i=1}^n(x_i-\overline{x}{)}^2}$$

$r$을 이용한 $\sigma$의 추정량(estimator)을 ${\hat{\sigma }}_r$, $s$를 이용한 $\sigma$의 추정량을 ${\hat{\sigma }}_s$라고 하자. 그렇다면 n의 크기에 따른 각 추정량의 성질이 어떻게 변하길래 그러한 기준이 생긴 것일까? $X\sim N(0,{\sigma }^2)$를 가정하고 진행해보자.

2. 어떤 추정량이 더 좋은 추정량인가?

2.1. 불편(unbiased)성 비교

가장 먼저 떠오르는 것은 $r$과 $s$의 기대값이 무엇이냐는 것이다. 만약 Bias가 있다면, correction이 가능한지도 궁금했다. $E(r)$를 구하기 전에 필요한 두가지 공식을 정리하자. 

임의의 Random Variable $X$의 누적분포함수(Cumulative Distribution Function; cdf)를 $F(x)$, 확률밀도함수(Probability Density function; pdf)를 $f_X(x)$라고 하면 0이 아닌 실수 $m$에 대하여 아래의 식이 성립한다.

$$\begin{array}{}\text{(1)}& E(X^m)={\int }_0^{\mathrm{\infty }}m{x}^{m-1}(1-F(x))dx+{\int }_{-\mathrm{\infty }}^{0}m{x}^{m-1}F(x)dx\end{array}$$

$Y$를 $X$의 positive part 즉, $max(X,0)$라고 하고 $Z$를 $X$의 negative part, $min(X,0)$이라고 하자. 단, $I$는 indicator function이다.

$$\begin{array}{rcl}{X}^m& =& (Y+Z{)}^m\\ & =& Y^m+Z^m\because YZ=0\\ & =& {\int }_0^{Y}m{x}^{m-1}dx+{\int }_Z^{0}m{x}^{m-1}dx\\ & =& {\int }_0^{\mathrm{\infty }}m{x}^{m-1}I(x\le Y)dx+{\int }_{-\mathrm{\infty }}^{0}m{x}^{m-1}I(x\ge Z)dx\end{array}$$

따라서, $E(X^m)$을 다음과 같이 쓸 수 있다.

$$\begin{array}{rcl}E(X^m)& =& E(Y^m)+E(Z^m)\\ & =& {\int }_0^{\mathrm{\infty }}{\int }_0^{\mathrm{\infty }}m{x}^{m-1}I(x\le y)dx{f}_X(y)dy+{\int }_{-\mathrm{\infty }}^0{\int }_{-\mathrm{\infty }}^{0}m{x}^{m-1}I(x\ge z)dx{f}_X(z)dz\\ & =& {\int }_0^{\mathrm{\infty }}{\int }_0^{\mathrm{\infty }}m{x}^{m-1}I(x\le y)f_X(y)dydx+{\int }_{-\mathrm{\infty }}^0{\int }_{-\mathrm{\infty }}^{0}m{x}^{m-1}I(x\ge z)f_X(z)dzdx\\ & =& {\int }_0^{\mathrm{\infty }}m{x}^{m-1}{\int }_0^{\mathrm{\infty }}I(x\le y)f_X(y)dydx+{\int }_{-\mathrm{\infty }}^{0}m{x}^{m-1}{\int }_{-\mathrm{\infty }}^{0}I(x\ge z)f_X(z)dzdx\\ & =& {\int }_0^{\mathrm{\infty }}m{x}^{m-1}{\int }_x^{\mathrm{\infty }}{f}_X(y)dydx+{\int }_{-\mathrm{\infty }}^{0}m{x}^{m-1}{\int }_{-\mathrm{\infty }}^x{f}_X(z)dzdx\\ & =& {\int }_0^{\mathrm{\infty }}m{x}^{m-1}(1-F(x))dx+{\int }_{-\mathrm{\infty }}^{0}m{x}^{m-1}F(x)dx◼\end{array}$$

$r$의 cdf $F_R(r)$는 아래와 같다.

$$\begin{array}{}\text{(2)}& F_R(r)=n{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}(F(v+r)-F(v){)}^{n-1}{f}_X(v)dv\end{array}$$

먼저 $r$의 pdf를 먼저 구하자. order statistics를 이용하면 최소값과 최대값의 joint pdf를 아래와 같이 구할 수 있다.

$f_{x_{(1)},x_{(n)}}(x,y)=n(n-1)(F(y)-F(x){)}^{n-2}{f}_X(x)f_X(y)$

그리고 variable change($r=y-x,v=x$) 및 $v$에 관한 적분을 통하여 $r$의 pdf를 다음과 같이 구할 수 있다.

$f_R(r)={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}n(n-1)(F(v+r)-F(v){)}^{n-2}{f}_X(v)f_X(v+r)dv$

이제 $r$에 관하여 적분하자.

$$\begin{array}{rcl}{F}_R(r)& =& {\int }_0^r{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}n(n-1)(F(v+w)-F(v){)}^{n-2}{f}_X(v)f_X(v+w)dvdw\\ & =& {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}n(n-1){\int }_0^r(F(v+w)-F(v){)}^{n-2}{f}_X(v+w)dw{f}_X(v)dv\\ & =& {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}n\{(F(v+w)-F(v){)}^{n-1}{|}_{w=0}^{w=r}\}{f}_X(v)dv\\ & =& n{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}(F(v+r)-F(v){)}^{n-1}{f}_X(v)dv◼\end{array}$$

(1)과 (2)를 이용하면 $E(r)$을 구할 수 있다. 단, $\mathrm{\Phi }$는 표준정규분포의 cdf이다.

$$\begin{array}{rcl}E(r)& =& \sigma {\int }_0^{\mathrm{\infty }}1-n{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}(\mathrm{\Phi }(v+r)-\mathrm{\Phi }(v){)}^{n-1}{\mathrm{\Phi }}^{\prime }(v)dvdr\\ & \equiv & d_2(n)\sigma \end{array}$$

$E(s)$는 정규분포와 ${\chi }^2-$ 분포의 관계를 이용하여 쉽게 구할 수 있다. 

$$\begin{array}{rcl}E(s)& =& \sqrt{\frac{{\sigma }^2}{n-1}}E\left(\sqrt{\frac{(n-1)s^2}{{\sigma }^2}}\right)\\ & =& \frac{\sigma }{\sqrt{n-1}}{\int }_0^{\mathrm{\infty }}\sqrt{x}\frac{1}{2^{(n-1)/2}\mathrm{\Gamma }((n-1)/2)}{x}^{(n-1)/2-1}\exp (-x/2)dx\\ & =& \frac{\sigma }{\sqrt{n-1}}\frac{2^{n/2}\mathrm{\Gamma }(n/2)}{2^{(n-1)/2}\mathrm{\Gamma }((n-1)/2)}{\int }_0^{\mathrm{\infty }}\sqrt{x}\frac{1}{2^{n/2}\mathrm{\Gamma }(n/2)}{x}^{n/2-1}\exp (-x/2)dx\\ & =& \sigma \sqrt{\frac{2}{n-1}}\frac{\mathrm{\Gamma }(n/2)}{\mathrm{\Gamma }((n-1)/2)}\\ & \equiv & c_4(n)\sigma \end{array}$$

$c_4(n)$와 $d_2(n)$ 모두 $n$에 의해서만 결정되는 상수이므로 ${\hat{\sigma }}_r=r/d_2(n)$, ${\hat{\sigma }}_s=s/c_4(n)$라고 하면 ${\hat{\sigma }}_r$, ${\hat{\sigma }}_s$ 모두 불편추정량이며, 서로 우열이 없음을 알 수 있다. 

2.2. 분산 비교

그 다음으로 생각난 것은 분산이다. 특정 n을 기준으로 ${\hat{\sigma }}_r$과 ${\hat{\sigma }}_s$의 분산이 어떻게 달라지는가? 먼저 식 (1)과 식 (2)를 이용하여 $Var({\hat{\sigma }}_r)$를 구한다.

$$\begin{array}{rcl}Var({\hat{\sigma }}_r)& =& Var\left(\frac{r}{d_2(n)}\right)\\ & =& \frac{1}{d_2^2(n)}\{E(r^2)-E^2(r)\}\\ & =& \frac{1}{d_2^2(n)}\{{\int }_0^{\mathrm{\infty }}2r(1-n{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}(\mathrm{\Phi }(v+r)-\mathrm{\Phi }(v){)}^{n-1}{\mathrm{\Phi }}^{\prime }(v)dv)dr-d_2^2(n)\}{\sigma }^2\\ & \equiv & d_3(n){\sigma }^2\end{array}$$

$Var({\hat{\sigma }}_r)$은 정규분포와 ${\chi }^2-$분포의 관계를 이용하여 구한다.

$$\begin{array}{rcl}Var({\hat{\sigma }}_r)& =& Var\left(\frac{s}{c_4(n)}\right)\\ & =& \frac{1}{c_4^2(n)}\{E(s^2)-E^2(s)\}\\ & =& \frac{1}{c_4^2(n)}(1-c_4^2(n)){\sigma }^2\\ & \equiv & c_5(n){\sigma }^2\end{array}$$

$d_3(n)$과 $c_5(n)$ 모두 대소를 구분할 정도의 근사값은 구할 수 있다. R로 정수 $n\in [2,50]$에 대하여 $d_3(n)/c_5(n)$의 그래프를 그려보자. 결과는 아래와 같다. 단, $\sigma =1$로 하며 그래프를 그리기 위한 코드를 첨부한다. 

range_issue.r

$n=2$일때는 분산이 서로 같고 3 이상의 $n$에 대하여 ${\hat{\sigma }}_s$가 더 우월하다는 것을 알 수 있다. 작은 $n$에 대하여 $r$이 어떠한 이점을 가지고 있는지 알기 어려웠다. 

3. $r$이 $s$보다 좋은 이유

계속해서 이유를 찾고 찾다가 한 책^[각주:3]에서 그 이유를 발견했다.

Traditionally, quality engineers have preferred the R chart to the s chart because of the simplicity of calculating R from each sample. The availability of hand-held calculators with automatic calculation of s and computers at workstations to implement control charts on site have eliminated any computational difficulty.

사실 $r$이 $s$보다 계산하기 쉽다는 것은 이미 알고 있었고 위키피디아의 Xbar-R chart 등에서 꾸준이 본 내용이기도 했다. 하지만 나는 SPC를 시스템으로 접근하고 있었기 때문에 계산이 쉽다는 말 자체를 이해하지 못했다. 그런데 Traditionally라는 단어를 보자마자 SPC가 1920년대에 아이디어가 나왔다는 것을 떠올렸다. 그 당시에는 $r$이 가지고 있는 계산상의 이점이 통계적인 약점보다 두드러졌을 것이다. 

결국 $r$을 사용하는 이유는 관습 때문이었다. 이젠 더이상 쓸 필요가 없는 것이다. 누군가가 말한대로 컨텍스트는 가고 텍스트만 남은 꼴이다. 

4. 마치며

$E(r)$ 및 $Var(r)$은 위의 식보다 더 간소하게 표현할 수 있다. Barbosa 등(2013)^[각주:4] 및 Tippet (1925)^[각주:5]에 그 내용이 있다. 안타깝게도 나는 En taro Adun하여 이해하지 못했으나 누군가는 할 수 있을 것이다.

허무한 결론이지만 어찌됬든 이유를 알았다는 점은 만족스러웠다. 또한 이유를 알아내는 과정에서 많은 공식을 증명하면서 수식 전개를 연습했다는 것도 좋았다. 그냥 버리기는 아까우니 증명한 공식은 차후에 따로 포스팅을 할 예정이다. 

1. 출처는 http://www.sixsigmatrainingfree.com/statistical-process-control-charts.html이며 가독성을 위하여 붉은 상자를 추가하였다. [본문으로]
2. 각 용어에 대한 설명도 이번 포스팅에 대한 주제가 아니므로 링크로 대체한다. [본문으로]
3. Introduction to Statistical Quality Control, Wiley; 6 edition, Douglas C. Montgomery, 2008, 254p [본문으로]
4. Range Control Charts Revisited: Simpler Tippett-like Formulae, Its Practical Implementation, and the Study of False Alarm, 2013, Emanuel Pimentel Barbosa 등, Communication in Statistics- Simulation and Computation 42(2) [본문으로]
5. On the Extreme Individuals and the Range of Samples from a Normal Population. 1925, Tippett, Biometrika, vol 17, 364-387 [본문으로]