1장 Introduction to Statistical Computing + Monte-Carlo Simulation

• Computing Basics

• Log-Exp Sum Trick

• Random Number Generators

• Monte-Carlo Integral

• Riemann-Stieltjes Integral

Computing Basics

아날로그 컴퓨터

아날로그 컴퓨터들의 경우 전기, 기계 등을 활용하여서 연속적으로 바뀌는 것을 활용합니다. 이것에 대한 대표적인 예시로는 아날로그 시계가 있습니다. 시계를 보면 $0 \rightarrow 1$ 으로 바로 쓱 올라가는게 아니라 연속적으로 올라갑니다. 하지만 이는 아날로그 컴퓨터들이며 요즘의 거의 모든 디지털 컴퓨터들은 이렇게 작동하지 않습니다.

디지털 컴퓨터

디지털 / 현대 컴퓨터들은 대부분 이진법 (Binary Language) 를 사용합니다 $(0 \text{ or }1)$. 컴퓨터의 CPU 에 따라서 이게 결정이 되는데, 우선은 이진법에 대해서 간단히 다루겠습니다.

이진법

$0$ 과 $1$ 으로만 숫자를 표현하는 방식이며 이것은 컴퓨터의 전기 신호가 켜졌는지 $(1)$ 아니면 꺼졌는지 $(0)$ 에 따라서 정보를 다루게 됩니다. 하나의 신호를 보통은 $\texttt{bit}$ 이라고 하며, 다음 예시를 보면 이해하기 편할겁니다.

표현할 숫자 (10-진법 기준) : $5 = 1 + 4 = 2^0 \cdot 1 + 2^1 \cdot 0 + 2^2 \cdot 1$

지금의 경우에는 $3 \texttt{ bit}$ 만 사용했기에 $2^3 -1= 8 - 1=7$ 까지 숫자를 표현하는게 가능합니다.

디지털 컴퓨터들의 경우 다음과 같은 숫자 표현들을 사용한다고 합니다. (128-bit 는 사용하기에는 비효율적)

그리고 Capacity 는 다음과 같습니다.

Capacity 를 넘게 된다면 Overflow / Underflow 문제가 발생하게 됩니다. 이를 방지하기 위한 방법들이 굉장히 많지만, 우선은 아래 내용을 진행하겠습니다.

그렇다면 소수들은 이 $\texttt{bit}$ 를 이용해서 어떻게 표현할까요? $\rightarrow$ Exponent & Mantissa

Mantissa 표현법

Fixed-Point Representation (고정소수점) 이 있으며 Floating-Point Representation (부동소수점) 이 있습니다.

• Fixed-Point Representation 의 경우 매 숫자마다 소수점을 미리 정해진 위치에 고정을 시켜줍니다.

• Floating-Point Representation 의 경우 매 숫자마다 소수점 위치를 지수에 따라서 변형시켜줍니다.

예시로 숫자 $75.31$ 이 있다고 하겠습니다.

• 이 숫자를 Fixed-Point Representation 으로 표현하면 똑같이 $75.31$ 이 됩니다.

• 반면에 Floating-Point Representation 으로 표현하면 $0.7531 \operatorname{E}02 = 0.7531 \cdot 10^2$ 으로 변형됩니다.

둘다 똑같은 값이며 그저 표현의 차이가 있을 뿐입니다. 그러면 부호 그리고 Mantissa 는 뭘까요?

• 부호의 경우에는 $+$ 아니면 $-$ 으로 둘 중에 하나이며, 이것도 Binary 으로 하나의 $\texttt{bit}$ 으로 표현해주면 됩니다.

• 그러면 64-bit 기준에서 1개의 bit 는 부호로 사용되기에 63 개로 이제는 표현해줘야 합니다. 이때 저희는 Floating-Point Representation 을 따르기에 먼저 숫자를 $0.7531 \operatorname{E} 02$ 으로 바꿔줍니다.

• $\operatorname{E}$ 뒤의 값들은 Exponent 으로 곱해줘야 하기에 이 값을 표현하기 위해서 비트를 $11$ 개를 또 할당시켜줍니다.
  • 이때 함정이 있다면 실제로 곱해지는 것은 $10^{\operatorname{e}}$ 가 아니라 $2^{\operatorname{e}}$ 입니다.
  • 추가적으로 $\operatorname{e}$ 의 값이 음수가 될 수도 있기에 원래는 $11$ 개로 $0 - 2047$ 까지 표현할 수 있지만, Bias Term 을 중앙값 $(1023)$ 으로 넣어줘서 $-1023 \sim 1023$ 까지 표현해줄 수 있게 됩니다.
  • Bias 를 안쓰고 $\operatorname{e}$ 의 첫번째 비트를 부등호로 바꿔서 위와 같이 똑같이 표현해주는 방식도 존재합니다.

• 그러면 이제는 $52$ 개의 비트가 남으며 이것으로 $\operatorname{E}$ 앞의 숫자를 표현해주면 됩니다.

정확하게는 다음과 같습니다.

소수는 어떻게 처리할까요?

사실 위의 경우에는 숫자를 어떻게 표현하냐에 대해서이며, 이는 이진법 10-진법 둘다 적용이 가능합니다. 하지만 컴퓨터는 이진법을 사용하기에 소수 부분만 갖고올 때도 10-진법에서 이진법으로 바꿔줘야 합니다. 그럼 이것을 어떻게 바꿀까요?

예시로 $0.7$ 와 같은 간단한 소수값을 가져오겠습니다. 이것을 다음과 같이 표현하는게 가능합니다.

$\begin{aligned} 0.7 &= (0.5) + 0.2 = (0.5 + 0 + 0.125) + 0.075 \\&= (0.5 + 0 + 0.125 + 0.0625) + 0.0125 {}_{10} \\&= (0.5 + 0 + 0.125 + 0.0625 + 0 + 0 + 0.0078125) + 0.0046875 {}_{10} = \cdots{}_{10} \\&= \left(1 \cdot 2^{-1} + 0 \cdot 2^{-2} + 1 \cdot 2^{-3} + 1 \cdot 2^{-4} + 0 \cdot 2^{-5} + 0 \cdot 2^{-6} + 1 \cdot 2^{-7} + \cdots \right) {}_{10} \\&= 0.1011001\cdots_2 \end{aligned}$

이것을 결국에는 이진법으로 처리하면 $0.1011001100110011001100110\cdots_2$ 으로 나오게 됩니다.

그러면 예시로 Mantissa 표현법에 대해서 조금 더 자세히 설명하겠습니다.

예시 1 $(63.745_{10})$ :

• 양수이기에 Sign 은 $0$

• 정수부분 :

  $63_{10} = 2^6 \cdot 1 + 2^5 \cdot 1 + 2^4 \cdot 1 + 2^3 \cdot 1 + 2^2 \cdot 1 + 2^1 \cdot 1 + 2^0 \cdot 1 \; {}_{10}= 1111111_{2}$

• 소수 부분 :

  $\begin{aligned} 0.745_{10} &= (0.5 + 0 + 0.125 + 0.0625 + 0.03125 + 0.015625 + \cdots) {}_{10} \\&= \left( 1 \cdot 2^{-1} + 0 \cdot 2^{-2} + 1 \cdot 2^{-2} + 1 \cdot 2^{-3} + 1 \cdot 2^{-4} + 1 \cdot 2^{-5} + 1 \cdot 2^{-6} + \cdots \right)_{10} \\&= 0.1011111\cdots_2 \end{aligned}$

• 합하면 :

  $\begin{aligned}63.745_{10} &= 111111.1011111010111000010\cdots_2 \\&= 1.111111011111010111000010\cdots_2 \times 2^5 \end{aligned}$

• 그러면 Mantissa 부분은 $111111011111010111000010\cdots$

• Exponent 는 $5$ 이며 Bias 는 $1023$ 이기에 $\operatorname{e} = 1023 + 5 = 1028$

• 최종 :

  

Errors in Computation

• Representation Error : 컴퓨터는 사실 유리수를 잘 표현하지 못합니다. 그렇기에 늘 근사를 합니다. $0.1$ 와 같은 단순한 숫자도 64-bit 기준으로는 $0.000110011\cdots_2$ 으로 표현이 됩니다. 그럼에도 불구하고 이 숫자를 표현하기 위해서는 끝자락에서 반올림을 시켜서 진행을 합니다. 즉 늘 실제 값이랑 컴퓨터가 표현한 값에는 Representation Error 가 존재한다는 것입니다 (딱 떨어지지 않는 이상).
  • 예시로 $0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 + 0.1 == 1$ 을 $\texttt{R}$ 코드에서 치면 $\texttt{FALSE}$ 가 나오게 됩니다.

• Truncation Error : 유리수도 표현 못하면 무리수도 당연하게 표현을 못합니다. 무리수들의 경우에는 그래서 Taylor Series 와 같은 것을 활용시키며, 이 값은 어느 시점 부터 Truncate 이 됩니다. 이때 Truncate 을 잘못하게 된다면 실제 값이랑 오차가 생기는 수밖에 없습니다.

Examples of Computation Error in $\texttt{R}$

Log-Exp Sum Trick

우선은 Overflow 와 Underflow 문제가 뭔지부터 파악해야됩니다. 앞서 저희는 Mantissa 표현이 어떻게 되는지에 대해서 알아보았으며, 이를 활용한다면 Overflow 와 Underflow 는 다음과 같이 정의가 가능합니다.

• Overflow : Mantissa 으로 표현할 수 있는 최대 숫자보다 값이 클 때 생기는 에러

• Underflow : Mantissa 으로 표현할 수 있는 최대 숫자보다 값이 작을 때 생기는 에러

생각보다 자주 일어나지는 않지만, 일어나더라도 해결할 수 있는 방안들이 많으며, 이 중에서 하나가 바로 Log-Exp Sum Trick 입니다. 로그와 exp 성질을 활용하면 다음과 같이 쓸 수 있습니다.

근데 이게 유용한 점이, 로그 안에서 Summation 값을 구하는데에도 유용해지게 됩니다.

이때 보통 $c$ 의 값은 $c := \max_{i\in [m]} x_i$ 으로 정해지게 됩니다.

예시로는 다음과 같이 활용할 수 있게 됩니다.

Random Number Generators

랜덤한 숫자를 생성하기 위해서는 어떻게 해야할까요?

물리적인 세상에서의 예시로는 주사위, 로또, 카드 썪기 등등의 방법들을 생각할 수가 있습니다. 하지만 모든게 계산으로 이뤄지는 전산에서는 결국에는 알고리즘을 활용해야됩니다.

이에 대한 예시를 Middle-Square Method, Linear Congruential Method 들이 존재합니다. 이때 RNG 가 가져야할 조건들이 다음과 같습니다.

Middle-Square Method

컴퓨터의 아버지이신 완전 말도 안되는 천재 John von Neumann 님이 만드신 너무나도 간단한 방법입니다. 그렇기에 사실 현대 사회에서는 아쉽게도 사용되기 애매합니다.

2 자리 RNG using Middle-Square Method

• 두자리 초깃값을 $x$ 라 하기.

• $x$ 를 제곱시켜서 4자리 수 만들기. 이때 앞자리가 빈다면 $0$ 으로 채워주기.

• $x$ 의 중간 숫자 2개를 다음 랜덤 숫자로 만들기.

x <- 어떤 2-Digit 초깃값
x <- x ** 2 # 제곱을 시켜주며 4-Digit 으로 바꿔줌, 4자리수가 안된다면 0 을 붙여준다.
x <- x[2:4] # 근데 사실 이건 R 코드에서 안되기에 % 를 활용하면서 이거를 해줘야함.
cat("새로운 숫자는 ", x)

4 자리로 늘렸을 때 나오는 Histogram 입니다.

뭐가 문제점일까요? 바로 Loop 입니다. 예시로 $3317$ 의 값을 초깃값으로 주었다고 해보겠습니다.

$3317 ^ 2 = 11002489 \rightarrow 0024^ 2 = 00000576 \rightarrow 0005 ^2 = 00000025 \rightarrow 0000^2 = 00000000$

이러한 문제점이 이 Middle-Square Method 에 존재하며, 초깃값에 따라서 어떨 때는 RNG 의 역할을 못할 때도, 어떨 때는 괜찮게 수행할 때도 있습니다.

Linear Congruential Method

1958 년도에 만들어진 RNG 이며 가장 영향력이 큰 RNG 기법 중에서 하나입니다. 방식은 다음과 같습니다.

1. $X_0, A, C, M$ 을 먼저 고른다.

2. $X_t \leftarrow (A X_{t-1} + C) \operatorname{mod} M$

3. 2 번을 계속해서 반복해준다.

이때 이 방법의 특징들이 있는데, 이는 다음과 같다.

• 이론적으로 최대 Period 는 $M$ 이다 (값들만 잘 정해준다면).

• $X_t$ 에 추가시키기 전에 값을 $M-1$ 으로 나눠주면 $0 \sim 1$ 사이의 값이 나오게 된다.

• 처음 값들을 잘 설정만 해준다면 균일분포처럼 나오게 된다.

그럼 예시 결과를 한번 보겠습니다.

보이시는 바와 같이 값들을 잘 선정하게 된다면 오른쪽과 같이 거의 균일분포의 모습을 띄우게 됩니다.

그러면 저희가 구현한 RNG 가 진정한 RNG 인 $\operatorname{Unif}(0,1)$ 과 같은 분포를 갖게 된다는 확인을 어떻게 할까요?

Statistical Tests for Uniformity

검정 기법들

• $\chi^2$-goodness-of-fit test

• Kolmogorov-Smirnov test

• Double-Root Residual

그래프 기법들

• 히스토그램 (Histogram)

• Probability Plot

• QQ-Plot

$\texttt{runif}(\cdot)$ 도 사실은 완전히 랜덤한 RNG 는 아니지만, 사용하기에 부족함이 없다.

Monte-Carlo Integral

적분의 경우 전산적으로 어떻게 계산할까요? 그것을 알기 위해서는 적분에 대해서 잠깐 얘기해보겠습니다.

위와 같은 적분 값이 있다고 하겠습니다. 그러면 이때 $f(x)$ 가 어떻게 생겼는지에 따라서 저희는 그저 손으로 적분도 할 수 있으며 너무 어렵다면 적분이 손으로는 현실적으로 불가능한 경우도 존재합니다. 그렇기에 저희는 이 적분 값을 근사 시켜주면 되며, Riemann Integral 를 활용하는 Numerical Integration 방법도 있지만 이 강의의 대부분은 Monte-Carlo Integration 방법을 사용한다고 보시면 됩니다.

위의 적분값을 다음과 같이 조금 다르게 해석해보겠습니다.

수식적으로 틀린 것은 지금까지 아무 것도 없습니다. 즉 다시말해 근사를 한 것은 아무것도 없습니다. 하지만 이때 여기서부터 저희는 Monte-Carlo Integration 을 적용시켜줍니다.

Law of Large Numbers

대수의 법칙으로도 알려져 있으며, 이것을 알기 위해서는 우선은 $X_1, X_2, \cdots, X_n \stackrel{\text{i.i.d}}{\sim} f(x;\theta)$ 이라고 정의하며, 특징으로는 $\mathbb{E}[X_i] = \mu, \quad \operatorname{Var}(X_i) = \sigma^2$ 이라고 하겠습니다.

추정량을 $\bar{X}_n =\frac{1}{n} \sum_{i=1}^n X_i$ 으로 정의하겠습니다. 그러면 다음과 같이 말할 수가 있습니다.

• Weak Law of Large Numbers :

  $${}^{\forall}\varepsilon > 0, \quad \lim_{n \rightarrow \infty} \mathbb{P}\left( \left| \bar{X}_n - \mu \right| > \varepsilon \right) = 0$$

  • 증명은 Chebyshev’s Inequality 사용하면 됩니다.

    Chebyshev’s Inequality 에 의하면,

    

    $$\operatorname{Var}(\bar{X}_n) = \frac{\sigma^2}{n}, \quad \text{By Chebyshev's Inequality, } k:= \frac{\sqrt{n} \varepsilon}{\sigma} \\~\\ {}^{\forall}\varepsilon > 0, \quad \lim_{n \rightarrow \infty} \mathbb{P}\left( \left| \bar{X}_n - \mu \right| \ge \varepsilon \right) < \frac{\sigma^2}{n \varepsilon^2} \stackrel{n\rightarrow \infty}{\rightarrow} 0$$

• Strong Law of Large Numbers :

  $$\mathbb P\left( \lim_{n \rightarrow \infty} \bar{X}_n = \mu \right) = 1$$

  즉 확률 1으로 $\bar{X}_n$ 은 $\mu$ 으로 수렴하게 된다는 말입니다.

이것을 활용하게 된다면, $\mathbb{E}[X] = \mathbb{E}[\bar{X}_n] = \mu \approx \bar{X}_n$ 으로 추정해도 된다는 뜻입니다.

Monte-Carlo Integration

아까 보았던 예시로 돌아가겠습니다.

이때 $U$ 는 그저 균일분포의 확률변수, 즉 RNG 와 같은 분포라고 보시면 됩니다. LLN 을 이용하면 다음과 같이 적을 수 있게 됩니다.

ㄹ이 방식을 이용해서 적분값을 추정하는 것이 바로 Monte-Carlo Integral 입니다. 당연하게도 지금의 경우에는 $x \in (0,1)$ 이라는 아주 간편한 구간 내에 있기 때문에 이렇게 쉽게 이해가 될 수 있지만, 이것이 $x \in (a,b)$, $x\in (0, \infty)$, $x\in (-\infty, \infty)$ 등의 구간에서도 활용이 충분히 가능합니다.

• $x \in (a,b)$ 인 케이스

  Change of Variables 를 적용시키면 됩니다. 예시로 위와 같은 Notation 으로 다음과 같은 상황을 고려해보겠습니다.

  $$y = \int_a^b f(x) dx$$

  이때 저희가 $u := \frac{x - a}{b - a}, \; z \in (0, 1)$ 으로 두며 Jacobian 은 $|J| = \left| \frac{dx}{du} \right| = |b-a|$ 으로 세팅한다면 다음과 같이 적분식이 변하게 됩니다. p.s. $x = a + (b-a) \cdot u$,

  $$y = \int_a^b f(x) dx = \int_0^1 f(a + (b-a) u) \cdot |J| du = \int_0^1 f(a + (b-a)u) (b-a)du \\~\\ = (b-a) \cdot \int_0^1 f(a+(b-a)u)du \approx (b-a) \cdot \frac{1}{n}\sum_{i=1}^n f(a + (b-a) U_i ), \\ U_i \stackrel{\text{i.i.d}}{\sim} \operatorname{Unif}(0, 1)$$

• $x \in (0,\infty)$ 인 케이스

  $(0, \infty) \rightarrow (0, 1)$ 으로 바꿔주는 함수가 무엇이 있을까? (이걸 만족하는 아무 함수 써도 되긴 함)

  

  • 우리는 $u = \frac{1}{x+1}$ 을 고려해볼 수가 있다. 준비를 해보자면,
    • $u := \frac{1}{x+1}, \; u \in (0, 1)$
    • $x = \frac{1}{u} - 1, \; dx = -\frac{1}{u^2} du$
    • $|J| = \left|\frac{dx}{du}\right| = \left| -\frac{1}{u^2} \right| = \frac{1}{u^2}$

  $$y = \int_0^\infty f(x) dx = \int_0^1 f\left(\frac{1}{u} - 1\right) |J| du = \int_0^1 \frac{1}{u^2} f\left(\frac{1}{u} - 1\right) du\\~\\ = \mathbb{E}\left[ \frac{1}{U_i^2} f\left(\frac{1}{U_i} - 1\right) \right] \stackrel{\text{LLN}}{\approx} \frac{1}{n}\sum_{i=1}^n {U_i^2} f\left(\frac{1}{U_i} - 1\right), \quad U_i\stackrel{\text{i.i.d}}{\sim} \operatorname{Unif}(0,1)$$

  • 혹은 $u = e^{-x}$ 도 고려해볼 수가 있다.

    

    • $u:= \exp(-x), u \in (0, 1)$
    • $x = -\log(u), \; dx = -\frac{du}{u}$
    • $|J| = \left| \frac{dx}{du} \right| = \left| -\frac{1}{u} \right| = \frac 1 u$

    $$y = \int_0^\infty f(x) dx = \int_0^1 f\left(-\log u \right) |J| du = \int_0^1 \frac{1}{u} f(-\log u)du \\~\\ = \mathbb{E}\left[ \frac{1}{u} f(-\log u) \right] \stackrel{\text{LLN}}{\approx} \frac{1}{n}\sum_{i=1}^n \frac{1}{U_i} f(-\log U_i), \quad U_i \stackrel{\text{i.i.d}}{\sim}\operatorname{Unif}(0, 1)$$

• $x \in (-\infty,\infty)$ 인 케이스

  $(-\infty, \infty) \rightarrow (0, 1)$ 으로 바꿔주는 함수가 무엇이 있을까? 대표적으로는 Sigmoid 함수가 있다 (Logit Transformation). 즉 $u := \frac{1}{1+e^{-x}}$ 를 고려해볼 수가 있다.

  

  • $u:= \frac{1}{1+e^{-x}}, \; u \in (0, 1)$
  • $x = \log\left( \frac{u}{1-u} \right), \; dx = \frac{1}{u(1-u)}du$
  • $|J| = \left| \frac{dx}{du} \right| = \left| \frac{1}{u (1-u)} \right| = \frac{1}{u(1-u)}$

  $$y = \int_{-\infty}^\infty f(x) dx = \int_0^1 f\left(\log\frac{u}{1-u}\right) |J| du \\~\\ = \int_0^1 \frac{1}{u(1-u)} f\left(\log\frac{u}{1-u}\right) du = \mathbb{E}\left[\frac{1}{U_i(1-U_i)} f\left(\log\frac{U_i}{1-U_i}\right) \right] \\~\\ \stackrel{\text{LLN}}{\approx} \frac{1}{n} \sum_{i=1}^n \frac{1}{U_i(1-U_i)} f\left(\log\frac{U_i}{1-U_i}\right), \quad U_i\stackrel{\text{i.i.d}}{\sim} \operatorname{Unif}(0, 1)$$

그렇다면 실제로 적용시킬만한 좋은 예시가 무엇이 있을까요? 바로 $\pi$ 값을 구하는 것입니다. 동그라미의 넓이는 초등학생도 알다시피 $A = \pi r^2$ 입니다. 그렇다면 이 동그라미를 $(-1,1)^2$ 에 그려보면 어떻게 나올까요?

보라색 구역의 경우 $x^2 + y^2 \le 1$ 의 구역을 그대로 그린 것입니다. 동그라미 넓이에 포함되는 구역을 구했으니, 이제는 균일분포를 활용시켜야합니다. 이때 저희가 알 수 있는 것은 원을 둘러싸는 정사각형이 $(-1, 1)^2$ 전체이며, 정사각형 내에서 아무 곳을 샘플하려면 그저 독립된 균일분포에서 한번 샘플하면 끝납니다.

그림에서의 원의 넓이는 그러면 $r = 1$ 이기에 $A_{\text{circle}} = \pi r^2 = \pi$ 이며, 정사각형의 넓이는 $A_{\text{square}} = (1 - (-1))^2 = 2^2 = 4$ 입니다. 그러면 정사각형 아무데나에서 샘플을 했을 때 좌표가 동그라미 안에 들어갈 비유는 다음과 같습니다.

그러면 수식적으로 적용을 시켜보자면,

• $(X,Y) \in \left(-1, 1\right)^2, \; X_i \stackrel{\text{i.i.d}}{\sim} \operatorname{Unif}(-1,1), \; Y_i \stackrel{\text{i.i.d}}{\sim} \operatorname{Unif}(-1,1)$

• $\mathbb{P}\left\{ (X,Y) \in \text{원} \right\} = \mathbb{P}\left\{ X^2 + Y^2 \le 1 \right\} = \frac{\pi}{4}$

저희가 알고 싶은 것은 사실 그러면 이 확률 값인 $\mathbb{P}\left\{ X^2 + Y^2 \le 1 \right\}$ 이며, 이를 추정하기 위해서 다음과 같이 세워주면 됩니다.

그러면 이것을 Monte-Carlo Integration 으로 바로 해보겠습니다.

그저 원과 함께 그린 결과 :

Trace Plot 결과 :

실제 값 결과 (또 다시 돌려서 위의 Trace Plot 과는 다를 수도 있음) :

• $\pi = 3.141593$

• $N = 1000, \; \hat \pi = 3.176, \; |\hat{\pi} - \pi| = 0.034407$

• $N = 10000, \; \hat \pi = 3.136, \; |\hat \pi - \pi | = 0.005592$

• $N = 100000, \; \hat \pi = 3.14004, \; |\hat \pi - \pi | = 0.001552$

• $N = 1000, \; \hat \pi = 3.142148, |\hat \pi - \pi | = 0.000555$

Riemann-Stieltjes Integral

적분이 있다면 이 적분을 구하는데 다양한 방법들이 존재하지만, 그 중 Riemann Integral 이 있습니다.

Riemann Integral

리만이란 수학자는 “함수 $f$ 가 정의된 구간 $[a,b]$ 에서 넓이를 어떻게 수학적으로 정의할까” 에 대한 질문을 답하기 위한 방법 중 가장 기본적인 적분법이 Riemann Integral 입니다.

아주 간단하게 $[a,b]$ 사이의 구간을 $n$ 개의 똑같은 길이의 작은 구간들로 쪼개갰습니다. 각 구간들의 값을 $x_i$ 라고 하며, 이는 그러면 다음과 같은 Partition 으로 표현할 수 있습니다.

작은 구간들은 그럼 총 $n$ 개로 $[x_i, x_{i-1}] ,\; i \in 1, \cdots, n$ 으로 표현할 수 있습니다.

이때 저희는 각각의 작은구간 $[x_i, x_{i-1}]$ 내에서 한 점 $\xi_i$ 를 골라줍니다. 그러면 Riemann Sum 을 다음과 같이 정의할 수 있습니다.

이때 가장 넓은 작은구간의 길이를 $\| P \| = \max_i \Delta x_i$ 으로 Mesh 이라고 정의할 수 있습니다.

만약에 이 작은구간들을 훨씬 더 촘촘히, 즉 $\| P \| \rightarrow 0$ 으로 두게 되어서 나온 Riemann Sum 이 $I$ 이라고 한다면, 다음과 같이 Riemann Integral 을 정의하는 것이 가능해집니다.

다시 말하자면 Partition $P$ 와 $\xi_i$ 를 골랐을 때 Riemann Sum 이 같은 값으로 수렴하면 Riemann 적분이 가능하다는 뜻입니다.

Riemann-Stieltjes Integral

아까 보았던 Riemann-Integral 은 사실 그저 기본 폭을 더해주는 것이었습니다. 즉 다시 말하자면 $x$-축의 폭에 대해서 나눈거나 다름 없습니다. 그렇다면 폭 대신에 어떤 함수 $h$ 의 증가분으로 바꿔주면 어떻게 될까요?

즉 다시 말하자면 $\Delta h_i = h(x_i) - h(x_{i-1})$ 을 사용하면 어떻게 될까요? 이게 사실은 Riemann-Integral 의 Generalized Version 입니다.

시작하자면, 이제는 분할을 $h(x)$ 를 기준으로 하기 때문에 저희가 구하고자 하는 적분은 다음과 같습니다.

작은구간들의 경우에는 우선은 Partition $P$ 는 다음과 같이 똑같이 정의해줍니다.

Riemann-Stieltjes Integral 의 정의로는 다음과 같습니다.

그러면 이에 해당하는 특수 케이스는 뭐가 있을까요?

• $g(x) = x$ 이면 $g(x_i) - g(x_{i-1}) = x_i - x_{i-1}$ 이기에 Riemann Integral 이랑 똑같아진다.

• $g(x)$ 가 계단 함수이면 $\text{점에서의 Jump }\times f \text{값}$ 들의 합으로 표현할 수 있습니다.

그렇다면 이를 기댓값에 적용시키면 어떻게 될까요?

Riemann-Stieltjes Integral on Expectation

Non-decreasing 하며 Right-Semicontinuous 한 함수 $f$ 에 대해서 저희는 다음과 같이 정의할 수 있게 됩니다.

여기까지는 그저 Riemann-Stieltjes Integral 과 똑같습니다. 하지만 기댓값을 이제는 적용시켜보겠습니다.

$X \sim F$ 이며 $F$ 가 CDF 이라고 할 때, 저희는 CDF 의 Non-decreasing 하면서 Right-Semicontinuous 하다는 특징을 활용하여 다음과 같이 임의의 함수 $f$ 에 대해서 $f(X)$ 에 대한 기댓값을 정의할 수 있습니다.

이때 사용한 가정은 바로 모든 $x$ 에 대해서 $F(x)$ 가 미분이 가능하다는 것이며, 미분이 불가능한 구간을 사실은 포함할 수도 있습니다. (물론 학부 수준에서는 요구가 안되긴 합니다)

그리고 이것은 $X$ 가 Continuous 한 경우에서 저렇게 전개되며 만약에 $X$ 가 Discrete 하다면 다음과 같이 전개가 됩니다.

Riemann-Stieltjes Integral 에 르베그 분해법을 활용한다면, 다음과 같이도 전개가 가능합니다.

• $F_{ac}$ : 절대적으로 연속인 구간에 대한 CDF

• $F_d$ : 이산적인 구간의 CDF

• $F_{sc}$ : 특이연속인 구간에 대한 CDF

이렇게 정의하게 된다면 엄밀히 따졌을 때 기댓값을 다음과 같이 정의할 수 있습니다.

어쨌든 본론으로 돌아오자면 저희가 CDF 을 Measure 으로 사용해서 얻는 가장 큰 이점은 바로 균일 분포에서의 샘플링을 통해서 CDF 를 구할 수가 있다는 것입니다.

예시로 $X_i \stackrel{\text{i.i.d}}{\sim} F, \; i=[n]$ 이라고 하겠습니다. 그러면 $F$ 를 추정하기 위해서 Empirical CDF $F_e$ 를 다음과 같이 정의하겠습니다.

이때 저희는 그러면 Monte-Carlo Integration 에서 평균을 사용하듯이 CDF 를 추정하면 되며, 이에 대한 이유는 다음과 같습니다.

이유 1 :

이유 2 :

Beyond Integral

사실 적분 뿐만 아니라 Monte Carlo Simulation 은 꽤나 다양하게 사용될 수가 있다.

• 기댓값을 활용한 다른 값들, 예시로 분산.

• 어떠한 분포에 대한 Sample Quantiles.

• 어떠한 확률변수의 분포밀도.