# 積分方法第二部
{%hackmd @themes/orangeheart %}
<style>
.likecoin-button {
position: relative;
width: 100%;
max-width: 485px;
max-height: 240px;
margin: 0 auto;
}
.likecoin-button > div {
padding-top: 49.48454%;
}
.likecoin-button > iframe {
position: absolute;
top: 0;
left: 0;
width: 100%;
height: 100%;
}
</style>
###### tags: `Numerical Methods`

在[積分方法第一部](https://hackmd.io/@yueswater/integration_I)中，我們提到了最直覺的中點法、梯形法與辛普森法 $\cdots$ 等等；且最後提及了有關數值積分的近似與計算方法，透過取節點與權重的方式進行積分計算。而在經濟領域上，我們常常使用積分計算期望值。根據不自覺統計學家法則(law of the unconscious statistician, LOTUS)，給定一個隨機變數的機率密度函數 $f_X$，如果要計算期望值，計算方式為

$$
\begin{gathered}
&\mathbb{E}[g(X)] = \sum_{x} g(x) f_{X}(x) \\
&\mathbb{E}[g(X)] = \int_{-\infty}^{\infty} g(x) f_{X}(x) \mathrm{d} x
\end{gathered}
$$

那上述的方法要如何使用數值積分的方式呈現呢？我們就來看看。

## Monte Carlo 積分法
考慮一個函數 $g(X)$ 的積分，其中 $X$ 為一個隨機變數

$$
I = \int^a_b g(X) dx
$$

其又可以表達為

$$
(b - a) \int^a_b g(X) \frac{1}{b - a} dx = (b - a) \mathbb{E}[g(X)]
$$

其中 $\frac{1}{b - a}$ 是一個服從上界為 $b$、下界為 $a$ 之均勻分配的 PDF。在數值分析領域上， $(b - a)$ 常常被稱作 **volume**。因此，上述的積分可以被看作是一個服從均勻分配的期望值。更精確地來說，我們可以將上述積分改為下列形式：
$$
I = \int_a^b \left[\frac{g(X)}{p(x)}\right] p(x)dx = (b-a) \int_a^b g(X) p(x) dx = (b-a) E[g(X)]
$$

其中 $p(x) = \frac{1}{b - a}$。而這個方法與所謂的**重要抽樣(importance sampling)** 類似。在實務上，我們通常會將上、下界轉為 $[0, 1]$，並運用轉換規則與 Jacobian 項的方式進行積分。

$$\begin{aligned}
I = \int_a^b g(X) dx = \int_0^1 f(t) dt = (1-0) \int_0^1 f(t) \frac{1}{1-0} dt = \mathbb{E}[f(t)].
\end{aligned}$$

因此我們要找到一個轉換規則 $x = \rho(t)$，使得 $\rho^{-1}(b) = 1$，以及 $\rho^{-1}(a) = 0$，並得到 $f(t) = g(\rho(t))\rho'(t)$，其中 Jacobian 項為 $\rho'(t)$。故可以將上述的積分轉換為：

$$
I = (b-a) \mathbb{E}[g(X)] = \mathbb{E}[f(t)]
$$
其中隨機變數 $x \sim U(a, b)$，$t \sim U(0, 1)$。以下提供一些轉換規則。[^1]

$$ x 值域 $$ | 轉換規則 | $t$ 值域 | Jacobian 項
:---: | :---: | :---: | :---:
$$[a, b]$$ | $$x = a + (b-a)t$$ | $$[0,1]$$ | $$b - a$$
$$[-\infty, \infty]$$ | $$x = \frac{2t-1}{t-t^2}$$ | $$[0,1]$$ | $$\frac{2t^2 - 2t+1}{(t^2 -t)^2}$$
$$[a, \infty]$$ | $$x = a + \frac{t}{1-t}$$ | $$[0,1]$$ | $$\frac{1}{(t-1)^2}$$
$$[-\infty, b]$$ | $$x = b + \frac{t-1}{t}$$ | $$[0,1]$$ | $$\frac{1}{t^2}$$

那要怎麼計算 $\mathbb{E}[f(t)]$ 呢？其實只要根據大數法則(law of large number, LLN)，抽取足夠大的樣本 $t$，樣本平均 $f(t_i)$ 便可以機率收斂到母體平均，即

$$
\begin{aligned}
\mathbb{E}[f(t)] \approx \bar{f}_n = \frac{1}{n} \sum_{i=1}^n f(t_i).
\end{aligned}$$

因此，

$$\begin{aligned}
I = \int_a^b g(X) dx = \int_0^1 f(t)dt \approx \bar{f}_n = \frac{1}{n} \sum_{i=1}^n f(t_i),
\end{aligned}$$

其中 $\{t_i\}$, $i=1,\cdots,n$ 服從 $U(0, 1)$。

:::success
實作：用 Monte Carlo 積分法估計平均。
:::

考慮以下的積分：
$$\begin{aligned}
I = \int_{-\infty}^\infty \exp\left(-\frac{1}{3}x^2\right)\sqrt{1+x^2} dx.
\end{aligned}$$

## 引用套件
using Random, Plots, Interact, WebIO, Random, QuadGK, LaTeXStrings

## 設定樣本數
n_min = 1; n_max = 1000

@manipulate for n=n_min:n_max
    g(x) = exp(-x^2/3) * sqrt(1+x^2) ## 設定被積函數(integrand)

    function monteInt(g::Function, n::Int)
        trans(x) = (2x - 1) / (x-x^2)
        jacobian(x) = (2x^2-2x+1) / (x^2-x)^2
        f(x) = g(trans(x)) * jacobian(x)
        sample = rand(n)
        result = (1 / n) * sum(f.(sample))
        return result
    end
answer = monteInt(g, 100000)
    
## 繪圖
nofn = n_max - n_min + 1  ## 設定儲存陣列長度

ans_monte = zeros(nofn) ## 建立抽樣儲存陣列
true_ans = zeros(nofn) ## 建立實際積分數值陣列

for n in n_min:n_max
    ans_monte[n-n_min+1] = monteInt(g, n) ## 抽樣計算
    true_ans[n-n_min+1] = quadgk(x->exp(-x^2/3) * sqrt(1+x^2), -Inf, Inf)[1] ##使用 QuadGK 套件計算真實數值
end

## 繪製積分收斂情況
plot(n_min:n_max,
     ans_monte,
     xlim = (n_min, n),
     ylim = (0,10),
     xlabel = L"$n$ of sample",
     ylabel = "Value",
     label = "Monte Carlo Integration"
)

## 繪製真實數值(為一常數)
plot!(n_min:n_max, 
      true_ans,
      linewidth = 1.5,
      label = L"\texttt{quadgk()}"
)
end

<iframe width="900" height="800" frameborder="0" scrolling="no" src="//plotly.com/~xiaolong70701/26.embed"></iframe>

可以看到，儘管我們用了大約 $1000$ 個樣本點進行抽樣，其收斂速度還是很慢，在術語上我們就說它不夠有「效率」(efficient)。或許還有其他方法可以讓它收斂速度更快！

[^1]:由國立臺灣大學經濟學研究所學長陳明鴻整理。

<div class="likecoin-embed likecoin-button">
<div></div>
<iframe scrolling="no" frameborder="0" src="https://button.like.co/in/embed/xiaolong70701/button?referrer=hackmd.io"></iframe>
</div>