몬테카를로 적분

섬네일|몬테카를로 적분 그림. 이 예시에서 영역 D는 안쪽 원이고 영역 E는 정사각형이다. 정사각형의 면적(4)은 쉽게 계산할 수 있으므로, 원의 면적(π*1.0<sup>2</sup>)은 원 안의 점(40)과 총 점 수(50)의 비율(0.8)로 추정할 수 있으며, 원의 면적은 4*0.8 = 3.2 ≈ π로 근사된다.

수학에서 **몬테카를로 적분**(Monte Carlo integration)은 난수를 사용하여 수치 적분을 수행하는 기술이다. 이는 정적분을 수치적으로 계산하는 특정 몬테카를로 방법이다. 다른 알고리즘은 보통 규칙적인 격자에서 피적분 함수를 평가하지만, 몬테카를로는 피적분 함수를 평가할 지점을 무작위로 선택한다. 이 방법은 특히 고차원 적분에 유용하다.

몬테카를로 적분을 수행하는 다양한 방법이 있으며, 예를 들어 균등 표집, 층화표집, 중요도 샘플링, 순차적 몬테카를로(파티클 필터라고도 함), 그리고 평균장 입자 방법 등이 있다.

## 개요
수치 적분에서 사다리꼴 공식과 같은 방법은 결정론적 알고리즘을 사용한다. 반면 몬테카를로 적분은 비결정론적 접근 방식을 사용한다. 각 실현은 다른 결과를 제공한다. 몬테카를로에서 최종 결과는 오차 막대를 포함한 정확한 값의 근사치이며, 정확한 값은 이 오차 막대 안에 있을 가능성이 높다.

몬테카를로 적분이 다루는 문제는 중적분의 계산이다.
I = \int_{\Omega}f(\overline{\mathbf{x}}) \, d\overline{\mathbf{x}}
여기서 Ω는 \mathbb{R}^m의 부분집합이며, 부피는 다음과 같다.
V = \int_{\Omega}d\overline{\mathbf{x}}

단순 몬테카를로 접근 방식은 Ω에서 점들을 균일하게 표집하는 것이다. N개의 균일 표본이 주어지면,
\overline{\mathbf{x}}_1, \cdots, \overline{\mathbf{x}}_N\in \Omega,

I는 다음과 같이 근사될 수 있다.
 I \approx Q_N \equiv V \frac{1}{N} \sum_{i=1}^N f(\overline{\mathbf{x}}_i) = V \langle f\rangle.

이는 큰 수의 법칙에 의해 다음이 보장되기 때문이다.
 \lim_{N \to \infty} Q_N = I.

QN으로부터 I의 추정치가 주어지면, QN의 오차 막대는 분산의 비편향 추정량을 사용하여 표본 분산으로 추정할 수 있다.

 \mathrm{Var}(f)=\mathrm{E}(\sigma_N^2) \equiv \frac{1}{N-1} \sum_{i=1}^N \mathrm{E}\left [\left (f(\overline{\mathbf{x}}_i) - \langle f \rangle \right )^2\right ]. 
이는 다음으로 이어진다.
 \mathrm{Var}(Q_N) = \frac{V^2}{N^2} \sum_{i=1}^N \mathrm{Var}(f) = V^2\frac{\mathrm{Var}(f)}{N} = V^2\frac{\mathrm{E}(\sigma_N^2)}{N}.

수열
 \left \{ \mathrm{E}(\sigma_1^2), \mathrm{E}(\sigma_2^2), \mathrm{E}(\sigma_3^2), \ldots \right \} 
이 Var(f)와 동일하여 유한하다고 가정되는 한, 이 분산은 1/N에 따라 점근적으로 0으로 감소한다. 따라서 QN의 오차 추정치는 다음과 같다.
\delta Q_N\approx\sqrt{\mathrm{Var}(Q_N)}=V\frac{\sqrt{\mathrm{Var}(f)}}{\sqrt{N}},
이는 \tfrac{1}{\sqrt{N}}로 감소한다. 이것은 V를 곱한 표준 오차이다.
이 결과는 적분의 차원 수에 의존하지 않는데, 이는 몬테카를로 적분이 차원에 지수적으로 의존하는 대부분의 결정론적 방법에 비해 약속된 장점이다.
결정론적 방법과는 달리 오차 추정치가 엄격한 오차 한계가 아니라는 점에 주목하는 것이 중요하다. 무작위 표집은 오차를 과소평가할 수 있는 피적분 함수의 중요한 특징을 모두 밝혀내지 못할 수도 있기 때문이다.

단순 몬테카를로는 간단한 예시에는 작동하지만, 결정론적 알고리즘보다 개선될 수 있는 것은 문제별 표집 분포를 사용하는 알고리즘에서만 가능하다.
적절한 표본 분포를 사용하면 거의 모든 고차원 피적분 함수가 매우 국지적이며 작은 부분 공간만이 적분에 현저하게 기여한다는 사실을 활용할 수 있다.
몬테카를로 문헌의 상당 부분은 오차 추정치를 개선하는 전략을 개발하는 데 할애된다. 특히, 층화표집(영역을 부분 영역으로 나누는 것)과 중요도 샘플링(비균일 분포에서 샘플링하는 것)이 이러한 기술의 두 가지 예시이다.

### 예시
몬테카를로 적분의 전형적인 예시는 π의 추정이다. 다음 함수를 고려해 보자.
H\left(x,y\right)=\begin{cases}
1 & \text{if }x^{2}+y^{2}\leq1\\
0 & \text{else}
\end{cases}
그리고 Ω = [−1,1] × [−1,1]이고 V = 4인 집합을 고려한다. 다음을 주목하자.
I_\pi = \int_\Omega H(x,y) dx dy = \pi.

따라서 몬테카를로 적분으로 π 값을 계산하는 대략적인 방법은 Ω에서 N개의 난수를 선택하고 다음을 계산하는 것이다.
Q_N = 4 \frac{1}{N}\sum_{i=1}^N H(x_{i},y_{i})

오른쪽 그림에서 상대 오차 \tfrac{Q_N-\pi}{\pi}는 N의 함수로 측정되며, \tfrac{1}{\sqrt{N}}를 확인시켜준다.

### C/C++ 예시

```c
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
int main() {
    // Initialize the number of counts to 0, setting the total number to be 100000 in the loop.
    int throws = 99999, insideCircle = 0;
    double randX, randY, pi;

    srand(time(NULL));

    // Checks for each random pair of x and y if they are inside circle of radius 1.
    for (int i = 0; i < throws; i++) {
        randX = rand() / (double) RAND_MAX;
        randY = rand() / (double) RAND_MAX;
        if (randX * randX + randY * randY < 1) {
            insideCircle++;
        }
    }

    // Calculating pi and printing.
    pi = 4.0 * insideCircle / throws;
    printf("%lf\n", pi);
}
```


### 파이썬 예시
파이썬으로 작성되었다.

```numpy
import numpy as np

rng = np.random.default_rng(0)

throws = 2000
radius = 1

# Choose random X and Y data centered around 0,0
x = rng.uniform(-radius, radius, throws)
y = rng.uniform(-radius, radius, throws)

# Count the times (x, y) is inside the circle,
# which happens when sqrt(x^2 + y^2) <= radius.
inside_circle = np.count_nonzero(np.hypot(x, y) <= radius)

# Calculate area and print; should be closer to Pi with increasing number of throws
area = (2 * radius)**2 * inside_circle / throws
print(area)
```


### Wolfram Mathematica 예시
아래 코드는 매스매티카에서 몬테카를로 방법을 사용하여
f(x) = \frac{1}{1+\sinh(2x)\log(x)^2}
함수를 0.8 범위에서 적분하는 과정을 설명한다.


```mathematica
func[x_] := 1/(1 + Sinh[2*x]*(Log[x])^2);

(*Sample from truncated normal distribution to speed up convergence*)
Distrib[x_, average_, var_] :=   PDF[NormalDistribution[average, var], 1.1*x - 0.1];
n = 10;
RV = RandomVariate[TruncatedDistribution[{0.8, 3}, NormalDistribution[1, 0.399]], n];
Int = 1/n Total[func[RV]/Distrib[RV, 1, 0.399]]*Integrate[Distrib[x, 1, 0.399], {x, 0.8, 3}]

NIntegrate[func[x], {x, 0.8, 3}] (*Compare with real answer*)
```


## 재귀적 층화표집

**재귀적 층화표집**은 1차원 적응형 구적법을 다차원 적분으로 일반화한 것이다. 각 재귀 단계에서 일반 몬테카를로 알고리즘을 사용하여 적분과 오차를 추정한다. 오차 추정치가 필요한 정확도보다 크면 적분 부피를 부분 부피로 나누고 그 부분 부피에 절차를 재귀적으로 적용한다.

일반적인 '절반으로 나누기' 전략은 다차원에서는 너무 빨리 부분 부피의 수가 증가하여 추적하기 어렵기 때문에 작동하지 않는다. 대신 어떤 차원을 따라 분할해야 가장 많은 이득을 가져올지 추정하고 이 차원을 따라 부피를 분할한다.

층화표집 알고리즘은 함수의 분산이 가장 큰 영역에 샘플링 지점을 집중시켜 전체 분산을 줄이고 샘플링을 더 효과적으로 만든다. 이는 그림에서 볼 수 있다.

널리 사용되는 MISER 루틴은 유사한 알고리즘을 구현한다.

### MISER 몬테카를로
MISER 알고리즘은 재귀적 층화표집에 기반한다. 이 기술은 가장 높은 분산을 가진 영역에 적분 지점을 집중시켜 전체 적분 오차를 줄이는 것을 목표로 한다.

층화표집의 아이디어는 서로소 집합인 두 영역 a와 b에 대해 몬테카를로 적분 추정치 E_a(f)와 E_b(f), 그리고 분산 \sigma_a^2(f)와 \sigma_b^2(f)가 있을 때, 결합된 추정치의 분산 Var(f)가 다음과 같다는 관찰에서 시작한다.
E(f) = \tfrac{1}{2} \left (E_a(f) + E_b(f) \right )
다음과 같이 주어진다.
\mathrm{Var}(f) = \frac{\sigma_a^2(f)}{4 N_a} + \frac{\sigma_b^2(f)}{4 N_b}

이 분산은 점들을 다음과 같이 분포시킬 때 최소화됨을 보여줄 수 있다.
\frac{N_a}{N_a + N_b} = \frac{\sigma_a}{\sigma_a + \sigma_b}

따라서 가장 작은 오차 추정치는 각 하위 영역에서 함수의 표준 편차에 비례하여 표본 지점을 할당함으로써 얻어진다.

MISER 알고리즘은 각 단계에서 하나의 좌표축을 따라 적분 영역을 이등분하여 두 개의 하위 영역을 만든다. 방향은 가능한 모든 d개의 이등분을 검토하고 두 하위 영역의 결합된 분산을 최소화할 것을 선택하여 결정된다. 하위 영역의 분산은 현재 단계에서 사용 가능한 총 점 수의 일부를 샘플링하여 추정된다. 동일한 절차는 최적의 이등분에서 나온 두 반공간 각각에 대해 재귀적으로 반복된다. 나머지 표본 점들은 Na와 Nb 공식에 따라 하위 영역에 할당된다. 적분 점의 이러한 재귀적 할당은 사용자가 지정한 깊이까지 계속되며, 각 하위 영역은 일반 몬테카를로 추정치를 사용하여 적분된다. 이 개별 값과 그 오차 추정치는 위로 결합되어 전체 결과와 그 오차 추정치를 제공한다.

## 중요도 샘플링

다양한 중요도 샘플링 알고리즘이 있다.
### 중요도 샘플링 알고리즘

중요도 샘플링은 몬테카를로 적분을 수행하는 매우 중요한 도구를 제공한다. 이 방법에서 중요도 샘플링의 주요 결과는 \overline{\mathbf{x}}의 균일 샘플링이 더 일반적인 선택의 특정 경우이며, 이 경우 샘플은 임의의 분포 p(\overline{\mathbf{x}})에서 추출된다는 것이다. 아이디어는 p(\overline{\mathbf{x}})를 QN 측정의 분산을 감소시키도록 선택할 수 있다는 것이다.

중심이 0이고 σ = 1인 가우스 함수를 −1000에서 1000까지 수치적으로 적분하려는 다음 예시를 고려해 보자. 당연히 [−1000, 1000] 구간에서 샘플이 균일하게 추출되면, 그 중 극히 일부만이 적분에 유의미할 것이다. 이는 샘플이 선택되는 다른 분포를 선택함으로써 개선될 수 있는데, 예를 들어 중심이 0이고 σ = 1인 가우스 분포에 따라 샘플링하는 것이다. 물론 "올바른" 선택은 피적분 함수에 크게 의존한다.

정식으로, 분포에서 선택된 샘플 집합이 주어지면
p(\overline{\mathbf{x}}) : \qquad \overline{\mathbf{x}}_1, \cdots, \overline{\mathbf{x}}_N \in V, 
I에 대한 추정량은 다음과 같다.
 Q_N \equiv \frac{1}{N} \sum_{i=1}^N \frac{f(\overline{\mathbf{x}}_i)}{p(\overline{\mathbf{x}}_i)}

직관적으로, 이는 특정 샘플을 다른 샘플보다 두 배 많이 선택하면, 다른 샘플보다 절반의 가중치를 부여한다는 것을 의미한다. 이 추정량은 p(\overline{\mathbf{x}})가 상수인 경우인 균일 샘플링에 자연스럽게 유효하다.

메트로폴리스-헤이스팅스 알고리즘은 p(\overline{\mathbf{x}})에서 \overline{\mathbf{x}}를 생성하는 데 가장 많이 사용되는 알고리즘 중 하나이며, 따라서 적분을 효율적으로 계산하는 방법을 제공한다.

### VEGAS 몬테카를로

VEGAS 알고리즘은 적분 영역을 여러 번 통과하여 함수 f의 히스토그램을 생성함으로써 정확한 분포를 근사한다. 각 히스토그램은 다음 통과를 위한 표본 분포를 정의하는 데 사용된다. 점근적으로 이 절차는 원하는 분포로 수렴한다. 히스토그램 빈의 수가 Kd와 같이 증가하는 것을 피하기 위해, 확률 분포는 분리 가능한 함수로 근사된다.
g(x_1, x_2, \ldots) = g_1(x_1) g_2(x_2) \ldots 
따라서 필요한 빈의 수는 Kd개에 불과하다. 이는 피적분 함수를 좌표축에 투영하여 함수의 극대점을 찾는 것과 동일하다. VEGAS의 효율성은 이 가정의 유효성에 달려 있다. 피적분 함수의 극대점이 잘 국지화되어 있을 때 가장 효율적이다. 피적분 함수를 대략적으로 분리 가능한 형태로 다시 작성할 수 있다면 VEGAS를 사용한 적분 효율성이 증가할 것이다. VEGAS는 여러 추가 기능을 포함하며 층화표집과 중요도 샘플링을 모두 결합한다.

## 같이 보기
* 준몬테카를로 방법
* 보조장 몬테카를로
* 통계물리학의 몬테카를로 방법
* 몬테카를로 방법
* 분산 줄이기

## 내용주

## 각주

* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 

## 외부 링크
* 카페 수학: 몬테카를로 적분 : 몬테카를로 적분을 설명하는 블로그 기사 (원리, 가설, 신뢰 구간)
* Boost.Math : 단순 몬테카를로 적분: C++ 단순 몬테카를로 루틴 문서
* 통계 물리학 문제에 적용된 몬테카를로 애플릿

분류:몬테카를로 방법
분류:코드 예시에 관한 문서