판별식을 시뮬레이션으로?

 

2차방정식과 관련된 좀 유명한 문제가 있습니다.  다음과 같습니다.

$a,b,c$ 가 $(0,1)$ 구간에 균등하게 분포되어 있는 난수일 때, 2차 방정식
$$ax^2 + bx+ c=0$$
이 실수근(real root)을 가질 확률을 구하여라.

몸에 새길정도의 중요성을 자랑하는 근의 공식

이 문제를 uniform random number를 생성하여 풀어보도록 하겠습니다.

2차 방정식 $ax^2+bx+c=0$ 이 실근을 갖기 위해서는 이것의 판별식(discreminant)이 $0$ 이상이어야 하죠. 즉

$$b^2 -4ac \geq 0$$

을 만족해야 합니다.

따라서 코딩은 아래와 같이 하면 될 것 같습니다.

 

1. $[0,1]$ 구간의 uniform random number를 3개 생성하고 각각을 $a,b,c$에 대입한다.
2. $b^2\geq 4ac$를 만족하면 카운팅 하고, 아니면 그냥 넘어간다.
3. 1~2 과정을 아주 많이 반복한다.
4. 반복이 끝나고, 카운팅 된 숫자를 전체 반복 숫자로 나누면 그것이 바로 확률이다.

아래는 간단히 python으로 구현해 본 것입니다. 전체 코드는 아래와 같습니다.

import math
import numpy as np

def rootTest():

    np.random.seed(1)
    nSimulation = 1000000
    cnt = 0
    for i in range(nSimulation):
        rdn = np.random.uniform(low=0, high=1, size=3)
        if rdn[1]**2 >= rdn[0]*rdn[2]*4:
            cnt += 1
    prob = cnt/nSimulation
    real_value = 5/36 + np.log(2)/6
    print('Simulation probability: {:.4f}'.format(prob))
    print('Theoretic  probability: {:.4f}'.format(real_value))

if __name__ == '__main__':
    rootTest()​

중요한 부분만 살펴보면

 

 

    nSimulation = 1000000

전체 반복 횟수는 100만 번으로 설정합니다.

 

        rdn = np.random.uniform(low=0, high=1, size=3)​

numpy 모듈의 random.uniform 을 사용합니다. $[0,1]$사이의 uniform number로서 size =3 개를 한 번에 뽑습니다.

$a,b,c$는 각각 rdn[0], rdn[1], rdn[2] 라 생각합니다.

 

        if rdn[1]**2 >= rdn[0]*rdn[2]*4:
            cnt += 1

판별식을 구하여 $b^2\geq 4ac$를 만족하면 cnt 변수를 증가시킵니다.

 

    prob = cnt/nSimulation

반복이 끝난 후 카운팅 변수인 cnt을 전체 반복 횟수로 나누어 확률을 구합니다.

 

    real_value = 5/36 + np.log(2)/6

사실, 이 문제의 정확한 답이 있습니다. 그것을 계산하여 prob와 비교해 볼 참입니다. 이제 결과를 보시죠.

 

---

Simulation probability: 0.2545
Theoretic  probability: 0.2544

Process finished with exit code 0

어떻습니까? 거의 정확한 답과 시뮬레이션을 통한 추정값이 비슷하죠? 이러한 시뮬레이션을 MonteCarlo Simulation이라 합니다. MonteCarlo Simulation에 대해서는 나중에 상당한 비중의 글을 통해 다뤄볼 생각입니다.

---

그렇다면 이 문제의 정확한 답은 어떻게 구할까요? 

$(a,b,c)$의 joint probability density function $f(a,b,c) =1$입니다. 따라서 영역 $\mathcal{R}$을

$$\mathcal{R}=\{(a,b,c)\in\mathbb{R}^3 | b^2\geq 4ac \}$$ 라 정의한다면 우리가 구하는 확률은

$$\int_{\mathcal{R}} f(a,b,c) dadbdc = \int_{\mathcal{R}} dadbdc\tag{1}$$

입니다. 영역 $\mathcal{R}$을 좀 더 알기 쉽게 써 보면,

$$\mathcal{R} = \{(a,b,c)\in\mathbb{R}| b\geq 2\sqrt{ac} ~,~ c\leq \min\left( \textstyle\frac{1}{4a},1\right) ~ , ~0\leq a\leq 1 \}$$

왜일까요? $b\geq2\sqrt{ac}$는 판별식에서 나왔고, $0\leq b\leq 1$이므로 $1\geq 2\sqrt{ac}$, 즉, $\frac1{4a} \geq c$ 인 것입니다. 그런데 $c$는 $c\leq 1$을 만족하므로 결국 $c\leq \min(\frac{1}{4a},1)$을 만족하는 것이죠.

현기증이 좀 나네여

자 이제 적분만 하면 끝입니다. 수식(1)은

$$\begin{align} \int_0^1 &  \int_0^{\min(\frac{1}{4a},1)} \int_{2\sqrt{ac}}^1 dbdcda \\ & = \int_0^1 \int_0^{\min(\frac{1}{4a},1)} (1-2\sqrt{ac}) dc da \\ & =\int_0^1 \left[ c- \frac43 \sqrt{a} c^{3/2} \right]_0^{\min(\frac{1}{4a},1)} da \\ & = \int_0^1 \left( \min(\frac{1}{4a},1) - \frac43\sqrt{a} \min(\frac{1}{4a},1)^{3/2} \right) da \\ & = \int_0^{\frac14} \left(1-\frac43\sqrt{a}\right)da + \int_{\frac14}^1 \left( \frac{1}{4a} - \frac43\sqrt{a}\left(\frac1{4a}\right)^{3/2}\right) da\\ & = \left[ a- \frac89 a^{3/2}\right]_0^{\frac14} + \int_{\frac14}^1 \frac{1}{12a} da \\ & = \frac5{36} + \frac1{12} \left[\ln a\right]_{\frac14}^1 =\frac5{36}+\frac16\ln 2 \end{align}$$

Do Math!

 

별거 아니죠? [영화 '이상한 나라의 수학자' 캡쳐, 제공=쇼박스]