무엇이 더 정규분포인가? Jarque-Bera Test

 

 

이번 글은 

2022.06.17 - [수학의 재미/확률분포] - 정균분포 난수를 만들어 주세요 #6 : random.normal

정균분포 난수를 만들어 주세요 #6 : random.normal

에서 이어집니다.  저번글까지 normal random number를 만드는 방법에 대해서 다뤘습니다. 이론적인 방법 3개에 python 함수까지 다뤘었죠.

방법1	중심극한 정리를 이용하여, uniform random number의 표본평균 구하기  [여기]
방법2	Normal cumulative distribution function 을 이용하기 [여기]에서 [여기]까지
방법3	Box Muller 방법으로 정규분포 난수 2개를 한방에 구하기 [여기]
방법4	Pythom 함수: numpy.random.normal [여기]

이제 이러한 방법들이 정규분포 난수를 안정적으로 만들어 내는지, 더 나아가 어떤 방법이 가장 정규분포를 따르는 난수를 만들어 내는지 통계학적으로 분석해 보도록 하겠습니다. 이른바 정규성 검정입니다.

무슨 방법이 쵝오일까요~?

 

복잡한 이론은 일단 제외하고 어떤 통계량 하나 소개하겠습니다. Jarque-Bera test 통계량이라고, 줄여서 JB라 부르겠습니다.

$n$개의 통계데이터 $x_1,x_2,\cdots, x_n$에 대하여
$$ {\rm{JB}}~=~\frac{n}6 \cdot \left( S^2 + \frac{(K-3)^2}{4}\right) $$
라 정의한다. 여기서,
$$S=\frac{\hat{\mu_3}}{\hat{\mu_2}^{3/2}}~,~ K=\frac{\hat{\mu_4}}{{\hat{\mu_2}}^2}$$
이고
$$\hat{\mu_j} =\frac1n \sum_{i=1}^n (x_i-\bar{x})^j $$
이다. $\bar{x}$는 $n$개 데이터의 표본평균이다.

이럴 때, Jarque-Bera test는 다음과 같습니다(우선 결과만 알고 넘어갑시다.)

Jarque-Bera test

귀무가설 $H_0$ 를 
                                             $H_0$: $x_1,\cdots,x_n$이 표준정규분포를 따른다.

라 할 때, 어떤 유의 수준 $\alpha$에 대하여 
$$JB\geq \chi^2_{1-\alpha, 2}$$
라면 $H_0$를 기각한다.

만일 유의 수준을 일반적인 0.05라 한다면, $\chi^2_{1-\alpha,2} = 5.99$ 이므로, 

            " 어떤 데이터 ${x_i}$들의 JB가 5.99를 넘는다면, 이 데이터들은 정규분포 데이터라 볼 수 없다."
라는 뜻이다.

입니다. 만일 유의 수준을 0.01이라 한다면, $\chi^2_{1-\alpha,2} = 9.21$ 이 되겠습니다.

 

자 그럼 테스트를 진행해보겠습니다. 다음의 5가지 경우를 분석해보겠습니다.

번외	균등분포 난수의 JB 구하기 numpy.random.uniform 사용, JB가 5.99는 넘을 것이라 예상됨
방법1	python 의 numpy.random.normal 사용
방법2	central limit theorem을 이용한 방법
방법3	normal cdf를 이용한 방법
방법4	Box-Muller 방법

번외는 비교 분석을 위해 추가한 테스트입니다.

 

먼지 JB 구하는 python code입니다.

def JB(data):
    avr = np.average(data)
    n = len(data)
    mu2 = np.average((data - avr) ** 2)
    mu3 = np.average((data - avr) ** 3)
    mu4 = np.average((data - avr) ** 4)
    skew = mu3 / mu2 ** (3 / 2)
    kurt = mu4 / mu2 ** 2
    JB = n / 6 * (skew ** 2 + (kurt - 3) ** 2 / 4)
    return JB

수식에 충실한 식입니다.

 

그다음 테스트 code입니다.

import math
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
import scipy.stats as stats

def Jarque_Bera_test():
    nSize = 100000

    # uniform 난수를 발생시켜 JB test를 진행한다.
    data_uniform = np.random.rand(nSize)
    jb_uniform = JB(data_uniform)

    # python이 제공하는 random.normal 함수로 정규분포난수를 발생시킨다.
    dataSet1 = np.random.normal(0, 1, nSize)
    jb1 = JB(dataSet1)

    # Central limit theorem으로 정규분포 난수를 발생시킨다.
    dataSet2 = []
    for i in range(nSize):
        rand_unif = np.random.uniform(low=0, high=1, size=12)
        normal_data = np.sum(rand_unif) - 6
        dataSet2.append(normal_data)
    jb2 = JB(dataSet2)

    # 정규분포 cdf의 역함수를 이용하여 정규분포난수를 발생시킨다.
    rand_unif = np.random.uniform(low=0, high=1, size=nSize)
    dataSet3 = norm.ppf(rand_unif)
    jb3 = JB(dataSet3)

    # Box-Muller method로 정규분포난수를 발생시킨다.
    dataSet4 = box_muller(nSize)[0]
    jb4 = JB(dataSet4)

    print('-' * 50)
    print('번외 : uniform random number     JB : {:.4f}'.format(jb_uniform))
    print('방법1: random.normal function    JB : {:4f}'.format(jb1))
    print('방법2: central limit theorem     JB : {:4f}'.format(jb2))
    print('방법3: normal  cdf               JB : {:4f}'.format(jb3))
    print('방법4: Box-Muller method         JB : {:4f}'.format(jb4))
    print('-' * 50)

    print('Using python scipy.stats.jarque_bera function')
    res_uniform = stats.jarque_bera(data_uniform)
    res1 = stats.jarque_bera(dataSet1)
    res2 = stats.jarque_bera(dataSet2)
    res3 = stats.jarque_bera(dataSet3)
    res4 = stats.jarque_bera(dataSet4)

    print('번외의 결과', res_uniform)
    print('방법1의 결과', res1)
    print('방법2의 결과', res2)
    print('방법3의 결과', res3)
    print('방법4의 결과', res4)


def box_muller(n=1):
    u = np.random.rand(n)
    v = np.random.rand(n)

    r = np.sqrt(-2 * np.log(u))
    x = np.cos(2 * np.pi * v)
    y = np.sin(2 * np.pi * v)
    z1 = r * x
    z2 = r * y

    return z1, z2

if __name__ == '__main__':
    Jarque_Bera_test()

 

정규분포 난수를 생성하는 방법은 이전 글들에서 모두 다뤘으니 설명 없이 넘어가겠습니다. 그런데 code를 보니 아주 간편한 부분이 있습니다. 바로

 

    res_uniform = stats.jarque_bera(data_uniform)
    res1 = stats.jarque_bera(dataSet1)
    res2 = stats.jarque_bera(dataSet2)
    res3 = stats.jarque_bera(dataSet3)
    res4 = stats.jarque_bera(dataSet4)

이 부분이죠! scipy.stats라는 모듈은 jarque_bera라는 함수를 제공합니다.  간단히 설명을 하면

 

scipy.stats.jarque_bera(x)
1. Parameters
    - xarray_likeObservations of a random variable.
 2. Returns
    - jb_value (float): The test statistic
    - p (float) : The p-value for the hypothesis test.

어떤 data를 입력받아서 JB값과 p-value를 한 번에 제공하는 함수입니다. 이것만 있으면 저처럼 JB함수를 새롭게 coding 하거나 그럴 필요가 없겠죠. 어쨌든 결과를 보겠습니다.

 

--------------------------------------------------
번외 : uniform random number     JB : 5986.5002
방법1: random.normal function    JB : 1.868369
방법2: central limit theorem     JB : 28.687458
방법3: normal  cdf               JB : 0.910101
방법4: Box-Muller method         JB : 0.665952
--------------------------------------------------
Using python scipy.stats.jarque_bera function
번외의 결과 Jarque_beraResult(statistic=5986.500173417081, pvalue=0.0)
방법1의 결과 Jarque_beraResult(statistic=1.868368929818959, pvalue=0.39290616367033504)
방법2의 결과 Jarque_beraResult(statistic=28.6874578358127, pvalue=5.896545866024283e-07)
방법3의 결과 Jarque_beraResult(statistic=0.9101005870688575, pvalue=0.6344160601198944)
방법4의 결과 Jarque_beraResult(statistic=0.6659521848710391, pvalue=0.7167873305901779)

Process finished with exit code 0

역시 번외인 uniform random number는 JB값이 엄청 큽니다. p-value도 아예 0이네요.

 

반면에 방법 1~방법 4는 아주 안정적인 값을 제공하고 있습니다. 기준점인 5.99와 비교할 수 없을정도로 0에 가까운 값들이죠. 귀무가설을 채택 안할 만한 근거가 없다, 즉 정규 분포가 아니라고 결론낼 만한 근거가 없다는 뜻입니다. 말이 좀 복잡하지만, 어쨋든 통계적으로 봤을 때, 정규분포를 따르는 것 같다라는 얘기입니다.

 

저희가 지금까지 알아본 정규분포 난수를 구하는 모든 방법이 아주 훌륭한 방법들이라는 것을 Jarque-Bera test를 통해 알았습니다. 일반적으로 central limit theorem을 쓴 방법이 조금 효과가 떨어지고 Box Muller 방법이 좋네요. 

아울러 직접 만든 방법과 stats의 jarque-bera 함수를 이용한 결과가 똑같이 나오는 것도 확인할 수 있습니다.

 

앞으로의 글에서 정규분포 난수는 그냥 python numpy의 random.normal 함수를 이용해서 편리하게 얻어내도 될 것 같습니다.

균등난수테스트 빼고 다 좋은 결과임