본 포스팅에서는 수식을 포함하고 있습니다.
티스토리 피드에서는 수식이 제대로 표시되지 않을 수 있으니
PC 웹 브라우저 또는 모바일 웹 브라우저에서 보시기 바랍니다.
이번 포스팅에서는 각종 기대값의 근사치를 구하거나 순서 통계량의 누적 분포를 구할 때 필요한 Gaussian Quadrature 적분 관련 논문인 Steen 외 2명의 논문 'Gaussian Quadratures for the Integrals'을 소개하고 여기에서 제시한 적분 근사 방법을 파이썬으로 구현해보고자 한다.
여기서 다루는 내용은 다음과 같다.
2. Computation of Weights and Abscissae
Abstract
Gaussian-Hermite Quadrature 적분은 아래와 같이 $(-\infty, \infty)$에서의 적분을 특정 지점에서의 함수값과 가중치의 합으로 근사하는 방법이다.
$$\int_{-\infty}^{\infty}\exp (-x^2)f(x)dx \approx \sum_{i=1}^nw_if(x_i)$$
하지만 구간 $(-\infty,\infty)$ 뿐 아니라 $(0,\infty)$ 또는 $(0, b), \;b<\infty$와 같은 구간에서의 적분도 많이 활용된다. 따라서 본 논문에서는 마지막 2개의 구간에 대해서 적분값을 근사하는 방법에 대해서 소개한다. 여기서 $w_i$ 는 가중치(weight), $x_i$는 가로 좌표점(abscissa)이라 하는데 가로 좌표점의 경우에는 특정 포인트 혹은 특정 지점으로 용어를 바꿔 사용하겠다.
본 논문에서는 $(0,\infty)$와 $(0, 1)$에 대해서 $2\leq n \leq 15$에서의 가중치와 특정 지점을 계산하였으며 각각의 값은 15자리의 유효숫자로 표시하였다.
1. Introduction
아래와 같은 적분 계산은 여러 분야에서 활용된다.
$$\int_0^b\exp (-x^2) f(x) dx \tag{1}$$
또는
$$\int_0^{\infty}\exp (-x^2) f(x) dx \tag{2}$$
$f(x)$가 다항함수와 같이 쉬운 함수라면 위의 적분은 Analytical 하게 구할 수 있으므로 문제가 되지 않지만 $f(x)$가 $x<0$에서 정의되지 않은 상황에서 복잡한 함수라면 적분을 직접 계산하는 것보다는 근사를 통해 계산하는 것도 고려해볼 만하다. 따라서 본 논문에서는 (1), (2)의 적분을 근사하는 방법을 소개한다.
기존의 Gauss-Laguerre quadrature 적분 근사는 $u=x^2$와 같은 변환이 필요한데 이 경우에는 원점에서 Singular 한 문제가 발생하여 구할 수 없다. 또 다른 잘 알려진 방법으로 Gauss-Hermite quadrature 적분 근사가 있는데 이 경우 $f$에 대한 even extension을 이용하여 구할 수도 있다. $g$를 $f$에 대한 even extension이라고 하자. $g$가 원점에서 높은 차수의 미분이 되지 않는다면 다시 말해 충분히 smooth하지 않다면 (2)에 대한 근사값이 $n$이 증가함에 따라 진동한다고 한다. 저자는 $f(x) = x^k, \;\;k=1, 3$에 대해서 Gauss-Hermite Quadrature 적분 근사 값과 직접 적분한 값과의 상대적 차이를 비교해보았다. 그 결과는 아래 Table 1에서 볼 수 있다. 실제로 높은 차수의 미분이 되지 않는 경우($k=1$) 상대적 오차가 굉장히 크고 $n$이 증가함에 따라 부호가 바뀌는 것을 볼 수 있다.
만약 $f$가 even이라면 기존의 Gauss-Hermite quadrature 적분 근사 방법이 본 논문에서 제시한 방법보다 더 좋을 수 있다고 한다. 왜냐하면 $\exp (-x^2)f(x)$가 원점에 대하여 대칭이 되므로 한쪽 방향의 특정 지점만 구하면 되기 때문이다. 즉, 기존에 필요한 특정 지점 개수보다 절반만 있으면 되기 때문에 계산이 빨라진다.
2. Computation of Weights and Abscissae
먼저 $\gamma_N$을 다음과 같이 정의한다.
$$\gamma_N = \int_0^b\exp (-x^2)p_N(x)^2dx$$
여기서 $p_N(x)$는 다음과 같이 재귀적으로 정의된다.
$$\begin{align}p_0(x) &= 1 \\ p_1(x) &= x-[1-\exp (-b^2)]/(\sqrt{\pi}\text{erf} (b)) \\ p_{k+1}(x) &= (x+\alpha_k)p_k(x) + \beta_kp_{k-1}(x)\;\; k=1, 2, \ldots \end{align}$$
이때 $\alpha_k$, $\beta_k$는 다음과 같이 정의된 상수이다.
$$\alpha_k = -\gamma_k^{-1}\int_0^b \exp (-x^2)xp_k(x)^2dx, \; \beta_k = -\gamma_k/\gamma_{k-1}$$
이제 가중치 $w_j$와 특정 지점 $x_j$을 구하자. $x_j, j=1, \ldots, n$은 $p_n(x)$의 $n$개 근이며 가중치는 다음과 같이 구한다.
$$w_j = \frac{\gamma_n}{(p_n'(x_j))p_{n-1}(x_j)}$$
이렇게 가중치와 특정 지점을 구했다면 (1)은 다음과 같이 근사할 수 있다.
$$\int_0^b\exp (-x^2)f(x)dx \approx \sum_{i=1}^nw_if(x_i)$$
여기서는 $(0, b)$ 위에서의 적분이지만 앞의 결과에서 $b\rightarrow \infty$로 하고 가중치와 다항식 그리고 특정 지점을 구해준다면 $(0, \infty)$에서의 적분 근사도 할 수 있다.
이 이후에는 Orthogonality를 이용하여 계산을 좀 더 쉽게 할 수 있는 방법이 나온다.
3. Implementation with Python
이제 파이썬을 이용하여 $n$이 주어졌을 때 가중치와 특정 지점을 계산하는 함수를 만들 것이다.
필요한 모듈을 먼저 임포트한다.
import numpy as np
from scipy.special import erf
다음으로 다음의 적분값을 계산하는 함수를 구한다.
$$\int_0^b \exp(-x^2)x^kdx, \;\; k=0, 1, \ldots $$
def get_exp_xpower_int(k, b=None):
'''
if b == None, b == infty
'''
if b:
if k==0:
return (np.sqrt(np.pi)/2)*erf(b)
elif k==1:
return 0.5*(1-np.exp(-np.square(b)))
else:
return ((k-1)*0.5)*get_exp_xpower_int(k-2,b)-0.5*np.power(b,k-1)*np.exp(-np.square(b))
else:
if k==0:
return np.sqrt(np.pi)/2
elif k==1:
return 1/2
else:
return ((k-1)*0.5)*get_exp_xpower_int(k-2)
다음으로 다항식 $p(x) = a_0 + a_1x + \cdots + a_nx^n$에 대하여 다음의 적분값을 계산하는 함수를 정의한다.
$$\int_0^b \exp(-x^2)p(x)dx$$
def get_exp_poly_int(coef,b=None):
res = 0
for i in range(len(coef)):
res += coef[i]*get_exp_xpower_int(i,b)
return res
다음으로 Orthogonal 다항식 $p_n(x)$의 계수를 계산해주는 함수를 정의한다.
def get_polynomial(n,b=None):
if n == 0:
return np.array([1])
elif n == 1:
if b:
return np.array([1, -(1-np.exp(-np.square(b)))/(np.sqrt(np.pi)*erf(b))])
else:
return np.array([1, -(1/np.sqrt(np.pi))])
else:
first = get_polynomial(n-1,b)
## calculate gamma
p = np.poly1d(first)
p2 = p*p ## square polynomial
p2_coef = np.flip(p2.coef) ## reverse coefficient
gamma_first = get_exp_poly_int(p2_coef,b) ## gamma[n-1]
p2_coef_shifted = np.insert(p2_coef,0,0) ## x*polynomial
integral_component = get_exp_poly_int(p2_coef_shifted,b)
alpha = -((1/gamma_first))*integral_component
first_plus_order = np.insert(first, len(first), 0) ## x*p(x)
first_plus_order = np.poly1d(first_plus_order)
first_multiplied_const = alpha*first ## alpha*p(x)
first_multiplied_const = np.poly1d(first_multiplied_const)
second = get_polynomial(n-2,b) ## second term
ps = np.poly1d(second)
ps2 = ps*ps
ps2_coef = np.flip(ps2.coef)
gamma_second = get_exp_poly_int(ps2_coef,b)
beta = -gamma_first/gamma_second
second_multiplied_const = beta*second
second_multiplied_const = np.poly1d(second_multiplied_const)
return (first_plus_order+first_multiplied_const+second_multiplied_const).coef
이제 마지막으로 가중치와 특정 지점을 계산하는 함수를 정의한다.
def get_weight_and_poins(n,b=None):
assert n > 0 and isinstance(n,int), 'n must be integer greater than 0'
coef = get_polynomial(n,b)
coef_prev = get_polynomial(n-1,b)
points = np.roots(coef)
p = np.poly1d(coef)
p_prev = np.poly1d(coef_prev)
p2_prev = p_prev*p_prev
p2_prev_coef = np.flip(p2_prev.coef)
gamma = get_exp_poly_int(p2_prev_coef,b)
first = p.deriv()(points)
second = p_prev(points)
denominator = first * second
weights = gamma/denominator
weights
return (weights, points)
이제 실험을 해보자. 먼저 $n=3$인 경우 $(0, \infty)$에서의 적분을 근사하기 위한 가중치와 특정 지점을 구해보자.
get_weight_and_poins(3)
결과를 논문의 결과와 비교해보자.
정확하게 똑같이 나왔다. 또한 $n$이 4~10인 경우에도 정확하게 일치하였다. 이번엔 $n=15$인 경우를 구해보자.
get_weight_and_poins(15)
비교해보면 가중치와 특정 지점의 값이 차이가 난다. 처음엔 내 코드가 잘못된 건 줄 알았는데 논문에서 $n=15$인 경우 $\gamma_{15}$의 계산값이 부정확할 수 있다고 한다. 따라서 이러한 부정확성이 가중치와 특정 지점을 계산하는데 영향을 미친 것으로 보인다. 또한 $b=1$로 지정하여 $(0, 1)$ 위에서의 적분값 근사를 위한 가중치와 특정 지점을 계산해보았더니 $n$이 3~8까지는 거의 정확하게 나왔고 9에서 약간의 차이를 보였으며 10에서 특정 지점이 음수값이 나왔다(무조건 양수만 나와야한다).
이 논문은 사실 Dixon's Q Test에서 기각값을 계산할 때 semi-infinite interval에서의 적분을 근사해야 할 필요가 생겨서 읽어보게 된 것이다. 근데 이게 은근 유용하게 쓰일 수 있고 또한 파이썬에서 이러한 적분값을 계산해주는 라이브러리가 없는 것 같아서 내가 나중에 쓰려고 함수를 구현해보았다. 찾아보니까 scipy 라이브러리에서 제공하는 scipy.integrate.quad 함수가 semi-infinite interval에서의 적분값을 계산해주더라.... 그래도 적분을 근사하는 여러 가지 방법에 대해서 알게 된 것 같아 뜻깊은 시간이었다.
댓글