정균분포 난수를 만들어 주세요 #1 : 중심극한정리 이용

 

 

이번 글은 

2022.06.10 - [수학의 재미/확률분포] - 균등난수(uniform random number) 만들어보자 #3 : Halton Sequence

균등난수(uniform random number) 만들어보자 #3 : Halton Sequence

에서 이어집니다. 균등난수를 만드는 방법을 총 3회에 걸쳐서 소개했는데요. 사실 

정규분포 난수를 만들기 위해 균등분포 먼저 소개한 것

에 지나지 않습니다. 정규분포를 따르는 난수를 만드는데 균등분포가 반드시 필요하거든요.

한국영화계의 필수요소들

 

그럼 어떠한 원리로 균등분포(uniform random number)에서 정규분포 난수를 뽑는지 알아보겠습니다.

자자. 얼른 정규분포처럼 서봐 다들~[출처]위키피디아

 

잘 알려진 방법이 크게 세 가지가 있습니다.

방법1	중심극한 정리를 이용하여, uniform random number의 표본평균 구하기
방법2	Normal cumulative distribution function 을 이용하기
방법3	Box Muller 방법으로 정규분포 난수 2개를 한방에 구하기

 

방법 1. 중심 극한 정리 이용

서로 독립인 unifom random number $n$개를 준비합시다. uniform의 약자인 $u$를 써서

$u_1, u_2 , \cdots, u_n$ 이라 표시합니다. 

$$ u_i \sim \mathcal{U}[0,1]~,~ i=1,\cdots, n$$ 

이고 각각의 평균은 $\mu= \frac12$, 분산은 $\sigma^2= \frac1{12}$이므로 중심극한정리(링크참조)에 의하여

$n$이 커질수록

$$\bar{u}= \frac1n \sum_{i=1}^n u_i \sim \mathcal{N}(\mu, \frac{\sigma^2}n) =\mathcal{N}\left(\frac12, \frac1{12n}\right)$$

인 정규분포로 (분포) 수렴하게 됩니다. 이것을 표준화해보면

$$\frac{\bar{u}-\frac12}{\frac{1}{\sqrt{12n}}} \sim \mathcal{N}(0,1) $$

계산을 쉽게 하기 위해 

$$ n= 12$$ 

를 대입하여 정리하면

$$ \sum_{i=1}^{12} u_i - 6 \sim \mathcal{N}(0,1) \tag{1}$$

입니다.

 

python으로 확인해 보겠습니다.  확인할 사항은 수식(1)으로 뽑은 데이터들이

• 평균은 0, 분산은 1을 만족하는가?
• 데이터들을 histogram으로 그렸을 때, 정규분포 형태가 나오는가?

입니다. code를 보시죠.

 

import math
import numpy as np
import matplotlib.pyplot as plt

def normal_random_test_clt():
    nSize = 12
    nSimulation = 5000
    data = []
    for i in range(nSimulation):
        rand_unif = np.random.uniform(low=0, high=1, size=nSize)
        normal_data = np.sum(rand_unif) - 6
        data.append(normal_data)
    print('Average of data : {:.2f}'.format(np.average(data)))
    print('Variance of data : {:.2f}'.format(np.var(data)))
    plt.hist(data, bins=20)
    plt.show()

if __name__ == '__main__':
    normal_random_test_clt()

code 자체는 중심극한 정리에서 다룬 동전던지기, 주사위던지기와 비슷합니다. 

 

        rand_unif = np.random.uniform(low=0, high=1, size=nSize)

[0,1] 사이의 size가 nSize =12개인 난수 배열을 만들어 rand_unif에 할당합니다.

 

        normal_data = np.sum(rand_unif) - 6

수식(1)처럼 배열 내 원소의 합을 구해 6을 뺴줘서 정규분포를 따를만한 데이터를 만듭니다.

 

이제 결과를 보시죠.

AVerage of data : -0.01
Variance of data : 1.01

Process finished with exit code 0

평균은 0에 가깝고, 분산은 1에 가깝습니다. 또한 histogram은 다음과 같습니다.

 

정규분포처럼 평균 0을 중심으로 대칭 형태의 모양이 나오는 듯합니다.

---

그럼, 실제로 수식(1)의 난수가 어떤 distibutipn을 가지는 지 알아보도록 합니다.

$U_i, (1\leq i\leq 12) \sim \mathcal{U}[0,1]$ 을 i.i.d 확률변수라 하고

$$ Z = \sum_{i=1}^{12} U_i -6 $$

이라 정의합니다. 이때 $Z$의 cumulative density function $F$는

$$F(z) = \mathbb{P}(Z\leq z) = \mathbb{P}\left(\sum_{i=1}^{12} U_i \leq z+6 \right)\tag{2}$$

입니다.

그런데, $\sum_{i=1}^{12} U_i $ 는 바로 Irwin-Hall distribution을 따릅니다. Irwin-Hall cdf를 $F_{IH}$이라 쓴다면 수식(2)는

$$F(z) = F_{IH}(z+6)$$ 

이라 쓸 수 있죠.

 

우선 Irwin-Hall cdf를 coding 합시다.

def Irwin_Hall_distribution(n, x):
    sum_cnt = math.floor(x)
    temp_val = 0
    for i in range(sum_cnt):
        temp_val += (-1) ** i * math.comb(n, i) * (x - i) ** n

    return temp_val / math.factorial(n)

Irwin-Hall distribution에 있는 수식 그대로 옮겨 적은 것입니다.

 

    sum_cnt = math.floor(x)

floor(x)는 수식 기호로 $\lfloor x \rfloor$ 로, $x$를 넘지 않는 최대의 정수를 뜻합니다.

 

        temp_val += (-1) ** i * math.comb(n, i) * (x - i) ** n

math.comb(n,i)는 수식 기호로 $_{n}{\rm{C}}_i $을 뜻합니다.

 

---

이제 Irwin-Hall 함수와 normal cdf 함수를 비교해볼까요?

 

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

def Irwin_Hall_distribution(n, x):
    sum_cnt = math.floor(x)
    temp_val = 0
    for i in range(sum_cnt):
        temp_val += (-1) ** i * math.comb(n, i) * (x - i) ** n
    return temp_val / math.factorial(n)

if __name__ == '__main__':
    # normal_random_test_clt()

    x = np.arange(-5, 5, 0.1)
    x_shift = x + 6
    y_IH = []
    y_norm = []
    norm_mu, norm_std = 0, 1
    for i in range(len(x)):
        y_IH.append(Irwin_Hall_distribution(12, x_shift[i]))
        y_norm.append(norm.cdf(x[i], norm_mu, norm_std))

    y_diff = [abs(x - y) for x, y in zip(y_IH, y_norm)]
    error = np.max(y_diff)

    plt.plot(x, y_IH, 'b', label='Irwin-Hall')
    plt.plot(x, y_norm, 'g', label='N(0,1)')
    plt.legend()

    plt.show()
    print('Maximum difference between Norm.cdf and IR.cdf is {:.4f}'.format(error))

normal cdf 함수는 scipy.stats 모듈 안에 있습니다. 따라서

from scipy.stats import norm

이와 같이 import는 해 줍니다.

 

    x = np.arange(-5, 5, 0.1)

$x$를 $[-5,5]$ 구간에 그리고 싶은데 수식(2)에 따르면 $F_{IH}$에는 $x$를 $+6$만큼 shift 해서 비교해야 합니다.

    x_shift = x + 6

이렇게 하면 배열의 원소가 각각 +6이 됩니다.

 

    for i in range(len(x)):
        y_IH.append(Irwin_Hall_distribution(12, x_shift[i]))
        y_norm.append(norm.cdf(x[i], norm_mu, norm_std))​

각각의 함숫값을 계산합니다. normal cdf는 norm.cdf 함수를 써서 계산합니다.

 

    y_diff = [abs(x - y) for x, y in zip(y_IH, y_norm)]
    error = np.max(y_diff)

Irwin-Hall과 normal cdf의 함숫값의 차이를 각 $x$점에서 계산하여 그중 max를 찾습니다. 이것을 함수가 벌어진 차이라고 생각해보죠. list형끼리는 원소끼리 빼는 방법이 마땅찮아 위와 같이 zip 함수를 이용하여 계산해줍니다.

 

결과는 다음과 같습니다.

얼마나 비슷한지.. 그래프가 겹쳐서 하나로 보이네요. 두 함수의 차는

 

Maximum difference between Norm.cdf and IR.cdf is 0.0023

Process finished with exit code 0

불과 0.0023 정도 차이가 납니다. 상당히 정확하죠.

 

혹시 coding이 잘못되었나 싶어 shift를 6이 아닌 5.5를 시켜보겠습니다.

coding은 잘 된 것 같네요

 

정리하면,

균등난수 12개를 합하여 6을 빼주면 이것을 정규분포를 따르는 난수라고 생각할 수 있다

 

라는 결론을 얻습니다. 다음 글에서는 정규분포의 cdf를 이용하여 정규분포 난수 만드는 방법에 대해서 알아보겠습니다.