---
title: 數學軟體實作 - Monte Carlo integration
tags: 2020 Fall - 數學軟體實作
GA: G-77TT93X4N1
---

# Lecture - Monte Carlo integration

---

## One-dimensional integration

### Mean value theorem for definite integrals
Let $f : [a, b] \to R$ be a continuous function. Then there exists $c\in [a, b]$ such that
$$
\int^b_a f(x)\,dx = f(c)(b-a),
$$
and
$$
f(c) = \frac{1}{(b-a)}\,\int^b_a f(x)\,dx.
$$

---

$f(c)$ 事實上就是函數 $f$ 在 $[a. b]$ 區間的平均值.

因此, 如果我們知道一個函數在某個區間的平均值, 稱之為 $f^*$, 那將這平均值乘以區間大小, 我們就可求得這函數在區間內的積分.

要如何估算平均值?
* 取 $N$ 個 sample 求平均

給定一個函數 $f(x)$ 與區間 $[a, b]$. 我們在這區間內取點 $x_i\in[a, b]$. 接著我們估算 $f^*$:
$$
f^* \approx \frac{1}{N}\sum^N_{i=1} f(x_i)
$$
知道 $f^*$ 我們就可以算函數在區間內積分了, $\int^b_a f(x)\,dx = f^*(b-a).$

#### Example: Approximate $\int^1_0 \sin(x)\,dx$

%%
% Use Monte-carlo method to approximate 
% \int^b_a f(x)\,dx
% f(x) 定義在最後面
%

%%
N = 1000;                           % 取 N 個點
a = 0;                              % [a, b] 區間
b = 1;

sum_all = 0;
for ii=1:N                          % 以迴圈隨機取點做函數值加總
    x = rand(1)*(b-a)+a;            % rand 取值在在 [0,1] 區間
    sum_all = sum_all + f(x);       % 做個橫移使得 x 在 [a, b] 區間內
end
mean_val = sum_all/N;               % 估計 f^*
integral_val = mean_val*(b-a);      % 估計積分
disp(['積分值為 ', num2str(integral_val)])

function y = f(x);
    y = sin(x);
end

### Exercise 1

令上面這程式所算出的值為 $a_N$, 這樣我們可以定義出一個數列 $\{a_n\}$. 試寫一程式求此數列極限值.

#### Caution
事實上, 這並不是一個數列, 這是所謂的離散型隨機變數 (discrete random variable). 所以它的"收斂"應該要看標準差而不是以上的兩項差.

##### Key words: 平均值(mean)，標準差(standard deviation, std)，中央極限定理(central limit theorem, CLT)

---

## [Monte Carlo integration](https://en.wikipedia.org/wiki/Monte_Carlo_integration)

> 假設我們在 $[-1,1]\times[-1,1]$ 這方塊裡隨機取 $N$ 個點, 會有多少個點落在以原點為圓心半徑為 $1$ 的圓內?


已知:
* 方塊面積: $2\times 2 = 4$
* 圓面積: $1^2\times \pi = \pi$

因此, 隨機取點落在圓內的機率為 $\frac{\pi}{4}$.

因此, 若我們在 $[-1,1]\times[-1,1]$ 這方塊裡隨機取 $N$ 個點, 發現有 $M$ 個點落在圓內, 那我們知道
$$
\frac{M}{N} \approx \frac{\pi}{4} \,\to\,
\pi \approx 4\frac{M}{N}
$$
也就是說, 如果我們不知道 $\pi$ 的值, 我們可以用以上這方式來估算.

%%
%  \pi approximation

%%
N = 10^5;               % number of samples
M = 0;                  % counting initialize
for ii=1:N
    x = rand(2,1)*2-1;  % [-1,1]隨機取點
    if norm(x)<= 1
        M = M+1;
    end
end
disp(['pi 大約是 ', num2str(4*M/N)])

### Exercise 2

試以此方法估算橢圓內部面積:
$$
x^2 + 4y^2 \le 1.
$$

---

### Exercise 3

試寫一程式估算 $N$ 維單位球體積:
$$
\sum^N_{i=1} x_i^2 \le 1, \quad x=(x_1,\cdots,x_N)\in R^N
$$
此程式 `input` 為 $N$, 球體維度.

#### Reference:
* [Volume of an n-ball](https://en.wikipedia.org/wiki/Volume_of_an_n-ball)

---