데이타 처리를 하면서 평균과 표준편차를 구해야 하는 경우가 많이 있다. 이 때, 평균과 표준편차의 계산을 점화식을 이용해서 하면 반복문의 수를 반으로 줄여서 빠르게 계산을 할 수 있다.

• 증명따윈 필요없는 귀차니스트를 위한 식

평균: \(\mu_{n+1} = \displaystyle\frac{n}{n+1}\mu_n + \displaystyle\frac{a_{n+1}}{n+1},\quad(n\geq0). \)

표준편차: \(\sigma^2_{n+1}=\displaystyle\frac{n}{n+1}\sigma^2_{n}+\displaystyle\frac{(\mu_{n+1}-a_{n+1})^2}{n},\quad(n\geq1).\)

공변: \(\sigma(a,b)_{n+1}=\displaystyle\frac{n}{n+1}\sigma(a,b)_{n}+\displaystyle\frac{(\mu_{n+1}-a_{n+1})(\nu_{n+1}-b_{n+1})}{n},\quad(n\geq1).\)

• 평균 증명
  • 점화식을 이용해서 \(n+1\)번째 항까지의 평균을 구할 때 우리가 가지고 있는 정보는 다음과 같다.
    1. \(n\)번째 항까지의 평균 \(\mu_n\)
    2. \(n+1\)번째 항의 값 \(a_{n+1}\)
  • 이 정보들만 가지고 \(n+1\)번째 항까지의 평균을 구해내야 한다.

\[\begin{align}\mu_{n+1}&=\displaystyle\sum^{n+1}_{i=1}\displaystyle\frac{a_i}{n+1},\\&=\displaystyle\sum^{n}_{i=1}\displaystyle\frac{a_i}{n+1}+\displaystyle\frac{a_{n+1}}{n+1},\\&=\displaystyle\frac{a_{n+1}}{n+1}+\displaystyle\frac{n}{n+1}\displaystyle\sum^{n}_{i=1}\displaystyle\frac{a_i}{n}.\end{align}\]

\[\therefore \mu_{n+1}=\dfrac{n}{n+1}\mu_n+\dfrac{a_{n+1}}{n+1},\quad(n\geq0).\]

 

• 표준편차 증명
  • \(n+1\)번째 항까지의 표준편차를 구할 때에 우리가 가지고 있는 정보는 다음과 같다.
    1. \(n\)번째 항까지의 표준편차 \(\sigma^2_{n}\)
    2. \(n+1\)번째 항까지의 평균 \(\mu_{n+1}\)
    3. \(n+1\)번째 항의 값 \(a_{n+1}\)
  • 평균과 마찬가지로, 표준편차를 구할 때도 위 세 가지 값만 가지고 \(n+1\)번째 항까지의 표준편차를 구해야 한다.

\[\begin{align}\sigma^2_{n+1}&=\sum^{n+1}_{i=1}\dfrac{(a_i-\mu_{n+1})^2}{n+1},\\&=\sum^{n}_{i=1}\dfrac{(a_i-\mu_{n+1})^2}{n+1}+\frac{(a_{n+1}-\mu_{n+1})^2}{n+1}.\qquad (1)\end{align}\]

두번째 항은 알고있는 값들로만 돼 있으니 놔두고, 첫번째 항을 알고있는 값들로만 바꿔주면 된다.

\[\begin{align}\sum^{n}_{i=1}(a_i-\mu_{n+1})^2-\sum^{n}_{i=1}(a_i-\mu_n)^2&=\sum^{n}_{i=1}(\mu_n-\mu_{n+1})(2a_i-\mu_{n+1}-\mu_n),\\&=(\mu_n-\mu_{n+1})(2n\mu_n-n\mu_{n+1}-n\mu_n),\\&=n(\mu_{n+1}-\mu_n)^2.\qquad(2)\end{align}\]

이제 식 (2)를 식 (1)에 넣어 정리하면,

\[\begin{align}\sigma^2_{n+1}&=\sum^{n}_{i=1}\dfrac{(a_i-\mu_n)^2}{n+1}+\dfrac{n}{n+1}(\mu_{n+1}-\mu_n)^2+\dfrac{(a_{n+1}-\mu_{n+1})^2}{n+1},\\&=\dfrac{n}{n+1}\sigma^2_{n}+\dfrac{n}{n+1}(\mu_{n+1}-\mu_n)^2+\dfrac{(a_{n+1}-\mu_{n+1})^2}{n+1}.\qquad(3)\end{align}\]

\(n\)번째 항까지의 평균값은 \(n+1\)번째 항까지의 평균값을 구하면서 없어졌으므로, 평균의 점화식 공식을 이용해서 \(n\)번째 항까지의 평균값을 없애버리자. 점화식을 다시 정리하면

\[\mu_n=\dfrac{n+1}{n}\mu_{n+1}-\dfrac{1}{n}a_{n+1}\]

가 되므로, 이 식을 식 (3)의 두번째 항에 대입하면, 두번째 항은 다음과 같이 된다.

\[\begin{align}\dfrac{n}{n+1}(\mu_{n+1}-\mu_n)^2&=\dfrac{n}{n+1}\left(\mu_{n+1}-\dfrac{n+1}{n}\mu_{n+1}+\dfrac{1}{n}a_{n+1}\right)^2,\\&=\dfrac{1}{n(n+1)}(a_{n+1}-\mu_{n+1})^2\end{align}\]

\[\therefore \sigma^2_{n+1}=\dfrac{n}{n+1}\sigma^2_{n}+\dfrac{(\mu_{n+1}-a_{n+1})^2}{n},\quad(n\geq1).\]

• 공변 증명은 표준편차 증명과 똑같은 방식이므로 생략한다.

• 각 값의 가중치가 같지 않을 때의 평균과  공변 식

평균: \(\mu_{n+1} = \displaystyle\frac{W_n}{W_n+w_{n+1}}\mu_n + \displaystyle\frac{w_{n+1}a_{n+1}}{W_n+w_{n+1}},\quad(n\geq0). \)

표준편차: \(\sigma^2_{n+1}=\displaystyle\frac{W_n}{W_n+w_{n+1}}\sigma^2_{n}+\displaystyle\frac{w_{n+1}(\mu_{n+1}-a_{n+1})^2}{W_n},\quad(n\geq1).\)

공변: \(\sigma(a,b)_{n+1}=\displaystyle\frac{W_n}{W_n+w_{n+1}}\sigma(a,b)_{n}+\displaystyle\frac{w_{n+1}(\mu_{n+1}-a_{n+1})(\nu_{n+1}-b_{n+1})}{W_n},\quad(n\geq1).\)

위 식들에서 \(W_n=\displaystyle\sum^n_{i=1}w_i\)이고, \(w_n\)은 \(n\)번째 값의 가중치이다.

• Acknowledgement
  • 이정우님이 가중치가 있을 때의 계산 오류에 대해 지적해 주어 바로고침

우리말로 평균 제곱근이라고 하는 이 값은 수학에서는 다음과 같이 정의한다.

\[x_\text{rms}=\sqrt{\frac1n\sum_{i=1}^n x_i^2}\]

ROOT v5.34.04 버전 기준으로 TH1의 GetRMS()라는 메소드와, 특정 집단의 사람들이 RMS라고 말을 하면서 표준편차를 사용하는 경우가 있으니 주의하자. 참고로 GetRMS()에서 나오는 값은 표준편차이고 다음과 같이 정의된다.

\[\sigma=\sqrt{\frac1n\sum^n_{i=1}(x_i-\mu)^2}\]

물론 \(\mu\)는 분포의 평균이다. 자유도 글을 잘 생각해보면 \(n\)으로 나누는것이 아니라 \(n-1\)로 나누는게 더 좋은 값을 준다는 것을 알 수 있다.

어떤 것을 측정하는 실험을 한다고 할 때, 측정을 매우 많이 반복해서 얻을 수 있는 측정값의 분포는  가우시안(Gaussian) 분포, 이항(binomial) 분포 그리고 푸와송(Poisson) 분포로 크게 3가지가 있다.

실험을 한 후에 얻은 분포가 정말 맞는 분포인지 테스트 할 때 \(\chi^2\) 테스트를 사용한다. 측정값이 연속적인 값이든 불연속적인 값이든 상관없이 데이타를 한 빈에 5개 이상이 들어가고 빈 수가 4개 이상인 히스토그램으로 만든다. \(\chi^2\)는 다음과 같이 정의된다.

\[\chi^2=\sum^N_{k=1}\frac{(O_k-E_k)^2}{\sigma_k^2}.\]

위 식에서 \(E_k\)는 \(k\)번째 빈에 기대되는 데이타 갯수이고, \(O_k\)는 \(k\)번째 빈의 측정 데이타 갯수, \(\sigma_k\)는 테스트할 분포의 \(k\)번째 빈에서의 통계오차이다. 이렇게 구한 \(\chi^2\)값이 자유도 \(d\)근처에 있거나 작으면 기대했던 분포와 같은 결과임을 의미한다.

비교를 쉽게 하기 위해 \(\chi^2\)를 자유도 \(d\)로 나눠준 값을 reduced \(\chi^2\)라고 한다.

\[\tilde{\chi}^2=\frac{\chi^2}{d}.\]

이 때에는 \(\tilde{\chi}^2\)값이 1보다 많이 크지 않거나 작으면 우리가 예상했던 분포와 맞다는 것을 의미한다.

\(\tilde{\chi}^2\)값이 1보다 큰 경우, 우리가 실험에서 얻은 reduced \(\chi^2\)를 \(\tilde{\chi}^2_0\)라고 하면, 자유도 \(d\)에서의 \(\tilde{chi}^2\) 확률분포를 이용해서 \(\tilde{\chi}^2_0\)보다 큰 값이 나올 확률을 구하여 분포가 맞는지 틀리는지 여부를 결정한다.

\(N\)개의 데이타 포인트가 있을 때, 가장 일반적인 \(\chi^2\)의 정의는 다음과 같다.

\[\chi^2=\sum^N_{i=1}\left(\frac{y_i-f(x_i)}{\sigma_i}\right)^2.\]

\(f(x_i)\)는 \(x_i\)에서 기대되는 값, \(y_i\)는 실험값, \(\sigma_i\)는 표준편차이다.

예를들어 \(N\)개의 데이타 포인트를 가지고 이 점들을 가장 잘 설명하는 일차함수를 찾는다고 하자. 이 때 찾아야 할 미정계수는 절편값과 기울기로 2개이다.

여기서 자유도는 측정한 횟수에서 찾아야 할 미정계수를 뺀 수로 정의한다. 따라서, 위 예제의 자유도는 \(N-2\)가 된다.

(2013. 4. 6) 추가

\(N\)개의 데이타 포인트를 가지고 평균값을 구했다고 해보자. 이 때, 평균값을 가지고 있으면, \(N-1\)개의 데이타 포인트만 있으면 나머지 하나의 값은 구해낼 수 있기 때문에, 이 계의 자유도는 \(N-1\)이 된다.

자유도는 데이타 포인트의 갯수에서 제한자의 갯수를 빼준 수로 정의한다. 제한자는 데이타 포인트를 통해서 계산해낸 값을 말한다.

 

• 참고 문헌

1. John R. Taylor, An Introduction to Error Analysis 2nd Edition, University Science Books (1997).
   

최대 우도(maximum likelihood) 방법을 이용한 일련의 실험 포인트를 가장 잘 맞추는 직선을 찾아내는 해석적인 방법을 선형 회귀(linear regression) 또는 선의 최소 제곱 핏(least-squares fit for a line)이라고 한다.

 

• Linear fit

우리가 측정한 실험 데이타 포인트들이 1차 함수를 따라야 한다고 하고, 이 함수를 다음과 같이 정의하자.

\[y=ax+b\]

측정값의 \(x\)의 오차는 무시할 수 있다고 가정하고, \(y\)의 오차는 진실값에서 표준편차 \(\sigma_y\)를 가지는 가우시안 분포로 존재한다고 할 때, \(i\)번째 값 \(y_i\)를 얻을 확률은 다음과 같다.

\[P(y_i)\propto\frac1\sigma_ye^{-(y_i-ax_i-b)^2/2\sigma_y^2}\]

\(N\)개의 데이타 포인트를 얻을 확률은, 각각의 확률의 곱으로 나타낼 수 있으므로,

\[P(y_1,\cdots,y_N)\propto\frac1{\sigma_y^N}e^{-\sum_i(y_i-ax_i-b)^2/2\sigma_y^2}\]

와 같이 나타낼 수 있고, 여기서 지수함수의 지수부분의 합을 다음과 같이 \(\chi^2\)로 정의한다.

\[\chi^2=\sum^N_{i=1}\frac{(y_i-ax_i-b)^2}{\sigma^2_y}\]

미정계수 \(a\)와 \(b\)의 최적 추정값은 주어진 데이타 포인트들을 얻을 확률을 가장 크게 하는, 즉 \(\chi^2\) 값을 가장 작게 하는 방식으로 얻어낸다. \(\chi^2\) 값은 2차 함수이기 때문에 임의의 0이 아닌 임의의 값 \(a\)에 대해 최소값을 항상 가지고, 이 값은 \(a\)와  \(b\)로 \(\chi^2\) 값을 미분해서 얻어낼 수 있다.

미분값이 0이 되는 \(a\), \(b\)에 관한 두 식을 적어보면,

\[\begin{split}a\textstyle\sum x_i+bN&=\textstyle\sum y_i \\ a\textstyle\sum x_i^2+b\textstyle\sum x_i &=\textstyle\sum x_iy_i \end{split}\]

와 같고, \(a\), \(b\)의 값은 다음과 같다.

\[a=\frac{N\textstyle\sum x_iy_i-\textstyle\sum x_i\textstyle\sum y_i}{N\textstyle\sum x_i^2-(\textstyle\sum x_i)^2},\quad b=\frac{\textstyle\sum x_i^2\textstyle\sum y_i-\textstyle\sum x_i\textstyle\sum x_iy_i}{N\textstyle\sum x_i^2-(\textstyle\sum x_i)^2}.\]

이렇게 얻어진 선을 데이타의 최소 제곱 핏(least-square fit) 또는 x에 대한 y의 회귀선(line of regression)이라고 한다.

 

• 측정값 \(\mathbf{y}\)의 표준편차 \(\mathbf{\sigma}_\mathbf{y}\)

측정값 \(y\)는 참값 주위로 가우시안 분포를 가지고 있기 때문에, 최적 추정 핏으로 얻은 값 주위로 가우시안 분포를 가지고 있다고 볼 수 있다. 따라서 표준편차는 다음과 같이 구할 수 있다.

\[\sigma_y=\sqrt{\frac1N\sum(y_i-ax_i-b)^2}.\]

하지만, 데이타 포인트들을 가지고 미정계수 2개(\(a\), \(b\))를 구해냈기 때문에, \(N\)자리에는 이 계산에서의 자유도 \(N-2\)를 넣어 다음과 같이 표준편차를 구한다.

\[\sigma_y=\sqrt{\frac1{N-2}\sum(y_i-ax_i-b)^2}.\]

 

• 미정계수 \(\mathbf{a}\)와 \(\mathbf{b}\)의 오차

\[a=\frac{N\textstyle\sum x_iy_i-\textstyle\sum x_i\textstyle\sum y_i}{N\textstyle\sum x_i^2-(\textstyle\sum x_i)^2},\quad b=\frac{\textstyle\sum x_i^2\textstyle\sum y_i-\textstyle\sum x_i\textstyle\sum x_iy_i}{N\textstyle\sum x_i^2-(\textstyle\sum x_i)^2}.\]

위 식에서, 각각을 오차 전파식을 사용해서 구해내면 된다. 다음과 같다.

\[\sigma_a=\sigma_y\sqrt{\frac N{N\textstyle\sum x_i^2-(\textstyle\sum x_i)^2}},\quad\sigma_b=\sigma_y\sqrt{\frac {\textstyle\sum x_i^2}{N\textstyle\sum x_i^2-(\textstyle\sum x_i)^2}}\]

 

• 측정값 \(y\)의 각각의 표준편차가 다를 경우

가중 최소 제곱(weighted least-squares) 방법을 사용해야 한다. 별 다른건 없고, \(\chi^2\)의 \(\sigma_y\)를 각 항에 맞는 표준편차인 \(\sigma_i\)로 바꿔준 후 계산하면 된다.

 

• 실제 적용을 위한 Tip

위에서 얻은 두 식을 행렬로 적어보면,

\[\left[\begin{matrix}\textstyle\sum x_i & N \\ \textstyle\sum x_i^2 & \textstyle\sum x_i \end{matrix}\right]\left[\begin{matrix} a \\ b \end{matrix}\right]=\left[\begin{matrix}\textstyle\sum y_i \\ \textstyle\sum x_iy_i\end{matrix}\right]\]

와 같이 적을 수 있고, \(a\)와 \(b\)는 내부 함수로

\[\left[\begin{matrix} a \\ b \end{matrix}\right]=\left[\begin{matrix}\textstyle\sum x_i & N \\ \textstyle\sum x_i^2 & \textstyle\sum x_i \end{matrix}\right]^{-1}\left[\begin{matrix}\textstyle\sum y_i \\ \textstyle\sum x_iy_i\end{matrix}\right]\]

을 계산하면 쉽게 얻을 수 있다.

(2013. 4. 5) 라고 생각했지만 오차까지 계산할거면 행렬보다는 단순합이 더 나은 것 같다.

 

• 의문

1. 오차가 임의적일 때 측정값 \(y\)의 표준편차는 어떻게 구할 것인가?
   1. 오차가 임의적이면 실험을 잘못한거다. 어떤 값이 1차함수로 피팅돼야 한다면, 참값주위에 측정값이 많이 분포할 수 밖에 없다.

 

• 참고문헌

1. John R. Taylor, An Introduction to Error Analysis 2nd Edition, University Science Books (1997).