[확률및통계] 모멘트

우리가 통계를 배울 때 가장 먼저 배우는 개념은 평균과 분산(또는 표준편차)이다.

 

이 2개는 어떤 분포를 표현하는데에 사용되는 지표라고도 이해할 수 있다.

가우시안 분포의 경우에는 평균과 분산 두 개만 가지고도 표현할 수 있는 분포이기 때문에 많이 쓰이지만 실제 분포가 항상 가우시안일 것이라고 기대하기 힘들다.

 

가우시안 분포를 많이 쓰는 이유는 central limit theorem 때문이다.

개별 표본이 같은 분포(꼭 가우시안이 아니어도 됨)에서 추출되고 서로 독립일 때 이 표본들의 평균이 표본 수 증가함에 따라 가우시안 분포를 따를 것이라는 Theorem이다.

 

하여튼 이런 경우가 아니라면 분포를 표현하기 위해 여러 moment를 공부하게 된다.

 

 

$x$는 continuous random variable 이고

아래에서 다룰 내용은 $x$의 샘플로부터 구하는 Mean, Variance, Skew, Kurtosis이다. 

 

1. Mean

 

샘플 N개가 있을 때 일반적으로 평균을 구하는 방법은 다음과 같다.

$\operatorname{Mean}(x)=\mu_{x}=\mathbb{E}(x)=\dfrac{1}{N}\sum_{i=1}^{N}x_{i}$

 

확률 분포를 아는 경우에는 다음과 같이 구한다.

 

$\mu_{x}=\int_{x=-\infty}^{x=\infty} x P(x) dx$

 

평균을 구하는 것 자체가 확률 분포를 곱해서 적분한 것이므로 아래에서 나오는 $\mathbb{E}$는 모두 확률분포가 존재할 경우에 위와 같이 구한다고 보면 된다.

 

2. Variance

 

$\operatorname{Var}(x)=\mathbb[(x-\mu)^2]$

 

확률분포가 아니라 샘플로부터 분산을 구할 때는 아래와 구할 것처럼 생각할 수도 있다.

$\operatorname{Var}(x)=\dfrac{1}{N}\sum_{i=1}^{N}\mathbb[(x-\mu)^2]$

 

그러나 위의 식은 N이 큰 경우에만 쓸 만한 것이. variance 추정을 위한 unbiased estimator는

$\operatorname{Var}(x)=\dfrac{1}{N-1}\sum_{i=1}^{N}\mathbb[(x-\mu)^2]$

이므로 이 값을 사용한다.

 

$\sqrt(\operatorname{Var}(x))=\operatorname{Standard Deviation}$

 

3. Skew

Third standarized moment로 불린다.

 

$\operatorname{Skew}=\gamma_1:=\tilde{\mu}_3=\mathbb{E}\left[\left(\frac{X-\mu}{\sigma}\right)^3\right]=\frac{\mu_3}{\sigma^3}=\frac{\mathbb{E}\left[(X-\mu)^3\right]}{\left(\mathbb{E}\left[(X-\mu)^2\right]\right)^{3 / 2}}=\frac{\kappa_3}{\kappa_2^{3 / 2}}$

 

샘플에 대해서 구할 때는 다음과 같다.

 

$b_1=\frac{m_3}{s^3}=\frac{\frac{1}{n} \sum_{i=1}^n\left(x_i-\bar{x}\right)^3}{\left[\frac{1}{n-1} \sum_{i=1}^n\left(x_i-\bar{x}\right)^2\right]^{3 / 2}}$

 

$g_1=\frac{m_3}{s^3}=\frac{\frac{1}{n} \sum_{i=1}^n\left(x_i-\bar{x}\right)^3}{\left[\frac{1}{n} \sum_{i=1}^n\left(x_i-\bar{x}\right)^2\right]^{3 / 2}}$

 

마찬가지로 unbiased estimator로 계산하려면

$G_1=\frac{k_3}{k_2^{3 / 2}}=\frac{n^2}{(n-1)(n-2)} b_1=\frac{\sqrt{n(n-1)}}{n-2} g_1$

다음과 같이 계산하는데 우리가 알고 싶은 skewness는 $k_{3}$이다.

 

4. Kurtosis

kurtosis는 정의 상 다음과 같은데

$\operatorname{Kurt}[X]=\mathbb{E}\left[\left(\frac{X-\mu}{\sigma}\right)^4\right]=\frac{\mathbb{E}\left[(X-\mu)^4\right]}{\left(\mathrm{E}\left[(X-\mu)^2\right]\right)^2}=\frac{\mu_4}{\sigma^4}$

 

앞서 본 것처럼 n개의 샘플에 대해서 구하려면 처음에 드는 가장 간단하게 구하는 방법은 아래이다.

 

$g_2=\frac{m_4}{m_2^2}-3=\frac{\frac{1}{n} \sum_{i=1}^n\left(x_i-\bar{x}\right)^4}{\left[\frac{1}{n} \sum_{i=1}^n\left(x_i-\bar{x}\right)^2\right]^2}-3$

 

하지만 위 식은 biased estimator다. 따라서 아래와 같이 구한다.

 

$\begin{aligned}
G_2 & =\frac{k_4}{k_2^2} \\
& =\frac{n^2\left[(n+1) m_4-3(n-1) m_2^2\right]}{(n-1)(n-2)(n-3)} \frac{(n-1)^2}{n^2 m_2^2} \\
& =\frac{n-1}{(n-2)(n-3)}\left[(n+1) \frac{m_4}{m_2^2}-3(n-1)\right] \\
& =\frac{n-1}{(n-2)(n-3)}\left[(n+1) g_2+6\right] \\
& =\frac{(n+1) n(n-1)}{(n-2)(n-3)} \frac{\sum_{i=1}^n\left(x_i-\bar{x}\right)^4}{\left(\sum_{i=1}^n\left(x_i-\bar{x}\right)^2\right)^2}-3 \frac{(n-1)^2}{(n-2)(n-3)} \\
& =\frac{(n+1) n}{(n-1)(n-2)(n-3)} \frac{\sum_{i=1}^n\left(x_i-\bar{x}\right)^4}{k_2^2}-3 \frac{(n-1)^2}{(n-2)(n-3)}
\end{aligned}$

 

요약하면 

$\dfrac{k_{4}}{k_{2}^2}=\frac{n-1}{(n-2)(n-3)}\left[(n+1) g_2+6\right]$이다.

우리가 구하고자 하는 unbiased kurtosis estimator는 $k_{4}$이므로 구할 수 있다. ($k_{2}$는 위에서 구했던 unbiased estimate of sample variance이다.)

 

---

예시

 

각 지표가 실제 표본에서 어떤 의미를 갖는지 알기 위해 가우시안 분포를 갖는 랜덤 변수를 비선형 변환을 통하여 확률 분포를 바꿔보았다.

 

아래는 $x(k+1)=x(k)+ \left( f(x(k),u(k))+ g(x(k),u(k)) \right) \Delta t$ 인 시스템이고 u는 주어져 있지만 x는 주어진 값 + 가우시안 노이즈가 있다.

 

이 경우에 $x(k+1)$의 분포가 어떻게 되는지 보았다.

이렇게 한 이유는, skew와 kurtosis가 유의미한 값이 나오도록 하려면 비선형으로 변환해주는 게 편하기 때문이다.

 

MATLAB에서 skew와 kurtosis를 계산한 방식에 대해서는 다음과 같다.

 

Skewness

 

Kurtosis

 

% MATLAB CODE for visualizing the distribution
close all; clc;
rng('default') ;

x_tmp = -pi:2:pi;
u_tmp = -pi:2:pi;

f1_fun = @(x,u) sin(x)+cos(u);
f2_fun = @(x,u) cos(x)+sin(u);
M = 100000 ;

for k = 1:length(x_tmp)
    for j = 1:length(u_tmp)
        x_tf = zeros(M,1);
        x_corp = zeros(M,1);
        for i = 1:M
            x_temp = x_tmp(k) ;
            u_temp = u_tmp(j) ;

            x_corp(i) = x_temp + randn(1) ;
            dw = sqrt(Qs) * randn(1) ;
            f1 = f1_fun(x_corp(i),u_temp) ;
            f2 = f2_fun(x_corp(i),u_temp) ;
            dxdt = f1 + f2 * dw ;
            x_tf(i) = x_corp(i) + dxdt * dt ;
        end
        
        figure;
        hold on ;
        title('x='+string(x_temp)+' u='+string(u_temp));
        histogram(x_tf, 'Normalization','pdf') ;
        histogram(x_corp,'Normalization','pdf') ;
        xline(x_temp,'--r','LineWidth',2.0) ;
        dim = [.1 .1 .3 .3];
        str = "mean : "+num2str(mean_val,3)+ " var : "+num2str(var_val,3) +...
            " skew : " +num2str(skew,3) + " kurt : " + num2str(kurt,3);
        annotation('textbox',dim,'String',str,'FitBoxToText','on','BackgroundColor','w') ;
        legend('after','before','mean');
        pause(1);

        skew = skewness(x_tf) ;
        kurt = kurtosis(x_tf) ;
        mean_val = mean(x_tf) ;
        var_val = var(x_tf) ;

        fprintf("mean : %f   ", mean_val) ;
        fprintf("variance : %f   ", var_val) ;
        fprintf("skewness : %f   ", skew) ;
        fprintf("kurtosis : %f \n", kurt) ;
    end
end

 

 

(클릭하면 더 자세히 볼 수 있다)

 

아래 그림을 보면 빨간색 점선이 $x(k)$의 평균이고 주황색이 그 분포를 보여준다. 파란색이 변환 후의 분포를 말한다.

공정한 비교를 위해 정규분포의 값을 알려주면,

 

Mean	Variance	Skewness	Kurtosis
$\mu$	$\sigma^{2}$	0	3

위 표에서 Kurtosis가 정규분포일 때 3인데, 위의 설명을 보면 kurtosis 식에서 -3이 되어있는 것을 알 수 있다. 이는 정규분포가 3이므로 정규분포일 때 kurtosis가 0으로 나오게 하기 위해 사용된 정의라고 MATLAB에 설명이 나와있다.

 

여러 개의 분포를 보면 알겠지만 꼭 옆으로 치우친 분포라고 해서 skew의 절댓값이 큰 것은 아니다.

 

 

---

참고자료

위키피디아

https://kr.mathworks.com/help/stats/skewness.html

https://kr.mathworks.com/help/stats/kurtosis.html