e 를 시뮬레이션으로 구하기

 

테일러 전개(여기를 참고)에 따르면, 오일러 수(Euler's number)라고도 불리는 $e$는

$$e = 1+ 1+\frac1{2!} + \frac1{3!} +\cdots$$

로 구할 수 있다고 하였습니다.

 

Euler's number Kids T-Shirt [Designed and sold by Stargazerr] 얼마나 중요한 숫자인지 어린이 옷에 새겨넣었음~

 

$e$의 근사값을 구하는 방법은

2022.05.19 - [수학의 재미/아름다운 이론] - 테일러 전개 #2 : $e$에 대하여 에서 소개를 한 바 있는데요.

 

$e$의 값을 구하는 아주 흥미로운 방법이 있어서 소개해 보려 합니다.

 

$U_1, U_2, U_3, \cdots$를 $\mathcal{U}[0,1]$ 분포를 따르는 iid 확률변수라고 합시다. 즉, 0과 1에서 균등 분포를 따르는 서로 독립인 확률 변수입니다.
이때, $N$을
$$N = \min \{ n| U_1+U_2+\cdots +U_n >1 \}$$
인 확률 변수라 정의합시다. 즉 $N$은 균등 분포 변수를 더해 나가면서 최초로 1이 될 때 더해진 횟수를 뜻합니다. 그러면
$$\mathbb{E}(N) = e$$
가 됩니다.

 

증명을 위해 균등 분포의 합이 어떤 분포를 따르는지 떠올려 봅시다.

2022.05.25 - [수학의 재미/확률분포] - Irwin-Hall Distribution

Irwin-Hall Distribution

균등분포의 합이 따르는 분포를 Irwin-Hall distribution이라 합니다. cdf를 복습해보면

 

$F_X(t;n) = \mathbb{P}(U_1+U_2+\cdots +U_n \leq t)$ 는

 

$$F_X(t;n) = \frac{1}{n!}\sum_{i=0}^{\lfloor t \rfloor}  (-1)^i {n \choose i} (t-i)^n$$

 

입니다 여기서 $\lfloor t \rfloor$는 $t$를 넘지 않는 최대의 정수를 뜻합니다. 그리고 $x$는 $[0,n]$사이의 수입니다. 따라서

 

$$\mathbb{P}(U_1+U_2+\cdots +U_n \leq 1) = \frac{1}{n!}\sum_{i=0}^1 (-1)^i {n\choose i}(1-i)^n =\frac1{n!}$$

이죠.

 

 

 

$n$이 $U_1+U_2+\cdots+U_n>1$인 최소의 $n$이므로

$$\mathbb{P}(N=n ) = \mathbb{P}\left(\sum_{i=1}^n U_i >1 ~\cap~ \sum_{i=1}^{n-1}U_i \leq 1 \right)$$

입니다.

그런데

$$\begin{align} 1 & = \mathbb{P} \left(\sum^n U_i \leq 1 ~\cup~ \sum^n U_i  >1 \right)\\   & = \mathbb{P} \left(\sum^n U_i \leq 1\right)  + \mathbb{P} \left(\sum^n U_i >1 \right)\tag{1}\\ \end{align}$$

 

또한

$$\begin{align} \mathbb{P} \left(\sum^n U_i >1 \right) & = \mathbb{P} \left(\sum^n U_i > 1 ~\cap~  \sum^{n-1} U_i  \leq 1 \right) +\mathbb{P} \left(\sum^n U_i > 1 ~\cap~  \sum^{n-1} U_i  > 1 \right)\\ &= \mathbb{P} \left(\sum^n U_i > 1 ~\cap~  \sum^{n-1} U_i  \leq 1 \right) +\mathbb{P} \left( \sum^{n-1} U_i  > 1 \right)\\ &= \mathbb{P}(N=n)+\mathbb{P} \left( \sum^{n-1} U_i  > 1 \right)\tag{2} \end{align}$$

이므로 수식(1)과 (2)를 합하면

$$\begin{align} \mathbb{P}(N=n) & = 1-  \mathbb{P} \left(\sum^n U_i \leq 1\right) -\mathbb{P} \left( \sum^{n-1} U_i  > 1 \right)\\           &= \mathbb{P} \left(\sum^{n-1} U_i \leq 1\right) -\mathbb{P} \left(\sum^n U_i \leq 1\right) \\           & = \frac1{(n-1)^!} - \frac1{n!} \end{align}$$

이 됩니다.

 

결론적으로,

$$\mathbb{E}(N) = \sum_{n=2}^\infty n\cdot \mathbb{P}(N=n) = \sum_{n=1}^\infty n \left( \frac1{(n-1)^!} - \frac1{n!} \right)=\sum_{n=0}^\infty \frac1{n!} = e$$

를 만족합니다.

 

---

python code로 위의 정리를 확인해볼까요?

 

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

def value_of_e():
    res=[]
    nSimulation = 1000000

    for _ in range(nSimulation):
        rn = np.random.uniform(0,1,100)
        tmpSum = 0
        for i in range(100):
            tmpSum += rn[i]
            if tmpSum>1:
                break

        min_n = i+1
        res.append(min_n)

    print('average of n : {:.4f}'. format(np.average(res)))
    print('theoretic val: {:.4f}'. format(np.exp(1)))

    bins = np.zeros((20,3))
    for i in range(2,10):
        bins[i-2][0] = i
        bins[i-2][1] = round(res.count(i) /nSimulation, 4)
        bins[i-2][2] = 1/math.factorial(i-1) - 1/math.factorial(i)

	np.set_printoptions(suppress=True)    
    print(bins)

if __name__ == '__main__':
    value_of_e()

 

code를 간략하게 살펴볼까요?

 

    res=[]
    nSimulation = 1000000

○ res는 균등 분포의 합이 최초로 1이 넘었을 때 합한 개수를 저장하는 변수

○ 시뮬레이션 횟수는 100만 번 

 

 

    for _ in range(nSimulation):
        rn = np.random.uniform(0,1,100)
        tmpSum = 0
        for i in range(100):
            tmpSum += rn[i]
            if tmpSum>1:
                break

○ np.random.uniform 를 이용하여 0과 1 사이의 uniform random number를 100개 만든다.

○ 100개를 앞에서부터 차례로 더해가다가 1이 넘으면

○ break로 탈출

 

 

        min_n = i+1
        res.append(min_n)

○ i는 0부터 시작했으므로 +1을 해야 최소의 n이 됨

○ 결과를 res에 저장

 

 

    print('average of n : {:.4f}'. format(np.average(res)))
    print('theoretic val: {:.4f}'. format(np.exp(1)))

○ res의 기댓값을 np.average를 사용하여 구함

○ 이론값인 $e^1$과 비교 출력

 

 

    bins = np.zeros((20,3))
    for i in range(2,10):
        bins[i-2][0] = i
        bins[i-2][1] = round(res.count(i) /nSimulation, 4)
        bins[i-2][2] = 1/math.factorial(i-1) - 1/math.factorial(i)

○ $20\times 3$ 배열을 만들고

○ 첫째 열은 최소 횟수를 뜻하는 $n$,

○ 둘째 열은 res안의 최소 횟수 중복된 것을 카운팅 해서 시뮬레이션 횟수로 나눈, 즉, 시뮬레이션에 따른 확률을 저장

○ 셋째 열은 이론의 확률인 $\frac{1}{(n-1)!} - \frac{1}{n!}$ 을 계산하여 비교

 

 

다음은 결과입니다.

average of n : 2.7188
theoretic val: 2.7183
[[2.         0.5        0.5       ]
 [3.         0.3328     0.33333333]
 [4.         0.1256     0.125     ]
 [5.         0.0333     0.03333333]
 [6.         0.0069     0.00694444]
 [7.         0.0012     0.00119048]
 [8.         0.0002     0.00017361]
 [9.         0.         0.00002205]
 [0.         0.         0.        ]
 [0.         0.         0.        ]
 [0.         0.         0.        ]
 [0.         0.         0.        ]
 [0.         0.         0.        ]
 [0.         0.         0.        ]
 [0.         0.         0.        ]
 [0.         0.         0.        ]
 [0.         0.         0.        ]
 [0.         0.         0.        ]
 [0.         0.         0.        ]
 [0.         0.         0.        ]]

Process finished with exit code 0

어떻습니까? $e$의 근삿값이 잘 표현되죠? 그리고 시뮬레이션으로 구한 확률과 이론적으로 구한 확률도 비슷합니다.

 

이처럼 수많은 상황을 반복 발생시키고 그것의 기댓값을 구해  값의 근사치나 확률 등의 값을 알아내는 방법을

 

MonteCarlo Simulation(몬테카를로 시뮬레이션)

 

이라 합니다. 보통 $\pi$의 근사치를 구해보는 연습을 하는데요. 특이하게 $e$의 근사치를 MonteCarlo Simulation으로 할 수 있다기에 얼른 소개해 보았습니다. 내친김에 다음 글에서는 $\pi$를 구해보도록 하죠.