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

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-\bar{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 \)에 대하여 아래의 식이 성립한다.

$$ E(X^m)=\int_{0}^{\infty}mx^{m-1}(1-F(x))dx+\int_ {-\infty}^{0}mx^{m-1}F(x)dx \label{basic01}\tag{1} $$

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

\begin{eqnarray*} X^{m}&=&(Y+Z)^{m} \\&=& Y^{m}+Z^{m} \because YZ=0 \\&=& \int_{0}^{Y}mx^{m-1}dx+\int_{Z}^{0}mx^{m-1}dx \\&=& \int_{0}^{\infty}mx^{m-1}I(x\leq Y)dx+\int_{-\infty}^{0}mx^{m-1}I(x\geq Z)dx \end{eqnarray*}

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

\begin{eqnarray*} E(X^{m})&=&E(Y^{m}) + E(Z^{m}) \\ &=& \int_{0}^{\infty}\int_{0}^{\infty}mx^{m-1}I(x\leq y)dxf_X(y)dy+\int_{-\infty}^{0}\int_{-\infty}^{0}mx^{m-1}I(x\geq z)dxf_X(z)dz \\ &=& \int_{0}^{\infty}\int_{0}^{\infty}mx^{m-1}I(x\leq y)f_X(y)dydx+\int_{-\infty}^{0}\int_{-\infty}^{0}mx^{m-1}I(x\geq z)f_X(z)dzdx \\ &=& \int_{0}^{\infty}mx^{m-1}\int_{0}^{\infty}I(x\leq y)f_X(y)dydx+\int_{-\infty}^{0}mx^{m-1}\int_{-\infty}^{0}I(x\geq z)f_X(z)dzdx \\ &=& \int_{0}^{\infty}mx^{m-1}\int_{x}^{\infty}f_X(y)dydx+\int_{-\infty}^{0}mx^{m-1}\int_{-\infty}^{x}f_X(z)dzdx \\ &=& \int_{0}^{\infty}mx^{m-1}(1-F(x))dx+\int_ {-\infty}^{0}mx^{m-1}F(x)dx \quad\blacksquare \end{eqnarray*}

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

$$ F_R(r)=n\int_{-\infty}^{\infty}(F(v+r)-F(v))^{n-1}f_X(v)dv \label{basic02}\tag{2} $$

먼저 \(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_{-\infty}^{\infty}n(n-1)(F(v+r)-F(v))^{n-2}f_X(v)f_X(v+r)dv\)

이제 \(r\)에 관하여 적분하자.

\begin{eqnarray*} F_R(r) &=& \int_{0}^{r} \int_{-\infty}^{\infty}n(n-1)(F(v+w)-F(v))^{n-2}f_X(v)f_X(v+w)dvdw \\ &=& \int_{-\infty}^{\infty}n(n-1)\int_{0}^{r} (F(v+w)-F(v))^{n-2}f_X(v+w)dwf_X(v)dv \\ &=& \int_{-\infty}^{\infty}n\{(F(v+w)-F(v))^{n-1}|_{w=0}^{w=r}\}f_X(v)dv \\ &=& n\int_{-\infty}^{\infty}(F(v+r)-F(v))^{n-1}f_X(v)dv \quad\blacksquare \end{eqnarray*}

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

\begin{eqnarray*} E(r)&=&\sigma\int_{0}^{\infty}1-  n\int_{-\infty}^{\infty}(\Phi(v+r)-\Phi(v))^{n-1}\Phi'(v)dvdr \\&\equiv& d_2(n)\sigma \end{eqnarray*}

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

\begin{eqnarray*} E(s) &=& \sqrt{\frac{\sigma^2}{n-1}}E\left( \sqrt{\frac{(n-1)s^2}{\sigma^2}}\right) \\[3pt] &=& \frac{\sigma}{\sqrt{n-1}} \int_{0}^{\infty}\sqrt{x}\frac{1}{2^{(n-1)/2}\Gamma\left((n-1)/2\right)}x^{(n-1)/2-1}\exp(-x/2)dx \\[3pt] &=& \frac{\sigma}{\sqrt{n-1}} \frac{2^{n/2}\Gamma(n/2)}{2^{(n-1)/2}\Gamma((n-1)/2)}\int_{0}^{\infty}\sqrt{x}\frac{1}{2^{n/2}\Gamma(n/2)}x^{n/2-1}\exp(-x/2)dx \\[3pt] &=& \sigma\sqrt{\frac{2}{n-1}}\frac{\Gamma(n/2)}{\Gamma\left((n-1)/2\right)} \\[3pt] &\equiv& c_4(n)\sigma \end{eqnarray*}

\(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{eqnarray*} Var(\hat\sigma_r)&=& Var\left(\ \frac{r}{d_2(n)} \right) \\[3pt] &=& \frac{1}{d_2^2(n)} \left\{ E(r^2)-E^2(r) \right\} \\[3pt] &=& \frac{1}{d_2^2(n)} \left\{ \int_{0}^{\infty}2r \left(1- n\int_{-\infty}^{\infty}(\Phi(v+r)-\Phi(v))^{n-1}\Phi'(v)dv \right) dr -d_2^2(n)\right\}\sigma^2 \\[3pt] &\equiv& d_3(n)\sigma^2 \end{eqnarray*}

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

\begin{eqnarray*} Var(\hat\sigma_r)&=&Var\left( \frac{s}{c_4(n)} \right) \\[3pt] &=&\frac{1}{c_4^2(n)}\left\{ E(s^2)-E^2(s)\right\} \\[3pt] &=&\frac{1}{c_4^2(n)}(1-c_4^2(n))\sigma^2 \\[3pt] &\equiv&c_5(n)\sigma^2 \end{eqnarray*}

\(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 [본문으로]