정균분포 난수를 만들어 주세요 #3 : Normal CDF 이용

이번 글은

2022.06.13 - [수학의 재미/확률분포] - 정균 분포 난수를 만들어 주세요 #2 : Normal CDF 이용

정균분포 난수를 만들어 주세요 #2 : Normal CDF 이용

에서 이어집니다. 

저번글에서는 정규분포 cdf의 역함수를 사용하여 uniform random number에 대응한 정규분포 난수를 만들었습니다. 이 와중에 python 함수인 norm.cdf, norm.ppf 를 사용했었죠. 

norm.cdf 를 기호로 $\Phi$라 쓰므로 norm.ppf는 $\Phi^{-1}$이 됩니다. 실제로 cdf, ppf 두 함수가 역함수 관계에 있는지 볼까요?

 

역함수 관계를 알아보는 python code입니다.

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

def plot_normcdf_normppf():
    x = np.arange(-2, 2, 0.01)
    y = norm.cdf(x)
    plt.plot(x, y, label='cdf')

    x2 = np.arange(0, 1, 0.01)
    y2 = norm.ppf(x2)
    plt.plot(x2, y2, label='ppf')

    plt.plot(x, x, color='gray', linestyle='--', label='identity')
    plt.legend()
    plt.show()

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

$\Phi : (-\infty,\infty) \rightarrow (0,1)$  이지만, $\Phi$의 정의역을 $[-2,2]$ 구간으로 제한시켜 graph를 그려봅니다. 

 

 

    plt.plot(x, x, color='gray', linestyle='--', label='identity')

이 부분은 항등 함수를 그리는 부분입니다.

이제 결과를 보면,

위와 같이 항등함수를 기준으로 대칭인 두 graph 가 나오게 됩니다. 눈으로만 봐서는 실제적인 확인이 어려우므로 직접 계산을 해봅시다. 

$\Phi^{-1}(\Phi(x)) = x$$

가 맞는지 확인해 보는 것입니다.

 

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

def normcdf_ppf_test():
    x = np.arange(-5, 5, 0.01)
    y = norm.cdf(x, 0, 1)
    z = norm.ppf(y)
    data_diff = [abs(a - b) for a, b in zip(x, z)]
    print('Maximum difference is {:4,f}'.format(np.max(data_diff)))
    plt.plot(x, z)
    plt.show()

if __name__ == '__main__':
    normcdf_ppf_test()

$x$를 $[-5,5]$ 구간으로 설정, cdf 함숫값인 $y$를 먼저 구하기 다시 이것의 ppf 함숫값인 $z$를 구합니다.

이론적으로는 $x,z$가 동일해야 하므로, data_diff 를  써서 $x$와 $z$의 차이를 계산해 보았습니다. 결과를 볼까요?

 

항등함수가 나오네요

Maximum difference is 0.000000

Process finished with exit code 0

$x$, $z$사이에 차이도 없습니다.

이처럼 python은 norm.ppf 라는 함수를 제공하여 계산이 편리해지는 것이죠. 

 

---

그런데 만일! 이런 고급진 함수를 제공해 주지 않는다면 어떻게 할까요?  유용한 도구가 없으면 도구를 만들어 내면 됩니다.

도구가 없어 맨손으로 호두깨는 유투버

우리의 목적은 어떤 실수 $a\in (0,1)$ 에 대해

$$\Phi(x^\ast)= a\tag{1}$$

인 $x^\ast$를 찾아내는 것입니다. 물론 ppf 함수가 있으면 norm.ppf(a) 라는 식으로 구할 수 있죠. 

식(1)을 다시 써보면, 

함수 $f$를 $f(x) = \Phi(x)-a$라 할 때, $x^\ast$는 $f$의 해이다.

이렇게 쓸 수 있죠. 따라서 우리가 해 찾기 글에서 본 bisection method나 Newton-Raphson method를 쓸 수 있게 되는 것입니다. 우선 bisection method를 사용해 보죠

 

1. Bisection method로 해를 구하여 ppf 함수와 비교하기

Bisection method는 여기에서 다뤘습니다. 해당 글에서 썼던 code를 그대로 가져와 쓸 예정입니다.

 

python code를 보시죠.

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

def ppf_vs_bisection():
    x = np.arange(0.01, 0.99, 0.01)
    ppf_res=[]
    myfunc_res=[]
    data_diff=[]
    for i in range(len(x)):
        y=norm.ppf(x[i])
        z=find_root(x[i])
        ppf_res.append(y)
        myfunc_res.append(z)
        data_diff.append(y-z)
    print('Maximum difference is {:4,f}'.format(np.max(data_diff)))
    plt.plot(x, ppf_res, label='ppf graph')
    plt.plot(x, myfunc_res, label='besection graph')
    plt.legend()
    plt.show()

def find_root(u):
    f = lambda x: norm.cdf(x,0,1) - u
    xmin, xmax = -5, 5
    tol = 1e-8
    val = my_bisection(f,xmin,xmax,tol)
    return val

def my_bisection(f, a, b, tol):
    # approximates a root, R, of f bounded
    # by a and b to within tolerance
    # | f(m) | < tol with m the midpoint
    # between a and b Recursive implementation

    # check if a and b bound a root
    if np.sign(f(a)) == np.sign(f(b)):
        raise Exception(
            "The scalars a and b do not bound a root")

    # get midpoint
    m = (a + b) / 2

    if np.abs(f(m)) < tol:
        # stopping condition, report m as root
        return m
    elif np.sign(f(a)) == np.sign(f(m)):
        # case where m is an improvement on a.
        # Make recursive call with a = m
        return my_bisection(f, m, b, tol)
    elif np.sign(f(b)) == np.sign(f(m)):
        # case where m is an improvement on b.
        # Make recursive call with b = m
        return my_bisection(f, a, m, tol)

if __name__ == '__main__':
    ppf_vs_bisection()

간략히 설명하자면,

def find_root(u):
    f = lambda x: norm.cdf(x,0,1) - u
    xmin, xmax = -5, 5
    tol = 1e-8
    val = my_bisection(f,xmin,xmax,tol)
    return val

이 함수가 바로 $\Phi(x)=u$ 의 해를 bisection method로 구하는 part입니다.

 

그런데 사실 민감한 것이... $u=0$이거나 $u=1$일때 $x$를 구할 수 있을까요? 사실 $x$를 구할 수 없습니다.

$$\lim_{x\rightarrow \infty} \Phi(x) = 1 , \lim_{x\rightarrow -\infty}\Phi(x) =0$$

이지 딱 0과 1에 맞는 값을 찾을 수가 없는 상황이죠. 따라서 이 부분에서 bisection method가 에러가 납니다. 이것을 예외처리해서 제어할 수도 있지만.. 지금은 큰 틀에서 비교를 하는 중이니, $0$과 $1$은 분석에서 빼도록 합시다.

 

    x = np.arange(0.01, 0.99, 0.01)

이렇게요. 

code의 결과는 다음과 같습니다.

너무나 정확하여 아름답게 겹치는 두 차트

Maximum difference is 0.000000

Process finished with exit code 0

---

Newton-Raphson method를 안 볼 수 없겠죠? 해당 로직에 대해서는 여기에서 살펴보시면 됩니다.

 

2. Newton Raphson method로 해를 구하여 ppf 함수와 비교하기

 

python code는 아래와 같습니다.

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

def ppf_vs_NewtonRaphson():
    x = np.arange(0.01, 0.99, 0.01)
    ppf_res = []
    myfunc_res = []
    data_diff = []
    for i in range(len(x)):
        y = norm.ppf(x[i])
        z = find_root(x[i])
        ppf_res.append(y)
        myfunc_res.append(z)
        data_diff.append(y - z)
    print('Maximum difference is {:4,f}'.format(np.max(data_diff)))
    plt.plot(x, ppf_res, label='ppf graph')
    plt.plot(x, myfunc_res, label='NewtonRaphson graph')
    plt.legend()
    plt.show()

def find_root(u):
    f = lambda x: norm.cdf(x, 0, 1) - u
    df = lambda x: norm.pdf(x, 0, 1)
    x0 = 0
    tol = 1e-8
    val = my_newton(f, df, x0, tol)
    return val

def my_newton(f, df, x0, tol):
    # output is an estimation of the root of f
    # using the Newton Raphson method
    # recursive implementation
    if abs(f(x0)) < tol:
        return x0
    else:
        return my_newton(f, df, x0 - f(x0) / df(x0), tol)

if __name__ == '__main__':
    ppf_vs_NewtonRaphson()

my_newton 함수는 여기에서 그대로 가져왔습니다. 결과를 살펴볼까요?

아름다게 겹치는 두 그래프

Maximum difference is 0.000000

Process finished with exit code 0

 

이가 없으면 잇몸으로, 여러 아름다운 이론을 도입하여 해결할 수 있었습니다. 다음 글에서는 정규분포 난수 생성의 꽃과 같은 이론인 Box Muller Method를 다뤄보겠습니다.