이 글은
2022.10.25 - [수학의 재미/확률분포] - 이변량 정규분포(bivariate normal distribution)
이변량 정규분포(bivariate normal distribution)
어떤 이벤트가 두 개 이상의 확률변수에 의해서 결정될 때, 이벤트가 일어날 확률은 확률변수들 간의 결합 분포(joint distribution)에 의해 결정됩니다. 통계학적으로 널리 쓰이고 가장 중요한 정규
sine-qua-none.tistory.com
에서 이어집니다.
지난 글에서는 이변량 정규분포(bivariate normal distribution)의 개념과 pdf(probability density function), cdf(cumulative distribution function)의 정의를 알았습니다.
이제 이변량 정규분포의 pdf와 cdf 함숫값을 어떻게 구할 수 있는지 python code를 통해 알아보겠습니다.
Python Code (python 내장 함수)
우선 python의 scipy 라이브러리에서 제공해 주는 함수를 쓴 code입니다.
import numpy as np
from scipy.stats import multivariate_normal
from scipy.stats import norm
def PythonfuncVsNumerical_BND_CBND():
mean = np.zeros(2)
cov = np.identity(2)
rng = [0.5, 1]
rho = 0.3
cov[0, 1] = cov[1, 0] = rho
pyfunc = multivariate_normal(mean, cov)
print('Bivariate Normal PDF Result')
print('#1. Python function : {}'.format(pyfunc.pdf(rng)))
print('\n')
print('Bivariate Normal CDF Result')
print('#1. Python function : {}'.format(pyfunc.cdf(rng)))
if __name__ == '__main__':
PythonfuncVsNumerical_BND_CBND()
간단히 설명하자면,
from scipy.stats import multivariate_normal #다변량 정규분포 관련 여러 값들을 구하기위한 모듈
from scipy.stats import norm # 정규분포함수 관련 여러값들을 구하기 위한 모듈
# 이변량 정규분포의 값들을 구하기 위해 평균과 공분산을 설정
mean = np.zeros(2) # 2원소 배열의 평균
cov = np.identity(2) # 2*2사이즈의 모든 원소가 1인 행렬 지정
rng = [0.5, 1] # pdf, cdf 값을 구해보기 위한 2차원 좌료
rho = 0.3 # 두 변수의 상관계수 지정
cov[0, 1] = cov[1, 0] = rho # 상관계수 행렬 만듬○ $\mathbf{\mu}=(0,0), \Sigma =\textstyle{\begin{pmatrix} 1 &\rho \\ \rho& 1 \end{pmatrix}} $ 인 이변량 정규분포 설정합니다.
pyfunc = multivariate_normal(mean, cov)
print('Bivariate Normal PDF Result')
print('#1. Python function : {}'.format(pyfunc.pdf(rng)))
print('\n')
print('Bivariate Normal CDF Result')
print('#1. Python function : {}'.format(pyfunc.cdf(rng)))○ multivariate_normal(mean, cov)라는 클래스는 평균이 mean이고 공분산 행렬이 cov 인 다변량 정규분포를 설정한 클래스입니다.
○ multivariate_normal(mean, cov).pdf(point) 는설정된 다변량 정규분포의, 점 point에서의 확률밀도함숫값입니다.
○ multivariate_normal(mean, cov).cdf(point) 는설정된 다변량 정규분포의, 점 point에서의 누적분포함숫값입니다.
결과는 아래와 같이 나옵니다.
Bivariate Normal PDF Result
#1. Python function : 0.0989936325441632
Bivariate Normal CDF Result
#1. Python function : 0.6093086777999954
이 값은 당연히 참이겠지만, 다른 함수를 작성하여 이 값을 cross check 해보도록 하겠습니다.
Python Code : Numerical 버전
모든 확률변수의 cdf(누적분포함수)는 pdf(확률밀도함수)를 적분해서 나옵니다. 구체적으로, 확률변수 $X$의 pdf를 $f_X$라 하고, cdf를 $F_X$라 하면,
$$F_X(a) = \mathbb{P}(X\leq a) = \int_{-\infty}^{a} f_X(x)dx$$
입니다.
따라서 보통 pdf는 함수 형태로 주어지고 cdf는 적분 형태로 주어지는데, 적분을 코딩한다는 것은 구분구적법을 해야 하므로 시간도 많이 걸리고 또 값이 정확하지 않을 수도 있습니다. 이를 고민하던 학자들이 유명한 cdf들에 대해서는 적분을 하지 않고 수치해석적으로 cdf의 근삿값을 구하는 여러 결과를 도출해 냈는데요. 이러한 방법들은 이미 또 다른 여러 학자들에 의해 수치적인 오류가 없고 계산 소요시간도 훨씬 빠르다는 것이 보장되었습니다.
확률분포 중 가장 널리 쓰이고 매우 중요한 분포가 바로 정규분포이므로 이 정규분포의 cdf에 대한 수치해석적 접근이 없을 리 없겠죠. 여기 몇 가지 방법을 제시합니다.
다음의 방법들은 Haug의 "The Complete Guide to Option Pricing Formulas" 라는 책에서 발췌했습니다. 이 분의 프로그래밍 언어는 VBA였는데, 이를 Python에서 이용해 보기 위해 converting 해 보았습니다. 책에 대한 소개는 [금융공학] - 배리어옵션 -UOC 옵션 #1 수학공식(Closed form)에서도 한 바 있으니 한번 살펴보시기 바랍니다.
중요한 점은 밑에 소개할 3개의 cdf 함수들은 모두 이변량 표준정규분포(각 평균 0, 표준편차 1)에 대한 근사식입니다.
(python에서 제공하는 코드는 임의의 평균, 공분산 행렬에 대해 연산 가능)
Drezner-Wecolowsky Algorithm (1990)
def CBND2(a, b, rho):
# The cumulative bivariate normal distribution function
# Drezner-Wesolowsky 1990 simple algorithm
x = np.array([0.018854042, 0.038088059, 0.0452707394, 0.038088059, 0.018854042])
y = np.array([0.04691008, 0.23076534, 0.5, 0.76923466, 0.95308992])
sum = 0
for i in range(5):
P = y[i] * rho
g = 1 - P ** 2
sum += x[i] * np.exp((2 * a * b * P - a ** 2 - b ** 2) / g / 2) / np.sqrt(g)
res = rho * sum + norm.cdf(a) * norm.cdf(b)
return res
Drezner (1978)
def CBND3(a, b, rho):
# The cumulative bivariate normal distribution function
# Based on Drezner-1978
x = np.array([0.24840615, 0.39233107, 0.21141819, 0.03324666, 0.00082485334])
y = np.array([0.10024215, 0.48281397, 1.0609498, 1.7797294, 2.6697604])
a1 = a / np.sqrt(2 * (1 - rho ** 2))
b1 = b / np.sqrt(2 * (1 - rho ** 2))
if a <= 0 and b <= 0 and rho <= 0:
sum = 0
for i in range(5):
for j in range(5):
sum += x[i] * x[j] * np.exp(
a1 * (2 * y[i] - a1) + b1 * (2 * y[j] - b1) + 2 * rho * (y[i] - a1) * (y[j] - b1))
cbnd3_result = np.sqrt(1 - rho ** 2) / np.pi * sum
elif a <= 0 and b >= 0 and rho >= 0:
cbnd3_result = norm.cdf(a) - CBND3(a, -b, -rho)
elif a >= 0 and b <= 0 and rho >= 0:
cbnd3_result = norm.cdf(b) - CBND3(-a, b, -rho)
elif a >= 0 and b >= 0 and rho <= 0:
cbnd3_result = norm.cdf(a) + norm.cdf(b) - 1 + CBND3(-a, -b, rho)
elif a * b * rho > 0:
rho1 = (rho * a - b) * np.sign(a) / np.sqrt(a ** 2 - 2 * rho * a * b + b ** 2)
rho2 = (rho * b - a) * np.sign(b) / np.sqrt(a ** 2 - 2 * rho * a * b + b ** 2)
delta = (1 - np.sign(a) * np.sign(b)) / 4
cbnd3_result = CBND3(a, 0, rho1) + CBND3(b, 0, rho2) - delta
return cbnd3_result
Drezner & Wecolowsky Algorithm (Another version)
def CBND4(A, b, rho):
# 'modified/corrected from the second function in Drez & Wes paper pg. 105
# '0/0 case resolved by l'H rule
x = [0.04691008, 0.23076534, 0.5, 0.76923466, 0.95308992]
W = [0.018854042, 0.038088059, 0.0452707394, 0.038088059, 0.018854042]
h1 = A
h2 = b
h12 = (h1 * h1 + h2 * h2) / 2
if np.abs(rho) >= 0.7:
r2 = 1 - rho * rho
r3 = np.sqrt(r2)
if rho < 0:
h2 = -h2
h3 = h1 * h2
h7 = np.exp(-h3 / 2)
if np.abs(rho) < 1:
h6 = np.abs(h1 - h2)
h5 = h6 * h6 / 2
h6 = h6 / r3
AA = 0.5 - h3 / 8
ab = 3 - 2 * AA * h5
LH = 0.13298076 * h6 * ab * (1 - norm.cdf(h6)) - np.exp(-h5 / r2) * (ab + AA * r2) * 0.053051647
for i in range(5):
r1 = r3 * x[i]
rr = r1 * r1
r2 = np.sqrt(1 - rr)
if h7 == 0:
h8 = 0
else:
h8 = np.exp(-h3 / (1 + r2)) / r2 / h7
LH = LH - W[i] * np.exp(-h5 / rr) * (h8 - 1 - AA * rr)
res = LH * r3 * h7 + norm.cdf(min(h1, h2))
if rho < 0:
res = norm.cdf(h1) - res
else:
h3 = h1 * h2
LH = 0
if rho != 0:
for i in range(5):
r1 = rho * x[i]
r2 = 1 - r1 * r1
LH = LH + W[i] * np.exp((r1 * h3 - h12) / r2) / np.sqrt(r2)
res = norm.cdf(h1) * norm.cdf(h2) + rho * LH
return res
pdf, cdf 함숫값 크로스 체크
누적 분포 함수에 대해서는 여러 근삿값 구하는 식을 소개했으니, 이것들과 파이썬에서 제공하는 함숫값을 비교해 봐야겠습니다. 아울러
def BND(x, y, rho):
res = 1 / (2 * np.pi * np.sqrt(1 - rho ** 2)) * np.exp(
-1 / (2 * (1 - rho ** 2)) * (x ** 2 + y ** 2 - 2 * x * y * rho))
return res처럼, pdf 함수를 직접 구해서, 이것을 가지고 python 내장 함수와 비교도 해보도록 하겠습니다. 위 코드는
[수학의 재미/확률분포] - 이변량 정규분포(bivariate normal distribution)를 참고하시면 됩니다.
자, 크로스체킹 해보겠습니다. 제일 먼저 기술한 PythonfuncVsNumerical_BND_CBND() 코드만 아래처럼 바꿉니다.
def PythonfuncVsNumerical_BND_CBND():
mean = np.zeros(2)
cov = np.identity(2)
rng = [0.5, 1]
rho = 0.3
cov[0, 1] = cov[1, 0] = rho
pyfunc = multivariate_normal(mean, cov)
num_res3 = CBND3(rng[0], rng[1], rho)
num_res2 = CBND2(rng[0], rng[1], rho)
num_res4 = CBND4(rng[0], rng[1], rho)
print('Bivariate Normal PDF Result')
print('#1. Python function : {}'.format(pyfunc.pdf(rng)))
print('#2. Direct calculation: {}'.format(BND(rng[0],rng[1],rho)))
print('\n')
print('Bivariate Normal CDF Result')
print('#1. Python function : {}'.format(pyfunc.cdf(rng)))
print('#2. Numerical function2 : {}'.format(num_res2))
print('#3. Numerical function3 : {}'.format(num_res3))
print('#4. Numerical function4 : {}'.format(num_res4))
결과를 볼까요?
Bivariate Normal PDF Result
#1. Python function : 0.0989936325441632
#2. Direct calculation: 0.0989936325441632
Bivariate Normal CDF Result
#1. Python function : 0.6093086777999954
#2. Numerical function2 : 0.6093086775071506
#3. Numerical function3 : 0.6093084461096385
#4. Numerical function4 : 0.6093086775071506어떻습니까?
○ pdf 값은 완전 똑같은 걸로 봐서, 실제 위의 기술한 BND라는 함수를 쓴 거 같네요.
○ cdf 값도 무시무시할 정도록 비슷합니다. 하지만 뭔가 Drezner가 주도했던 방법과는 다른 알고리즘으로 파이썬 내장 함수가 설계된 것 같네요.
모든 수치해석적 방법이 사칙 연산 위주로 이루어져 구분구적법을 통한 적분보다는 비교할 수 없이 빠르고, 값 또한 정확하니, 위 세 가지 중 한 알고리즘을 골라 사용하면 무리가 없을 것 같습니다.
'수학의 재미 > 확률분포' 카테고리의 다른 글
| 이변량 정규분포(bivariate normal distribution) (0) | 2022.10.25 |
|---|---|
| 상관관계가 있는 두개의 표준정규분포 난수 구하기 (0) | 2022.09.16 |
| 확률측도를 바꿉시다: Girsanov Theorem (0) | 2022.06.28 |
| 술먹고 걷기(Random Walk) #2: 동서남북 내키는대로~ (0) | 2022.06.25 |
| 술먹고 걷기(Random Walk) #1: 길따라 걷기 (0) | 2022.06.25 |
댓글