How to solve $\int_{0}^{T} \int_{0}^{T} K(x, y) d z_{x} d z_{y}$

簡單的情形

如果 $K(x, y)$ 是這種 form $\sum_{i} f_{i}(x) \cdot g_{i}(y)$ ，則

$\begin{aligned}& \int_{0}^{T} \int_{0}^{T} K(x, y) d z_{x} d z_{y} \\& =\sum_{i} \int_{0}^{T} \int_{0}^{T} f_{i}(x) g_{i}(y) d z_{x} d z_{y} \\& =\sum_{i} \int_{0}^{T} f_{i}(x) d z_{x} \int_{0}^{T} g_{i}(y) d z_{y} \\& =\mu+\sum_{i}\left(\int_{0}^{T} f_{i}(x) d z_{x} \int_{0}^{T} g_{i}(y) d z_{y}-\int_{0}^{T} f_{i}(s) g_{i}(s) d s\right) \\& \mu \equiv \int_{0}^{T} K(s, s) d s=\sum_{i} \int_{0}^{T} f_{i}(s) g_{i}(s) d s\end{aligned}$

為一堆有相關性的 normal random variable 相乘再相加。如果 $f_{i}(x)$ 和 $g_{i}(x)$ 一樣，這樣計算積分或模擬機率變數就可以少一半。

如果 $K(x, y)$ 不是這種好用的 form ，就用這種好用的 form 去逼近。

誤差控制

取 $N$ 項計算的答案誤差的機率變數為：

$\sum_{i>N}\left(\int_{0}^{T} f_{i}(x) d z_{x} \int_{0}^{T} g_{i}(y) d z_{y}-\int_{0}^{T} f_{i}(s) g_{i}(s) d s\right)$

它的期望值為零。它的變異數為

$\begin{aligned}& E\left(\left(\sum_{i>N}\left(\int_{0}^{T} f_{i}(x) d z_{x} \int_{0}^{T} g_{i}(y) d z_{y}-\int_{0}^{T} f_{i}(s) g_{i}(s) d s\right)\right)^{2}\right) \\& =\sum_{i>N, j>N} E\left(\int_{0}^{T} \int_{0}^{T} \int_{0}^{T} \int_{0}^{T} f_{i}\left(x_{1}\right) g_{i}\left(x_{2}\right) f_{j}\left(x_{3}\right) g_{j}\left(x_{4}\right) d z_{x_{1}} d z_{x_{2}} d z_{x_{3}} d z_{x_{4}}\right) \\& -2 \sum_{i>N, j>N} E\left(\int_{0}^{T} f_{i}(x) d z_{x} \int_{0}^{T} g_{i}(y) d z_{y}\right) \int_{0}^{T} f_{j}(s) g_{j}(s) d s+\sum_{i>N, j>N}\left(\int_{0}^{T} f_{i}(s) g_{i}(s) d s\right)\left(\int_{0}^{T} f_{j}(s) g_{j}(s) d s\right) \\& =\sum_{i>N, j>N} E\left(\int_{0}^{T} \int_{0}^{T} \int_{0}^{T} \int_{0}^{T} f_{i}\left(x_{1}\right) g_{i}\left(x_{2}\right) f_{j}\left(x_{3}\right) g_{j}\left(x_{4}\right) d z_{x_{1}} d z_{x_{2}} d z_{x_{3}} d z_{x_{4}}\right)-\left(\sum_{i>N} \int_{0}^{T} f_{i}(s) g_{i}(s) d s\right)^{2}\end{aligned}$

summation 裡面第一項期望值不是零的項只有當 $x_{1}=x_{2},x_{3}=x_{4}$ 或 $x_{1}=x_{3}, x_{2}=x_{4}$ 或 $x_{1}=x_{4}, x_{2}=x_{3}$等於：

$\begin{aligned}& \int_{0}^{T} \int_{0}^{T} f_{i}\left(x_{2}\right) g_{i}\left(x_{2}\right) f_{j}\left(x_{4}\right) g_{j}\left(x_{4}\right) d x_{2} d x_{4} +\int_{0}^{T} \int_{0}^{T} f_{i}\left(x_{3}\right) g_{i}\left(x_{4}\right) f_{j}\left(x_{3}\right) g_{j}\left(x_{4}\right) d x_{3} d x_{4} +\int_{0}^{T} \int_{0}^{T} f_{i}\left(x_{4}\right) g_{i}\left(x_{3}\right) f_{j}\left(x_{3}\right) g_{j}\left(x_{4}\right) d x_{3} d x_{4} \\& =\left(\int_{0}^{T} f_{i}(s) g_{i}(s) d s\right)\left(\int_{0}^{T} f_{j}(s) g_{j}(s) d s\right)+\left(\int_{0}^{T} f_{i}(s) f_{j}(s) d s\right)\left(\int_{0}^{T} g_{i}(s) g_{j}(s) d s\right)+\left(\int_{0}^{T} f_{i}(s) g_{j}(s) d s\right)\left(\int_{0}^{T} f_{j}(s) g_{i}(s) d s\right)\end{aligned}$

所以誤差機率變數的變異數為 ：

$\begin{aligned}& \sum_{i>N, j>N}\left(\int_{0}^{T} f_{i}(s) g_{i}(s) d s\right)\left(\int_{0}^{T} f_{j}(s) g_{j}(s) d s\right) \\+ & \sum_{i>N, j>N}\left(\int_{0}^{T} f_{i}(s) f_{j}(s) d s \cdot \int_{0}^{T} g_{i}(s) g_{j}(s) d s+\int_{0}^{T} f_{i}(s) g_{j}(s) d s \cdot \int_{0}^{T} f_{j}(s) g_{i}(s) d s\right) \\- & \left(\sum_{i>N} \int_{0}^{T} f_{i}(s) g_{i}(s) d s\right)^{2} \\= & \sum_{i>N, j>N}\left(\int_{0}^{T} f_{i}(s) f_{j}(s) d s \cdot \int_{0}^{T} g_{i}(s) g_{j}(s) d s+\int_{0}^{T} f_{i}(s) g_{j}(s) d s \cdot \int_{0}^{T} f_{j}(s) g_{i}(s) d s\right)\end{aligned}$

以 $\eta_{i}$ 表示 $\sqrt{\int_{0}^{T} f_{i}(s)^{2} d s}$ 和 $\sqrt{\int_{0}^{T} g_{i}(s)^{2} d s}$ 的 upper bound，通常 $\eta_{i}$ 為一個隨著 $i$ 迅速遞減的數列，則：

$\begin{aligned}& \left| \sum_{i>N, j>N}\left(\int_{0}^{T} f_{i}(s) f_{j}(s) d s \cdot \int_{0}^{T} g_{i}(s) g_{j}(s) d s+\int_{0}^{T} f_{i}(s) g_{j}(s) d s \cdot \int_{0}^{T} f_{j}(s) g_{i}(s) d s\right) \right| \\& \leq \sum_{i>N, j>N}\left|\int_{0}^{T} f_{i}(s) f_{j}(s) d s \cdot \int_{0}^{T} g_{i}(s) g_{j}(s) d s\right| + \left|\int_{0}^{T} f_{i}(s) g_{j}(s) d s \cdot \int_{0}^{T} f_{j}(s) g_{i}(s) d s \right| \\& \leq \sum_{i>N, j>N} \sqrt{\int_{0}^{T} f_{i}(s)^{2} d s} \sqrt{\int_{0}^{T} f_{j}(s)^{2} d s} \sqrt{\int_{0}^{T} g_{i}(s)^{2} d s} \sqrt{\int_{0}^{T} g_{j}(s)^{2} d s}+\sqrt{\int_{0}^{T} f_{i}(s)^{2} d s} \sqrt{\int_{0}^{T} g_{j}(s)^{2} d s} \sqrt{\int_{0}^{T} f_{j}(s)^{2} d s} \sqrt{\int_{0}^{T} g_{i}(s)^{2} d s} \\& \leq \sum_{i>N, j>N} \eta_{i} \eta_{j} \eta_{i} \eta_{j}+\eta_{i} \eta_{j} \eta_{j} \eta_{i}=2 \sum_{i>N, j>N} \eta_{i}^{2} \eta_{j}^{2}=2\left(\sum_{i>N} \eta_{i}^2\right)^{2}\end{aligned}$

所以只要要求 $N$ 滿足 $\sqrt{2} \sum_{i>N} \eta_{i}^{2} \leq \epsilon$ ，就可以做到誤差的標準差小於 $\epsilon$ 。

例子

計算和模擬 $\int_{0}^{T} \int_{0}^{T} e^{-\beta|x-y|} d z_{x} d z_{y}$ 因為我們只在乎 $-T \leq a \leq T$ 的行為，所以把 $e^{-\beta|a|}$ 想成一個週期為 $2 T$ 由 $-T \leq a \leq T$ 所展開的函數。又因為它是偶函數，所以傅立葉展開為 $e^{-\beta|a|}=\frac{1}{2} a_{0}+\sum_{n} a_{n} \cos (n \omega a)$ 其中 $\omega=\frac{\pi}{T}$

$\begin{aligned}& a_{0}=\frac{1}{T} \int_{-T}^{T} e^{-\beta|a|} d a=\frac{2}{T} \int_{0}^{T} e^{-\beta a} d a=\frac{2\left(1-e^{-\beta T}\right)}{\beta T} \\& a_{n}=\frac{1}{T} \int_{-T}^{T} e^{-\beta|a|} \cos (n \omega a) d a=\frac{2}{T} \int_{0}^{T} e^{-\beta a} \cos (n \omega a) d a \\& =\frac{2}{T} \cdot \frac{\beta-e^{-\beta T}(\beta \cos (n \omega T)-n \omega \sin (n \omega T))}{\beta^{2}+n^{2} \omega^{2}}=\frac{2 \beta\left(1+e^{-\beta T}(-1)^{n+1}\right)}{\left(\beta^{2}+n^{2} \omega^{2}\right) T}\end{aligned}$

所以

$\begin{aligned}& \mu=T \\& e^{-\beta|x-y|}=\frac{1}{2} a_{0}+\sum_{n} a_{n} \cos (n \omega(x-y)) \\& =\frac{1}{2} a_{0}+\sum_{n} a_{n}(\cos (n \omega x) \cos (n \omega y)+\sin (n \omega x) \sin (n \omega y)) \\& =\frac{1}{2} a_{0}+\sum_{n} a_{n} \cos (n \omega x) \cos (n \omega y)+\sum_{n} a_{n} \sin (n \omega x) \sin (n \omega y)\end{aligned}$

所以可以設定

$\begin{aligned}f_{0}(x) & \equiv \sqrt{\frac{1-e^{-\beta T}}{\beta T}} \\g_{0}(y) & \equiv f_{0}(y)\end{aligned}$

其他的 $f_{i}(x)$ 和 $g_{i}(x)$ 為

$\begin{aligned}& f_{2 n-1}(x) \equiv \sqrt{a_{n}} \cos (n \omega x) \\& g_{2 n-1}(y) \equiv f_{2 n-1}(y) \\& f_{2 n}(x) \equiv \sqrt{a_{n}} \sin (n \omega x) \\& g_{2 n}(y) \equiv f_{2 n}(y)\end{aligned}$

所以

$\begin{aligned}& \int_{0}^{T} \int_{0}^{T} e^{-\beta|x-y|} d z_{x} d z_{y} \\& =T+\sum_{i} \int_{0}^{T} f_{i}(x) d z_{x} \int_{0}^{T} g_{i}(y) d z_{y}-\int_{0}^{T} f_{i}(s) g_{i}(s) d s \\& =T+\sum_{i}\left(\int_{0}^{T} f_{i}(x) d z_{x}\right)^{2}-\int_{0}^{T} f_{i}(s)^{2} d s \\& =T+\left(\int_{0}^{T} \sqrt{\frac{1-e^{-\beta T}}{\beta T}} d z_{x}\right)^{2}-\frac{1-e^{-\beta T}}{\beta} \\& +\sum_{n=1}^{\infty}\left(\int_{0}^{T} \sqrt{a_{n}} \cos (n \omega x) d z_{x}\right)^{2}+\left(\int_{0}^{T} \sqrt{a_{n}} \sin (n \omega x) d z_{x}\right)^{2}-a_{n} T \\& =T-\frac{1-e^{-\beta T}}{\beta}+\frac{1-e^{-\beta T}}{\beta T}\left(\int_{0}^{T} d z_{x}\right)^{2} \\& +\sum_{n=1}^{\infty} a_{n}\left(\left(\int_{0}^{T} \cos (n \omega x) d z_{x}\right)^{2}+\left(\int_{0}^{T} \sin (n \omega x) d z_{x}\right)^{2}-T\right)\end{aligned}$

定義

$\begin{aligned}& Z_{0} \equiv \int_{0}^{T} d z_{x} \\& Z_{2 n-1} \equiv \int_{0}^{T} \cos (n \omega x) d z_{x} \\& Z_{2 n} \equiv \int_{0}^{T} \sin (n \omega x) d z_{x}\end{aligned}$

模擬這群相關 normal 機率變數後即可用

$T\left(1-\sum_{n=1}^{N} a_{n}\right)-\frac{1-e^{-\beta T}}{\beta}+\frac{1-e^{-\beta T}}{\beta T} Z_{0}^{2}+\sum_{n=1}^{N} a_{n}\left(Z_{2 n-1}^{2}+Z_{2 n}^{2}\right)$

模擬

$\int_{0}^{T} \int_{0}^{T} e^{-\beta|x-y|} d z_{x} d z_{y}$

誤差控制。 $\eta_{i}=\sqrt{\frac{a_{i} T}{2}}$ 所以 $\sqrt{2} \sum_{i>N} \eta_{i}^{2} \leq \epsilon$ 為 $\frac{T}{\sqrt{2}} \sum_{n>N} a_{n} \leq \epsilon$

$\begin{aligned}& \frac{T}{\sqrt{2}} \sum_{n>N} a_{n} \\& =\frac{T}{\sqrt{2}} \sum_{n>N} \frac{2 \beta\left(1+e^{-\beta T}(-1)^{n+1}\right)}{\left(\beta^{2}+n^{2} \omega^{2}\right) T} \\& =\sqrt{2} \beta\left(\left(1-e^{-\beta T}\right) \sum_{n>N, \text { even }} \frac{1}{\beta^{2}+n^{2} \omega^{2}}+\left(1+e^{-\beta T}\right) \sum_{n>N, \text { odd }} \frac{1}{\beta^{2}+n^{2} \omega^{2}}\right) \\& <\sqrt{2} \beta\left(\left(1-e^{-\beta T}\right) \sum_{n>N, \text { even }} \frac{1}{2} \int_{n-2}^{n} \frac{1}{s^{2} \omega^{2}} d s+\left(1+e^{-\beta T}\right) \sum_{n>N, \text { odd }} \frac{1}{2} \int_{n-2}^{n} \frac{1}{s^{2} \omega^{2}} d s\right) \\& <\frac{\beta}{\sqrt{2}}\left(\left(1-e^{-\beta T}\right) \int_{N-1}^{\infty} \frac{1}{s^{2} \omega^{2}} d s+\left(1+e^{-\beta T}\right) \int_{N-1}^{\infty} \frac{1}{s^{2} \omega^{2}} d s\right) \\& =\frac{\sqrt{2} \beta}{\omega^{2}(N-1)}\end{aligned}$

所以如果 $N$ 滿足 $\frac{\sqrt{2} \beta}{\omega^{2}(N-1)}<\epsilon$ 則可以讓誤差機率變數的標準差小於 $\epsilon$ 。