# 蒙地卡羅模擬之資料分析方法

:::info
本文主旨為介紹常見蒙地卡羅之資料處理方法
目的為增加蒙地卡羅所產生資料之可靠性與獲得更好的分析結果
文中會提到的技術有自相關函數、分箱及拔靴法
本文所使用之程式為 Python, Julia
:::
:::warning
本文中的程式僅作為示範，實際分析請自行撰寫以得到更正確的結果或更好的效能
:::

## 1. 自相關函數(Autocorrelation function)

由於馬可夫鍊蒙地卡羅方法的特性，下一個狀態完全取決於上一個狀態，因此相鄰的兩個狀態之間會有很強的相關性，但在做統計時數據需要有獨立性，要解決這個問題我們可以計算求得之觀測量的自相關函數，並透過積分該函數來取得這筆資料的自相關時間，若資料中的兩個數據之時間(步數)間隔大於2倍以上的自相關時間，我們就可以說這兩個數據的相關性是可忽略的，這個操作稱作去相關(decorrelate)

自相關函數定義為
\begin{equation}
f_A(\tau) = {\Gamma (\tau ) \over \Gamma (0)}
\end{equation}
其中
\begin{equation}
\Gamma (\tau)={1 \over N- \tau} \sum^{N- \tau-1}_{i=0}[(O[i]- \bar{O}) \times (O[i+ \tau] - \bar{O})]
\end{equation}

其中$O$為觀測量，$N$為總觀測數，$i$表示第$i$個觀測量，而$\tau$為兩觀測量的間距

>[!Note]
> 在統計學中，共變異數(協方差)的定義為
> \begin{equation}
> \mathrm{cov}(X,Y)=\langle (X-\bar{X})(Y-\bar{Y}) \rangle
> \end{equation}
> 其中$\bar{X}$表示$X$的平均值，共變異數表示了兩個變數間線性相關的趨勢，若其中$Y=X$即為變異數
> ::: warning
> 注意：共變異數為0不表示兩變數無關，只能表示兩變數沒有線性關係，其有可能有非線性關係，如$\mathrm{cov}(X,X^2)=0$
> :::
> 在隨機過程$\{X_t\}$中，自共變異數定義為
> \begin{equation}
> K_{XX}(t_1,t_2)=\mathrm{cov}(X_{t_1},X_{t_2})=\langle (X_{t_1}-\overline{X_{t_1}})(X_{t_2}-\overline{X_{t_2}}) \rangle
> \end{equation}
> 若$\{X_t\}$為弱平穩(WSS)過程(平均值與自共變異數不隨時間變化)則有
> \begin{equation}
> \overline{X_{t_1}}=\overline{X_{t_2}}=\overline{X}
> \end{equation}
> 及
> \begin{equation}
> K_{XX}(t_1,t_2)=K_{XX}(t_2-t_1,0) \equiv K_{XX}(t_2-t_1) = K_{XX}(\tau)
> \end{equation}
> 而熱平衡後的靜態(目標分布不隨時間改變)蒙地卡羅通常為弱平穩，將弱平穩下的自共變異數歸一化即得到上面定義之自相關函數$f_A(\tau)$

一般情況下，自相關函數呈指數衰減(exponential decay)
\begin{equation}
f_A(\tau) ∼ e^{-\tau / \tau_A}
\end{equation}
其中$\tau_A$為自相關時間
自相關時間可以透過積分自相關函數算出
\begin{equation}
\int^\tau_0 d\tau'f_A(\tau') ∼ \int^\tau_0 d\tau' e^{-\tau / \tau_A} ∼ \tau_A \ \text{if} \ \tau \ ≫ \tau_A
\end{equation}
在數值計算中，自相關函數之積分定義為
\begin{equation}
\Theta (\tau)={1 \over 2}+\sum^{\tau}_{\tau'=1}f_A(\tau')
\end{equation}
當$\tau$足夠大至$\Theta$之值穩定時，我們即可得到自相關時間$\tau_A=\Theta$，在實際運算中，我們將$\tau$選在$\Gamma(\tau)$無法被正確計算前的點($\Gamma(\tau)$不再滿足指數衰減)

>[!Note]
>我們有弱平穩自共變異數的定義
>\begin{equation}
K(\tau)=\langle (O_i- \bar{O})(O_{i+ \tau} - \bar{O})\rangle
\end{equation}
>可以注意到$K(\tau)=K(-\tau)$，$-\infty < \tau<\infty$。因此自相關函數的積分
>\begin{align}
\Theta=\sum^{\infty}_{\tau'=-\infty}f_A(\tau') &= f_A(0)+2\sum^{\infty}_{\tau'=1}f_A(\tau')\\
&= 2 \left( {1 \over 2}f_A(0)+\sum^{\infty}_{\tau'=1}f_A(\tau') \right)
\end{align}
>其中$f_A(0)=1$。而我們只考慮$\tau >\tau' \ge 0$，因此有上述積分形式之定義

:::warning
不同的觀測量有不同的自相關時間，因此在處理數據時應分開處理所有觀測量，或直接以自相關時間最大的觀測量為標準進行數據處理
:::

在做去相關時要求至少間隔2倍以上的自相關時間是由於計算自相關時間會有很大的誤差，其誤差可以用Madras-Sokal公式估算
\begin{equation}
(\Delta \Theta(\tau))^2={4\tau+2 \over N}\Theta(\tau)^2
\end{equation}
由於篇幅較長，此文章中僅對自相關函數之計算做示範，並不包含誤差估算

### Python 分析實作

import numpy as np
#定義Γ函數
def Gamma(A: np.ndarray, tau: int) -> float:
    res = 0.0
    N = len(A)
    A_mean = np.mean(A)
    for i in range((N-tau)):
        res += ((A[i]-A_mean)*(A[i+tau]-A_mean))
    return res/(N-tau)

#計算自相關函數(用來判斷tau需要多大)
def autofunc(A: np.ndarray) -> np.ndarray:
    N = len(A)
    res = np.zeros(N)
    Gamma_0 = Gamma(A, 0)
    for i in range(N):
        res[i] = Gamma(A, i)/Gamma_0 # res[i] = f_A(i), 注意這裡包含 f_A(0)
    return res
        
#計算自相關時間
def autotime(autofunc: np.ndarray, tau: int) -> np.ndarray:
    return 0.5 + np.sum(autofunc[1:tau+1])

### Python 分析實作2(向量化版本)

import numpy as np
#向量化定義Γ函數
def Gamma_vec(A: np.ndarray, tau: int) -> float:
    N = len(A)
    A_mean = np.mean(A)
    return np.mean((A[:N-tau] - A_mean) * (A[tau:] - A_mean))

#向量化計算自相關函數
def autofunc_vec(A: np.ndarray) -> np.ndarray:
    N = len(A)
    Gamma_0 = Gamma_vec(A, 0)
    taus = np.arange(N)  # 計算 f_A(0) 至 f_A(N-1)
    return np.array([Gamma_vec(A, tau)/Gamma_0 for tau in taus])

### Julia 分析實作
:::warning
注意：Julia陣列編號由1開始
:::

using Statistics
#定義Γ函數
function Gamma(A::Vector, tau::Int)
    res = 0.0
    N = length(A)
    A_mean = mean(A)
    for i in 1:(N-tau)
        res += ((A[i]-A_mean)*(A[i+tau]-A_mean))
    end
    return res/(N-tau)
end

#計算自相關函數
function autofunc(A::Vector)
    N = length(A)
    res = zeros(N)
    Gamma_0 = Gamma(A, 0)
    for i in 1:N
        res[i] = Gamma(A, i)/Gamma_0 # res[i] = f_A(i), 注意這裡沒有 f_A(0)
    end
    return res
end

#計算自相關時間
function autotime(autofunc::Vector, tau::Int)
    return 0.5 + sum(autofunc[1:tau])
end

Julia本身擁有極為良好的編譯器，不用過度追求向量化

### 補充

>[!Tip]小技巧
>注意到離散化捲積有定義
>\begin{equation}
>y[n]=f[n] * g[n]=\sum^{M-1}_{m=0} f[n-m]g[m]
>\end{equation}
>對比$\Gamma$函數定義
>\begin{equation}
>\Gamma (\tau) = {1 \over N- \tau} \sum^{N- \tau-1}_{i=0}[(O[i]- \bar{O}) \times (O[i+ \tau] - \bar{O})]
>\end{equation}
>移項得
>\begin{equation}
>(N-\tau)\Gamma(\tau) = \sum^{N-\tau-1}_{i=0}O'[\tau+i]O'[i]
>\end{equation}
>其中$O'=O-\bar{O}$，可以看出其形式與離散化捲積類似，事實上我們僅需將左側$O'$做翻轉即與上述定義相同。而對於離散化捲積我們有傅立葉轉換
>\begin{equation}
>y[n]=f[n]*g[n] \leftrightarrow Y[f]=F[f]G[f]
>\end{equation}
>此兩函數在時域$(n)$做捲積相當於此兩函數的傅立葉轉換在頻域$(f)$相乘，因此我們有
>\begin{equation}
>y[n]=\mathcal{F^{-1}}(\mathcal{F}(f[n])\mathcal{F}(g[n]))
>\end{equation}
>而對於翻轉之$O'$我們可以將其想為$f(-x)$，而我們有
>\begin{equation}
>G(k)=\int^\infty_{-\infty}f(-x)e^{-2\pi ikx}dx=\int^\infty_{-\infty}f(u)e^{2\pi ikx}du=\overline{F(k)}
>\end{equation}
>因此自相關函數的計算簡化為
>\begin{equation}
>f_A(\tau)={1\over N \sigma_A^2}\mathcal{F^{-1}}(|\mathcal{F}(O')|^2)
>\end{equation}
>利用此方法可以將自相關函數計算的時間複雜度$\mathcal{O}(N^2)$減少至$\mathcal{O}(N\mathit{log}N)$，在資料量龐大時可以有效減少計算時間。
>:::warning
>要注意由於FFT方法有週期性邊界條件，且分母為$N$而非$N-\tau$，因此此方法在接近數據尾端時及$\tau$較大時有誤，如果有較大的資料量且$\tau_A \ll N$則偏差可以忽略，但仍建議使用補零法(Zero-padding)來避免此問題
>:::
>```julia=
>using FFTW
># 此程式未使用補零
>function autofunc_fft_large(data::Vector{Float64})
> N = length(data)
> x = data .- mean(data) # O - O_bar
> fx = fft(x)
> acf = ifft(abs2.(fx)) ./ N # FFT(O_rev) = FFT(O)*, normalize
> return real(acf[2:end] ./ acf[1]) # 自相關，跳過 tau=0
>end
>```

## 2. 分箱(Binning)
在做完去相關後，我們的數據量會大幅減少，會有許多數據被浪費掉，因此我們可以透過分箱的方式充分利用所有數據。分箱的方法為
* 若原資料$A$有自相關時間$\tau_A$，將連續的每$n\tau_A$個數據進行平均，其中$n \ge 2$
* 由上步得到新資料$O$，而$O$的數據量$N_O = {N_A \over n\tau_A}$

### Python 分析實作

import numpy as np
def bin(A: np.ndarray, ntau: int) -> np.ndarray:
    return A.reshape(-1, ntau).mean(axis=1)

### Julia 分析實作

using Statistics
function bin(A::Vector, ntau::Int)
    return vec(mean(reshape(A, ntau, :), dims=1))
end

## 3. 拔靴法(Bootstrap)

蒙地卡羅模擬通常需要大量計算資源。若缺乏高效算法或足夠的運算能力，且觀測量具有較大的自相關時間，則為了取得足夠的資料量，往往必須花費極長的計算時間。

為了在有限的資源下，更好的估計物理量的誤差以獲得更好的分析結果，除了再花時間增加更多的數據外，我們可以透過一些數學方法來對數據進行處理，常見的方法有拔靴法及刀切法(Jackknife)，這兩種方法都是對已有的資料進行重新抽樣，但運作原理不同，此文中僅對拔靴法進行示範

>[!Important]重點
>拔靴法能有效給出擁有更好標準差之資料，但並不會改變平均值，且即使此方法看似有數據造假之疑慮，但此方法並沒有憑空捏造數據，所有的數據皆由原始資料中產生。拔靴法僅透過增加數據量的方法改變資料的精確度(precision)----數據之間靠近的程度，而不會改變準確性(Accuracy)----數據與真實值接近的程度

拔靴法的操作流程如下
* 從總數據量$N$的原始資料$A$中均勻隨機且可重複的取出$N$個數據，其中$A=\{ a_1...a_N\}$，抽出後得$M=\{m_1...m_N\}$
* 重複上述步驟$\alpha N$次，其中$\alpha$通常為4~10間的整數 (抽樣次數的選擇非固定，視情況也可以改為如固定1000次等)
* 在上一步中，我們會得到$\alpha N$個資料集$\{M_1...M_{\alpha N}\}$，其中每個$M_i$皆含有$N$個數據
* 將每個$M_i$做平均得到$\alpha N$個數據$\{ \langle M_1 \rangle...\langle M_{\alpha N} \rangle \} \equiv G$
* 由$G$即可計算出拔靴後的標準差$\sigma_G$與平均值$\langle G \rangle$，應滿足$\langle G \rangle \approx \langle A \rangle$且$\sigma_G \lt \sigma_A$

>[!Note]
>若資料具有強關聯性或資料過少可能導致拔靴法結果不如預期，應先計算自相關函數確保拔靴法結果可靠性

其中$\alpha$的大小可以藉由持續加大$\alpha$看$\sigma_G$是否繼續變小，若沒有明顯差異則不需增加，減少計算量。若要在拔靴法中計算觀測量，則應利用步驟三之數據$M_i$進行計算得到$\alpha N$個觀測量$O_i$，並將每個$O_i$取平均得到有$\alpha N$個$\langle O_i \rangle$的資料集$G_O$，再進行第五步之分析

### Python 分析實作

def bootstrap(A: np.ndarray, alpha: int) -> np.ndarray:
    N = len(A)
    res = np.zeros(N*alpha)
    for i in range(N*alpha):
        r = np.random.choice(N, N)
        res[i]=np.mean(A[r])
    return np.array(res)

### Julia 分析實作

using Statistics
function bootstrap(A::Vector, alpha::Int)
    N = length(A)
    res = zeros(N*alpha)
    for i in 1:(N*alpha)
        res[i] = mean(rand(A, N))
    end
    return res
end

>[!Warning]注意
>如果資料為2維或更多，在做拔靴時應確保其他維度所抽出之數據皆源於同一筆資料，否則將破壞拔靴法的正確性
>例：實做中，我們通常有資料集
>\begin{equation}
>A(t)=\begin{pmatrix}
a_{1}(t_1) & a_{1}(t_2) & \cdots & a_{1}(t_n) \\
a_{2}(t_1) & a_{2}(t_2) & \cdots & a_{2}(t_n) \\
\vdots & \vdots & \ddots & \vdots\\
a_{m}(t_1) & \cdots & \cdots & a_{m}(t_n)
\end{pmatrix}
>\end{equation}
>其中下標$i$表示第$i$筆資料，每筆資料皆有$n$個數據$for \ t\in\{t_1...t_n\}$
>常見的做法為抽取時一次抽取整組數據$M_i=\{ a_{i}(t_1), a_{i}(t_2)...a_{i}(t_n)\}$

## 參考資料

[1] Hsiao Ho, _Lattice Studies of the Composite Higgs Model on $ \mathrm{Sp(4)} $ Gauge Theories._ (2024)
[2] C.-J. David Lin, Kenji Ogawa, Hiroshi Ohkic and Eigo Shintani, _Lattice study of infrared behaviour in $ \mathrm{SU(3)} $ gauge theory with twelve massless flavours_, JHEP08(2012)096