貝葉斯迴歸與高斯過程

§01學習重點

• 把第 2 章的「先驗→後驗」機制搬到迴歸權重上：親手算出權重後驗如何隨資料收縮，小樣本時如何被先驗拉住
• 說出貝葉斯預測分布與 OLS 點估計的本質差別：預測自帶變異數，而且離資料越遠、變異數越大
• 解釋「先驗＝正則化」這座橋：高斯先驗下的 MAP 估計就是 ridge 迴歸，懲罰強度由雜訊與先驗的變異數比值決定，並以數值驗證到浮點精度
• 完成從權重空間到函數空間的跳躍：說出高斯過程的定義（均值函數＋核函數）、為什麼它是「非參數式」模型
• 讀懂 GP 預測公式的兩個直覺：後驗均值是觀測值的核加權平均、後驗變異數是「先驗不確定性減去被資料解釋掉的部分」
• 用邊際似然這把自動奧卡姆剃刀調長度尺度，同時說出 \(O(N^3)\) 的計算帳單與可擴展 GP 的解法方向
• 描述 GP 代理定價模型的離線／線上工作流：為什麼它能替昂貴定價器建快取、Greeks 為什麼可以解析微分、預測帶爆寬時該做什麼

§02課程內容

一、把「先驗→後驗」搬到迴歸權重上

第 2 章練好的貝葉斯基本功——先驗信念、看資料、更新成後驗——到目前為止都用在「單一參數」身上。本章的第一步很簡單：把同一套機制原封不動搬到線性迴歸的權重向量上。

先看古典做法少了什麼。最小平方法（OLS）給你一組權重的點估計：一個數字，沒有附任何「這個數字多可靠」的資訊。實務上大家會補 t 值與信賴區間，但那套語言建立在第 1 章講過的統計建模文化假設上，而且信賴區間描述的是「重複抽樣下估計程序的行為」，不是「參數本身的機率分布」——這兩件事常被混為一談。貝葉斯迴歸的回答更直接：把權重當成隨機變數。給權重一個先驗分布 \(\theta \sim \mathcal{N}(\mu, \Sigma)\)，觀測到資料 \(\mathcal{D}=(X, \mathbf{y})\) 後，用貝葉斯定理算出權重的後驗分布。當雜訊是高斯 \(y_i = \theta^\top x_i + \epsilon_i\)、\(\epsilon_i \sim \mathcal{N}(0, \sigma_n^2)\)，而先驗也是高斯時，後驗有漂亮的閉式解：

逐項拆解：\(\Sigma^{-1}\) 是先驗共變異矩陣的反矩陣，叫精度矩陣——精度是變異數的倒數，「你原本有多確定」；\(\tfrac{1}{\sigma_n^2}X^\top X\) 是資料帶來的精度——資料越多（\(X\) 的列越多）、雜訊越小（\(\sigma_n^2\) 越小），這一項越大。後驗精度＝先驗精度＋資料精度，兩股確定性直接相加，所以資料進來後驗只會更窄、不會更寬。後驗均值 \(\mu'\) 則是把「先驗的意見 \(\Sigma^{-1}\mu\)」與「資料的意見 \(\tfrac{1}{\sigma_n^2}X^\top \mathbf{y}\)」各依精度加權平均：資料少時先驗說話大聲，資料多時資料輾壓先驗。

比喻： 這像美食指南的密探評鑑一家新餐廳。開幕之前，評估只能來自主廚的出身與菜單——這是先驗，而且是一個「範圍」：應該有水準，但沒把握。密探第一次造訪（少量資料）只會微調這個範圍，因為單獨一餐可能撞上廚房的好日子或壞日子；一季內造訪十次之後，評鑑幾乎完全由實吃決定，範圍也收得很窄。關鍵是指南自始至終記錄的都是範圍而非單一分數——要不要給星、給幾顆星，看的正是範圍的下緣（最差的一餐也不能失手），而不是中間值。

用一個金融場景實跑：估計某檔台股對大盤的 beta。先驗說「一般股票 beta 在 1 附近」（\(\mathcal{N}(1, 0.5^2)\)），真實 beta 是 1.3：

import numpy as np

rng = np.random.default_rng(3)

# 情境：估計某檔台股對大盤的 beta。真實 alpha=0.001、beta=1.3
true_theta = np.array([0.001, 1.3])
sigma_n = 0.012                      # 日報酬雜訊（殘差波動）

# 先驗：alpha 以 0 為中心、beta 以 1 為中心（「一般股票跟著大盤動」）
mu0 = np.array([0.0, 1.0])
Sigma0 = np.diag([0.02**2, 0.5**2])

def posterior(X, y):
    """貝葉斯線性迴歸權重後驗：N(mu_post, Sigma_post)"""
    A = np.linalg.inv(Sigma0) + (X.T @ X) / sigma_n**2
    Sigma_post = np.linalg.inv(A)
    mu_post = Sigma_post @ (np.linalg.inv(Sigma0) @ mu0 + (X.T @ y) / sigma_n**2)
    return mu_post, Sigma_post

# 模擬 250 個交易日的大盤報酬，逐步餵給後驗
x_mkt = rng.normal(0.0, 0.01, 250)
y_stk = true_theta[0] + true_theta[1] * x_mkt + rng.normal(0.0, sigma_n, 250)

for n in (5, 20, 250):
    X = np.column_stack([np.ones(n), x_mkt[:n]])
    mu_p, Sig_p = posterior(X, y_stk[:n])
    ols = np.linalg.lstsq(X, y_stk[:n], rcond=None)[0]
    print(f"n={n:>3d}｜beta 後驗 {mu_p[1]:.3f} ± {np.sqrt(Sig_p[1,1]):.3f}"
          f"｜OLS 點估計 {ols[1]:.3f}")
# 輸出：
# n=  5｜beta 後驗 1.388 ± 0.292｜OLS 點估計 1.594
# n= 20｜beta 後驗 1.226 ± 0.200｜OLS 點估計 1.270
# n=250｜beta 後驗 1.316 ± 0.074｜OLS 點估計 1.323

三行輸出講了三件事。只有 5 天資料時，OLS 被雜訊帶著跑，估出離譜的 1.594；後驗均值 1.388 被先驗往 1 拉回來，而且 ±0.292 的寬度老實承認「我還不太確定」。20 天時兩者接近，但後驗仍附帶寬度。250 天後，後驗 1.316 ± 0.074 與 OLS 的 1.323 幾乎重合——資料夠多時，貝葉斯與頻率學派殊途同歸；差別在資料少的時候，那正是金融裡最常見的處境：新上市的股票、剛換經理人的基金、沒有歷史的新商品。

二、預測分布：讓模型自己說「我多不確定」

估權重通常不是終點，預測才是。貝葉斯預測的做法跟點估計有根本差異：不挑出一組最好的權重來預測，而是讓所有可能的權重一起預測，再依後驗機率加權平均。在新輸入 \(x_*\) 上，模型輸出 \(f_* := f(x_*)\) 的預測分布是：

逐項拆解：積分符號做的事就是「加權平均」——每一組權重 \(\theta\) 給出自己的預測 \(p(f_*\mid x_*,\theta)\)，權重是它的後驗機率 \(p(\theta\mid\mathcal{D})\)，積分把權重這個變數整個邊緣化掉，最終的預測不依賴任何特定一組權重。高斯情形下答案還是高斯：預測均值 \(x_*^\top\mu'\) 就是用後驗均值做預測；預測變異數 \(x_*^\top\Sigma' x_*\) 把「權重的不確定性」透過輸入 \(x_*\) 傳導到輸出上——\(x_*\) 落在資料密集、後驗掐得緊的方向，變異數小；落在資料稀疏的方向，變異數大。若要預測含雜訊的觀測值 \(y_*\)，再加上一項 \(\sigma_n^2\)。

這就是貝葉斯方法在金融裡最值錢的性質：每個預測都自帶一支「誤差條」，而且誤差條的寬度隨輸入位置變化。風險管理者拿到的不是「這檔債券的違約率是 2.1%」，而是「2.1%，上下界是 1.4% 到 3.0%」——後者才撐得起部位限額與資本計提的決策。第 1 章說過機器學習用樣本外表現當裁判；貝葉斯再往前一步，讓模型在給答案的同時申報自己的可信度。

三、先驗＝正則化：ridge 的貝葉斯身世

如果你不想要整個後驗分布、只想要一組「最好的權重」，貝葉斯框架也有對應的點估計：取後驗分布的眾數，叫 MAP 估計（maximum a posteriori）。高斯後驗的眾數就是均值，所以 MAP 解等於上一節的 \(\mu'\)。有意思的是把 MAP 的最佳化問題攤開來看。取對數後驗（對數把連乘變連加，不改變最大值位置），先驗取零均值等向高斯 \(\theta \sim \mathcal{N}(0, \sigma_\theta^2 I)\)：

逐項拆解：第一項 \(\|\mathbf{y}-X\theta\|_2^2\) 是殘差平方和——對數概似取負號後剩下的東西，「模型貼資料貼多緊」；第二項 \(\lambda\|\theta\|_2^2\) 來自對數先驗——高斯先驗的對數就是負的權重平方和乘上常數，「權重離 0 多遠就罰多重」。這正是 ridge 迴歸的損失函數。懲罰強度 \(\lambda\) 不是憑空調出來的超參數，而是兩個變異數的比值：雜訊越大（資料越不可信）罰越重，先驗越寬（你越沒有立場堅持權重要小）罰越輕。極限 \(\lambda \to 0\) 對應「先驗寬到等於沒有意見」，MAP 退化回最大概似估計，也就是 OLS。

比喻： 高斯先驗像把每個迴歸係數用一條彈簧拴在原點。資料想把係數拉到它喜歡的位置，就得出力拉伸彈簧；彈簧越硬（先驗變異數越小、\(\lambda\) 越大），資料要拉動係數就越費勁，係數就越乖乖待在 0 附近。ridge 迴歸的「懲罰」不是什麼抽象的數學裝置——它就是這條彈簧的回復力。而 LASSO 只是換了一種彈簧（拉普拉斯先驗），特性是乾脆把拉不太動的係數直接釘死在 0。

這座橋值得記牢，因為它雙向通行。往一邊走：任何正則化手法都可以問「它隱含什麼先驗」，第 1 章說貝葉斯方法「天生內建防過擬合機制」，機制就在這——先驗收縮了模型的有效複雜度。往另一邊走：第 4 章講前饋神經網路時，正則化會以頻率學派的面貌重新登場（權重衰減、dropout），你手上已經有一副貝葉斯之眼可以翻譯它們。作業二會請你把「ridge 解＝MAP 解」驗證到浮點數精度。

四、從權重空間跳進函數空間：高斯過程登場

到目前為止，我們都在權重空間工作：先選定函數形式（線性），再對權重放先驗。但這有個彆扭之處——你真正關心的往往是函數本身長什麼樣（平滑嗎？擺動多大？），而「對權重的先驗」只能間接、迂迴地表達這些直覺。何況函數形式選錯，權重估得再好也沒用。

換個角度想。在貝葉斯線性迴歸裡，權重的先驗其實已經誘導出函數值的分布：任取兩個輸入點 \(x_i, x_j\)，函數值 \(f(x_i), f(x_j)\) 都是權重的線性組合，所以它們服從聯合高斯分布，共變異數由先驗決定。既然最後用到的只是「函數值之間的共變異結構」，何不跳過權重，直接對函數放先驗？這一步就是高斯過程（Gaussian process, GP）：

逐項拆解：\(\mu(\cdot)\) 是均值函數——還沒看資料前，你猜函數在每一點的高度（慣例上取 0，不失一般性，因為可以先把資料減去均值）；\(k(\cdot,\cdot)\) 是核函數（kernel）——輸入任兩點，回傳這兩點函數值的共變異數，「兩點靠多近，函數值就得多同步」。定義的意思是：不管你挑哪 n 個輸入點、n 多大，這些點上的函數值永遠服從一個以核矩陣 \(K_{X,X}\) 為共變異的多元高斯分布。一個 GP 就是「無限維的高斯分布」——你永遠只需要跟它的有限維切片打交道，而每個切片都是熟悉的多元高斯。

這回應了第 1 章埋的伏筆：參數式模型的參數個數固定，非參數式模型的複雜度隨資料成長——GP 是後者的代表。它不學一組固定的權重，而是學「給定觀測資料後，整個函數空間上的後驗分布」；資料越多，後驗掐得越緊，模型的表達也越細。

比喻： GP 先驗像一條攤在桌上的彈性長繩。核函數規定繩子的「質地」：長度尺度大，繩子硬挺，起伏舒緩；長度尺度小，繩子柔軟，皺褶密集。從先驗抽樣，就是隨手撥動這條繩子看它躺成什麼形狀——每次都不同，但質地不變。觀測資料則是把繩子在幾個位置用圖釘釘住：釘點之間繩子只能小幅擺動（後驗變異數小），離釘點遠的繩段依然可以大幅甩動（後驗變異數大）。GP 的後驗，就是「所有通過釘點附近的繩形」構成的分布。

圖 03-1 的右半就是這個機制的實跑結果：六個觀測點（圖釘）把後驗樣本壓成幾乎同一條曲線，±2 個標準差的信賴帶在資料密集處寬約 ±0.16，到了沒有資料的右端張開到約 ±1.08——差了將近七倍。模型不需要任何外掛，就把「哪裡可信、哪裡不可信」畫在臉上。

五、核函數：一顆調整「平滑度」的旋鈕

核函數是 GP 的靈魂：它把你對函數的所有先驗直覺——多平滑、擺動多大、有沒有週期性——編碼成一個二元函數。最常用的是 RBF 核（radial basis function，又叫平方指數核）：

逐項拆解：\(\|x-x'\|\) 是兩個輸入點的距離——RBF 核只看距離，不看絕對位置，這叫平穩性；\(\ell\) 是長度尺度（length-scale），白話講就是曲線的記性有多長——一筆觀測能替多大範圍的鄰居作證：兩點相距遠小於 \(\ell\) 時共變異接近滿格、函數值幾乎連動，相距一旦遠超過 \(\ell\)，共變異趨近 0，這筆觀測對那個位置就再也說不上話；\(\sigma_f^2\) 是訊號變異數，控制函數擺動的整體幅度——把它開大，先驗允許的曲線上下振幅就大。RBF 核產生無限次可微的極平滑函數；如果嫌它太平滑，Matern 核家族提供一顆額外的旋鈕來調函數的粗糙度，金融時間序列建模常用它。數學上核不能亂選：它必須是對稱半正定函數，這樣任何點集上的核矩陣才是合法的共變異矩陣。

長度尺度是實務上最要緊的超參數，它直接對應第 1 章的偏差–方差取捨：\(\ell\) 太小，模型柔軟到把雜訊當訊號，過擬合；\(\ell\) 太大，模型硬挺到抹平真實結構，欠擬合。圖 03-2 用同一批資料、三個 \(\ell\) 實跑對照——數字出自作業三，先看圖建立直覺，第七節再講怎麼讓資料自己挑 \(\ell\)。

六、GP 迴歸實戰：後驗均值與變異數的直覺

把觀測雜訊放進來：\(y_i = f(x_i) + \epsilon_i\)，\(\epsilon_i \sim \mathcal{N}(0, \sigma_n^2)\)。因為 GP 的任何有限維切片都是高斯，「觀測值向量」與「測試點上的函數值」合起來也是一個大高斯，用多元高斯的條件分布公式（機制同第一節的權重後驗）直接得到 GP 迴歸的預測分布——均值函數取 0 的慣例下：

逐項拆解：\(K_{X,X}\) 是觀測點之間的核矩陣，\(K_{*,X}\) 是測試點對觀測點的核向量，\(K_{*,*}\) 是測試點對自己的核值（＝先驗變異數 \(\sigma_f^2\)）。後驗均值可以讀成「觀測值的加權平均」：把 \(\alpha := (K_{X,X}+\sigma_n^2 I)^{-1}\mathbf{y}\) 先算好，預測就是 \(K_{*,X}\,\alpha\)——離測試點近的觀測值（核值大）拿高權重，遠的拿低權重。後驗變異數讀成「先驗不確定性，減去被資料解釋掉的部分」：第二項越大，代表測試點與觀測資料的關聯越強、砍掉的不確定性越多；測試點離所有資料都遠時，第二項趨近 0，變異數回升到先驗水準——模型攤手說「這裡我只剩下先驗可講」。注意一件事：變異數公式裡完全沒有 \(\mathbf{y}\)——不確定性只由「觀測點的位置」決定，跟觀測到什麼值無關，這是高斯假設的特色。

上代碼。場景：選擇權的隱含波動率微笑——同一到期日、不同履約價的隱含波動率連成的曲線。市場只在幾個履約價上有可靠報價，其他位置要靠模型補：

import numpy as np

rng = np.random.default_rng(11)

# 情境：選擇權隱含波動率微笑。x=價內外程度（moneyness），y=隱含波動率
def true_smile(m):
    return 0.20 + 0.35 * (m - 1.0) ** 2      # 真實微笑曲線（未知，僅用來造資料）

m_obs = np.array([0.80, 0.87, 0.93, 1.00, 1.06, 1.12, 1.20])   # 有報價的履約價
y_obs = true_smile(m_obs) + rng.normal(0.0, 0.004, m_obs.size) # 報價含雜訊

ell, sig_f, sig_n = 0.18, 0.10, 0.004        # 長度尺度、訊號幅度、雜訊標準差

def rbf(a, b):
    d = a[:, None] - b[None, :]
    return sig_f**2 * np.exp(-0.5 * (d / ell) ** 2)

K = rbf(m_obs, m_obs) + sig_n**2 * np.eye(m_obs.size)
alpha = np.linalg.solve(K, y_obs)

def gp_predict(m_new):
    Ks = rbf(m_new, m_obs)                    # 測試點與觀測點的核
    mu = Ks @ alpha                           # 後驗均值
    cov = rbf(m_new, m_new) - Ks @ np.linalg.solve(K, Ks.T)
    return mu, np.sqrt(np.clip(np.diag(cov), 0, None))

for m in (1.00, 1.16, 1.35, 1.60):
    mu, sd = gp_predict(np.array([m]))
    print(f"moneyness {m:.2f}｜GP 預測波動率 {mu[0]:.4f} ± {2*sd[0]:.4f}"
          f"｜真實值 {true_smile(m):.4f}")
# 輸出：
# moneyness 1.00｜GP 預測波動率 0.1986 ± 0.0059｜真實值 0.2000
# moneyness 1.16｜GP 預測波動率 0.2106 ± 0.0067｜真實值 0.2090
# moneyness 1.35｜GP 預測波動率 0.1827 ± 0.0869｜真實值 0.2429
# moneyness 1.60｜GP 預測波動率 0.0319 ± 0.1960｜真實值 0.3260

四行輸出是一堂完整的風險課。在報價覆蓋區內（1.00、1.16），預測貼著真實值，帶寬只有 ±0.006 上下。往外走到 1.35，帶寬爆增十幾倍到 ±0.087——真實值 0.2429 仍落在帶內，模型用帶寬保住了誠實。但走到 1.60，事情變了：預測均值一路滑向先驗均值 0（GP 在資料外會回歸先驗），連 ±0.196 的帶都罩不住真實值 0.326。教訓有兩層：第一層，帶寬爆寬本身就是「別用我」的訊號，任何下游系統都該設帶寬上限；第二層更深——不確定性量化本身也是模型的產物，先驗設錯（波動率曲面均值哪會是 0），誤差條在外推區一樣失真。務實解法包括把資料標準化、給均值函數合理的形狀，或改用能外推的核組合（例如線性核加 RBF 核，讓價外深處保持線性），原書對這一手有專門的實驗。

七、超參數與計算帳單：奧卡姆剃刀，以及 \(O(N^3)\)

長度尺度 \(\ell\)、訊號變異數 \(\sigma_f^2\)、雜訊變異數 \(\sigma_n^2\) 這些超參數怎麼定？GP 給了一條非常貝葉斯的路：最大化邊際似然（marginal likelihood，也叫證據 evidence）——「在這組超參數之下，模型生出這批觀測資料的機率」：

逐項拆解：\(\lambda\) 代表全部超參數（別跟 ridge 的懲罰係數混淆，這裡沿用慣例符號）；第一項是擬合項——資料在模型眼中的「意外程度」，模型解釋得越好，這項越大（負得越少）；第二項 \(-\tfrac{1}{2}\log\det K_y\) 是複雜度罰項——核矩陣的行列式衡量「這組超參數允許的函數花樣有多少」，允許的花樣越多，行列式越大，罰得越重；第三項是常數。擬合與複雜度在同一條式子裡自動討價還價，文獻稱之為自動奧卡姆剃刀：不必切驗證集，先驗自己會懲罰過度複雜的設定。

比喻： 邊際似然像賭場結算。每組超參數是一位賭客，它的先驗決定籌碼怎麼下注在「所有可能出現的資料集」上。長度尺度極小的賭客什麼曲線都覺得可能，籌碼攤得極薄，就算實際資料開出，它押在上面的籌碼也沒多少；長度尺度太大的賭客把籌碼全押在少數幾種平滑曲線上，資料一旦不配合就血本無歸；押注集中程度剛好、又真的押中的賭客，結算金額最高。邊際似然結算的，就是每位賭客押在「實際觀測到的這批資料」上的籌碼。

作業三會實跑這件事：三個長度尺度裡，邊際似然最高的那個（\(\ell=0.7\)，對數邊際似然 \(-3.76\)）恰好就是樣本外 MSE 最低的那個；把 \(\ell\) 放上格點掃描，邊際似然挑中 \(\ell \approx 0.60\)，測試 MSE 0.0180，跟事後諸葛選出的最優值幾乎一樣——而整個過程沒動用任何測試資料。

再說帳單。訓練 GP 要解 \(K_y^{-1}\mathbf{y}\) 並算 \(\log\det K_y\)，標準做法是對 \(N \times N\) 矩陣做 Cholesky 分解，成本 \(O(N^3)\)；預測時每個測試點的均值只要 \(O(N)\)、變異數 \(O(N^2)\)。\(N\) 上千還算輕鬆，上十萬就是災難——這是 GP 最著名的軟肋。解法方向是可擴展 GP：用 \(m \ll N\) 個誘導點架橋，配合結構化核插值（SKI）與 FFT 把訓練成本壓到近似線性、單點預測壓到 \(O(1)\)，原書以 MSGP 框架為例給了完整的實驗；概念上你只要記得「\(O(N^3)\) 有藥醫，代價是近似」。

八、金融應用：把昂貴定價器做成會標示信心的快取

進正題：GP 在金融裡最成熟的應用之一，是替昂貴的定價器建代理模型（surrogate model）——一個學會「模仿既有模型輸出」的模型。

先看痛點。奇異選擇權、可贖回債、結構型商品的定價常常沒有閉式解，要靠蒙地卡羅模擬或偏微分方程數值解，單次定價動輒數十毫秒到數秒。平常這不是問題，但三個場景會把它變成問題：一是風險系統要對整本投資組合、上千種市場情境重算價值（風險值、壓力測試、XVA 的巢狀模擬），呼叫次數是「部位數 × 情境數 × 時點數」的乘積；二是造市報價要在市場每一跳後刷新幾百檔商品的報價與避險參數——台灣的權證自營商對此再熟悉不過，上千檔權證掛牌、盤中每檔都要維持買賣報價與 delta 避險，定價引擎的每一毫秒都是錢；三是模型校準要在最佳化迴圈裡反覆定價。三個場景的共同結構是：同一個定價函數，被在略有不同的輸入上呼叫成千上萬次。

這正是快取的形狀，而 GP 是特別聰明的快取。工作流分兩段。離線：挑一組代表性的輸入網格（現貨價 × 波動率 × 利率……），呼叫昂貴定價器把這些點的價格算出來——這是唯一要付大錢的地方；然後用這些「輸入–價格」對訓練一個 GP，訓練成本 \(O(N^3)\) 一次付清。線上：任何新輸入進來，GP 後驗均值直接吐出價格，單次查詢是毫秒級以下的矩陣運算；同時後驗變異數免費附上「這個價格多可靠」。

比起查表插值，GP 快取多給兩樣東西。第一是解析的 Greeks。後驗均值對輸入的偏導數有閉式：

逐項拆解：\(\alpha\) 在訓練（定價）階段就算好了，是固定向量；要 delta，只需把「核函數對測試輸入的導數」（RBF 核有現成的解析式）跟 \(\alpha\) 做一次內積——不必對定價器做有限差分（省掉兩次額外定價與差分誤差），也不必訓練時見過任何 Greeks 資料。第二是信心標示：查詢點落在訓練網格覆蓋不到的區域時，變異數自動升高——上一節微笑曲線的教訓在這裡變成制度設計：帶寬超標，就別用快取，回頭呼叫原定價器補算，順手把新點加進訓練集。

實跑一個最小可行版本。用 Black–Scholes 閉式解扮演「昂貴定價器」（實務上換成 Heston、蒙地卡羅或 PDE 求解器，流程一字不改），離線只呼叫它 12 次：

import numpy as np
from math import erf, exp, log, sqrt

# 「昂貴定價器」：Black-Scholes 買權（實務上換成 Heston/MC/PDE 都一樣）
K_strike, r, vol, tau = 50.0, 0.01, 0.25, 0.5

def Phi(x):
    return 0.5 * (1.0 + erf(x / sqrt(2.0)))

def bs_call(S):
    d1 = (log(S / K_strike) + (r + 0.5 * vol**2) * tau) / (vol * sqrt(tau))
    d2 = d1 - vol * sqrt(tau)
    return S * Phi(d1) - K_strike * exp(-r * tau) * Phi(d2)

def bs_delta(S):
    d1 = (log(S / K_strike) + (r + 0.5 * vol**2) * tau) / (vol * sqrt(tau))
    return Phi(d1)

# 離線階段：只呼叫昂貴定價器 12 次，建訓練網格
S_grid = np.linspace(30.0, 70.0, 12)
C_grid = np.array([bs_call(s) for s in S_grid])

ell, sig_f = 9.0, 12.0                       # 核超參數（價格曲面很平滑，ell 取大）
def rbf(a, b):
    d = a[:, None] - b[None, :]
    return sig_f**2 * np.exp(-0.5 * (d / ell) ** 2)

K = rbf(S_grid, S_grid) + 1e-8 * np.eye(S_grid.size)   # 無雜訊：只加數值抖動
alpha = np.linalg.solve(K, C_grid)

def gp_price(S_new):
    return rbf(S_new, S_grid) @ alpha

def gp_delta(S_new):
    """解析微分：d mu/dS* = d k/dS* @ alpha，不必有限差分"""
    d = S_new[:, None] - S_grid[None, :]
    dK = -rbf(S_new, S_grid) * d / ell**2
    return dK @ alpha

# 線上階段：查詢任意現貨價（都不在訓練網格上）
S_test = np.array([37.3, 46.1, 50.0, 58.8, 66.4])
p_err = np.abs(gp_price(S_test) - np.array([bs_call(s) for s in S_test]))
d_err = np.abs(gp_delta(S_test) - np.array([bs_delta(s) for s in S_test]))
print("測試點：", S_test)
print("定價最大絕對誤差：", f"{p_err.max():.5f}", "美元")
print("Delta 最大絕對誤差：", f"{d_err.max():.5f}")
print(f"對照：S=58.8 時 GP delta {gp_delta(np.array([58.8]))[0]:.4f}"
      f"，BS 解析 delta {bs_delta(58.8):.4f}")
# 輸出：
# 測試點： [37.3 46.1 50.  58.8 66.4]
# 定價最大絕對誤差： 0.00013 美元
# Delta 最大絕對誤差： 0.00249
# 對照：S=58.8 時 GP delta 0.8496，BS 解析 delta 0.8494

十二個訓練點，換來全域定價誤差不超過 0.00013 美元；更妙的是 delta——GP 從頭到尾沒看過任何一筆 delta 資料，純靠對核函數解析微分，就把 delta 追到小數第三位（0.8496 對 0.8494）。要 gamma 就再微分一次，要 vega 就把波動率也放進輸入維度。這裡藏著一個值得誠實面對的哲學問題：代理模型學的是既有模型的輸出，不是市場本身——它繼承原模型全部的對與錯，只負責「更快」，不負責「更對」。這與第 1 章「機器學習從資料找規律」的敘事方向相反，卻是機器學習在計算金融裡吃下的一塊紮實地盤。

最後兩句話交代延伸方向。其一：投資組合層級的應用（例如信用估值調整 XVA）要同時預測一整籃商品的價格，這時單輸出 GP 各管各的會漏掉商品間的相關結構，多輸出 GP 用「輸入核 × 任務共變異矩陣」的張量結構把它補上，原書第六節有完整公式。其二：盤中定價模型重校準後代理要跟著更新，貝葉斯的「舊後驗當新先驗」讓 GP 能線上增量學習，不必每次重訓。

九、收官：貝葉斯之眼，下一站神經網路

本章走了一條完整的弧線：先驗→後驗搬上迴歸權重（點估計變成分布）、MAP 揭開 ridge 的身世（先驗＝正則化）、跳進函數空間（GP＝直接對函數放先驗）、核函數與邊際似然（把「模型該多複雜」交給資料裁決）、最後落地在代理定價（不確定性量化變成可營運的制度）。留一個對照給下一章：神經網路走的是另一條路——參數極多、點估計為主、正則化以顯式懲罰項登場。當你在第 4 章看到權重衰減時，記得回頭想這一章的彈簧；當你看到 dropout 時，想一想「對權重的分布取平均」——貝葉斯之眼會讓那些看似工程直覺的技巧，露出統計的骨架。

§03原書對照

本課以白話重組原書第三章的主幹，以下內容原書有、但本課未展開，按頁碼供進階讀者深入。其一，無雜訊情形的完整推演：從兩參數線性模型出發，證明先驗誘導的函數值共變異為 \(1 + x_i x_j\)，作為「權重先驗誘導函數分布」的最小示範（pp.82–83，式 3.1–3.6）。其二，權重後驗的嚴格推導：以自然參數化配方（completing the square）處理高斯密度連乘，得後驗矩的閉式（pp.84–85，式 3.12–3.21），並附轉載自 Bishop 的四格圖，展示後驗隨觀測數 0→1→2→20 逐步收縮的過程（p.86，圖 3.2）。其三，頻率學派對照組：最大概似估計的梯度推導、雜訊變異數的估計式，與 OLS 作為正交投影的幾何解讀（pp.86–88）。其四，Schur 恆等式：以分割共變異矩陣的反矩陣性質另闢一條推導預測分布的路，並給出區塊反矩陣的完整代數（pp.89–90，式 3.28–3.40）。其五，GP 與計量經濟學的銜接：AR(p) 過程是配特定 Matern 類核之 GP 的離散時間等價物，故 GP 可視為古典計量技術的貝葉斯非參數推廣；並回顧 kriging 文獻、期望短缺的 metamodeling（以 GP 免去巢狀模擬的內層重算）、CVA 計算等金融應用，及「多輸出 GP 才能刻畫組合層級預測誤差」的論證（pp.92–93）。其六，Matern 核的完整定義與參數意義（p.94，式 3.42），邊際似然對超參數的解析梯度（p.95，式 3.46–3.49）。其七，線上學習：定價模型盤中重校準時，把舊後驗當新先驗的增量更新式（p.96，式 3.50）。其八，可擴展 GP 全節：SKI 近似核、誘導點網格上 Kronecker 與 Toeplitz 結構如何配合 FFT 加速特徵分解，及與 SoR、FITC 等經典誘導點法的關係（pp.96–98，式 3.51–3.52）。其九，Heston 模型定價實驗：隨機波動率動態與參數表、買權賣權同時擬合、30×30 訓練網格對 40×40 測試網格的價格曲面重建與縮放到單位區間的數值穩定性註記（pp.98–100，圖 3.3）、線性核加 SE 核的組合如何讓深度價內外的外推保持線性（p.100，圖 3.4）、GP delta 與 vega 對解析值的重合（p.101，圖 3.5）、mesh-free 取樣的收斂與 SKI 訓練時間實驗（pp.101–103，圖 3.6–3.8）。其十，多輸出 GP 全節：矩陣變量高斯、輸入核與任務共變異的 Kronecker 結構、預測分布與負對數邊際似然（pp.103–105，式 3.57–3.62）。其十一，章末習題四題（權重後驗推導、常態共軛、GP 預測分布推導、GP 學習買權價格形狀的程式題）與三組選擇題詳解，及 CVA 建模 notebooks 的說明（pp.106–108）。原書第三章對應印刷頁 pp.81–109。

§04作業和解答

import numpy as np

rng = np.random.default_rng(21)

def f_true(x):
    return 0.03 * np.sin(2.0 * np.pi * x) + 0.01 * x   # 某月內季節性訊號

x_obs = rng.uniform(0.0, 1.0, 10)
y_obs = f_true(x_obs) + rng.normal(0.0, 0.008, 10)

def rbf(a, b, ell=0.25, sig_f=0.03):
    d = a[:, None] - b[None, :]
    return sig_f**2 * np.exp(-0.5 * (d / ell) ** 2)

def gp_fit_predict(sig_n, x_new):
    K = rbf(x_obs, x_obs) + max(sig_n**2, 1e-10) * np.eye(x_obs.size)
    L = np.linalg.cholesky(K)                       # Cholesky 分解，穩定求解
    alpha = np.linalg.solve(L.T, np.linalg.solve(L, y_obs))
    Ks = rbf(x_new, x_obs)
    mu = Ks @ alpha
    v = np.linalg.solve(L, Ks.T)
    var = np.clip(np.diag(rbf(x_new, x_new)) - np.sum(v**2, axis=0), 0, None)
    return mu, np.sqrt(var)

for sig_n in (0.008, 0.0):
    mu_at_obs, _ = gp_fit_predict(sig_n, x_obs)
    mse_train = np.mean((mu_at_obs - y_obs) ** 2)
    mu_far, sd_far = gp_fit_predict(sig_n, np.array([1.8]))
    mu_in, sd_in = gp_fit_predict(sig_n, np.array([0.5]))
    print(f"sigma_n={sig_n:.3f}｜訓練 MSE {mse_train:.2e}"
          f"｜x=0.5 帶寬 ±{2*sd_in[0]:.4f}｜x=1.8 帶寬 ±{2*sd_far[0]:.4f}")
# 輸出：
# sigma_n=0.008｜訓練 MSE 7.70e-05｜x=0.5 帶寬 ±0.0117｜x=1.8 帶寬 ±0.0600
# sigma_n=0.000｜訓練 MSE 8.55e-06｜x=0.5 帶寬 ±0.0002｜x=1.8 帶寬 ±0.0599

import numpy as np

rng = np.random.default_rng(8)

n, p = 40, 3
X = rng.normal(size=(n, p))
theta_true = np.array([0.6, -0.4, 0.2])
sigma_n = 0.5
y = X @ theta_true + rng.normal(0.0, sigma_n, n)

sigma_theta = 0.8                            # 等向高斯先驗 N(0, sigma_theta^2 I)
lam = sigma_n**2 / sigma_theta**2            # 對應的 ridge 懲罰強度

# 路線一：ridge 閉式解 (X^T X + lam I)^{-1} X^T y
theta_ridge = np.linalg.solve(X.T @ X + lam * np.eye(p), X.T @ y)

# 路線二：貝葉斯後驗均值 Sigma'(X^T y / sigma_n^2)，先驗均值為 0
Sigma_post = np.linalg.inv(np.eye(p) / sigma_theta**2 + (X.T @ X) / sigma_n**2)
theta_map = Sigma_post @ (X.T @ y) / sigma_n**2

print("ridge 解：", np.round(theta_ridge, 6))
print("MAP  解：", np.round(theta_map, 6))
print("最大差異：", np.abs(theta_ridge - theta_map).max())

# 加碼驗證：lam -> 0（先驗無限寬）時收斂到 OLS
theta_ols = np.linalg.lstsq(X, y, rcond=None)[0]
Sigma_wide = np.linalg.inv(np.eye(p) / 1e8**2 + (X.T @ X) / sigma_n**2)
theta_wide = Sigma_wide @ (X.T @ y) / sigma_n**2
print("先驗極寬 vs OLS 最大差異：", np.abs(theta_wide - theta_ols).max())
# 輸出：
# ridge 解： [ 0.556329 -0.405685  0.230894]
# MAP  解： [ 0.556329 -0.405685  0.230894]
# 最大差異： 2.7755575615628914e-17
# 先驗極寬 vs OLS 最大差異： 2.7755575615628914e-16

import numpy as np

rng = np.random.default_rng(33)

def f_true(x):
    return np.sin(3.0 * x)

x_tr = np.sort(rng.uniform(0.0, 3.0, 12))
y_tr = f_true(x_tr) + rng.normal(0.0, 0.15, 12)
x_te = np.linspace(0.0, 3.0, 200)
y_te = f_true(x_te)

sig_n = 0.15

def rbf(a, b, ell, sig_f=1.0):
    d = a[:, None] - b[None, :]
    return sig_f**2 * np.exp(-0.5 * (d / ell) ** 2)

def fit_eval(ell):
    K = rbf(x_tr, x_tr, ell) + sig_n**2 * np.eye(x_tr.size)
    L = np.linalg.cholesky(K)
    alpha = np.linalg.solve(L.T, np.linalg.solve(L, y_tr))
    mu = rbf(x_te, x_tr, ell) @ alpha
    mse = np.mean((mu - y_te) ** 2)
    # 對數邊際似然：擬合項 - 複雜度罰項 - 常數
    logml = (-0.5 * y_tr @ alpha - np.sum(np.log(np.diag(L)))
             - 0.5 * x_tr.size * np.log(2 * np.pi))
    return mse, logml

for ell in (0.08, 0.7, 3.0):
    mse, logml = fit_eval(ell)
    print(f"ell={ell:4.2f}｜測試 MSE {mse:.4f}｜對數邊際似然 {logml:7.2f}")

# 讓邊際似然自己挑：掃一遍格點
grid = np.linspace(0.05, 3.0, 60)
scores = [fit_eval(l)[1] for l in grid]
best = grid[int(np.argmax(scores))]
print(f"邊際似然最大處 ell ≈ {best:.2f}，其測試 MSE {fit_eval(best)[0]:.4f}")
# 輸出：
# ell=0.08｜測試 MSE 0.1769｜對數邊際似然  -11.44
# ell=0.70｜測試 MSE 0.0206｜對數邊際似然   -3.76
# ell=3.00｜測試 MSE 0.3454｜對數邊際似然  -95.57
# 邊際似然最大處 ell ≈ 0.60，其測試 MSE 0.0180

§05參考資料

• Rasmussen, C. E., & Williams, C. K. I. (2006). Gaussian Processes for Machine Learning. MIT Press（全書免費線上版） — GP 的聖經級教科書，本章預測公式、邊際似然與核函數理論的權威出處。
• Görtler, J., Kehlbeck, R., & Deussen, O. (2019). A Visual Exploration of Gaussian Processes. Distill — 互動式視覺化講解 GP 的先驗、後驗與核函數，建立幾何直覺的最佳入口。
• Wilson, A. G., & Nickisch, H. (2015). Kernel Interpolation for Scalable Structured Gaussian Processes (KISS-GP). arXiv:1503.01057 — 結構化核插值（SKI）的原始論文，第七節可擴展 GP 的方法源頭。
• Gardner, J., Pleiss, G., Wu, R., Weinberger, K., & Wilson, A. (2018). Product Kernel Interpolation for Scalable Gaussian Processes. arXiv:1802.08903 — 把 SKI 推進到高維輸入的 MSGP 框架，原書可擴展 GP 一節的直接參照。
• Crépey, S., & Dixon, M. F. (2019). Gaussian Process Regression for Derivative Portfolio Modeling and Application to CVA Computations. arXiv:1901.11081 — 多輸出 GP 用於衍生品組合與 CVA 計算，本章代理定價應用的進階延伸（作者含原書作者）。
• scikit-learn：Gaussian Processes 官方文件 — GP 迴歸與核函數的工程實作參考，可對照本章的 numpy 手寫版。
• Murphy, K. P. (2012). Machine Learning: A Probabilistic Perspective. MIT Press（官方頁面） — 貝葉斯線性迴歸與邊際似然（§15.2.4）的嚴格處理，原書 evidence 公式的直接引用來源。
• Dixon, M. F., Halperin, I., & Bilokon, P. (2020). Machine Learning in Finance: From Theory to Practice. Springer — 本課程對照之原書，第三章 pp.81–109。