10.3 數學演算法的實戰案例#

import numpy as np
from scipy.optimize import minimize
import random

rng = np.random.default_rng(7)

# ========== 案例一：房價回歸完整管線 ==========
n = 100
area = rng.uniform(15, 60, n)                      # 坪數
age = rng.uniform(0, 40, n)                        # 屋齡
price = 3.5*area - 0.8*age + 50 + rng.normal(0, 8, n)

A = np.column_stack([np.ones(n), area, age])
Q, R = np.linalg.qr(A)
coef = np.linalg.solve(R, Q.T @ price)
resid = price - A @ coef
se = np.sqrt(resid @ resid / (n - 3)) * np.sqrt(np.diag(np.linalg.inv(A.T @ A)))
print("=== 案例一：回歸 ===")
for name, c, s in zip(["截距", "坪數係數", "屋齡係數"], coef, se):
    lo, hi = c - 1.96*s, c + 1.96*s
    print(f"{name}: {c:.3f} ± 1.96·{s:.3f} → 95% CI [{lo:.3f}, {hi:.3f}]")

# ========== 案例二：Markowitz 投資組合 ==========
# 3 檔資產的期望報酬與共變異矩陣
mu = np.array([0.08, 0.12, 0.15])                  # 期望年報酬
Sigma = np.array([[0.010, 0.002, 0.001],
                  [0.002, 0.020, 0.004],
                  [0.001, 0.004, 0.035]])
target = 0.11                                      # 目標報酬 11%

def port_var(w):
    return w @ Sigma @ w

cons = [{"type": "eq", "fun": lambda w: w.sum() - 1},           # 滿倉
        {"type": "ineq", "fun": lambda w: w @ mu - target}]     # 報酬達標
res = minimize(port_var, np.ones(3)/3, constraints=cons,
               bounds=[(0, 1)]*3, method="SLSQP")
w = res.x
print("\n=== 案例二：投資組合 ===")
print(f"權重 = {np.round(w, 4)}, 年化波動 = {np.sqrt(port_var(w)):.4f}")
print(f"實際報酬 = {w @ mu:.4f}（目標 {target}）")

# 蒙地卡羅驗證：隨機權重找出的最小波動不應更好
best_var, best_w = np.inf, None
for _ in range(200_000):
    p = rng.dirichlet(np.ones(3))
    if p @ mu >= target and p @ Sigma @ p < best_var:
        best_var, best_w = p @ Sigma @ p, p
print(f"蒙地卡羅最佳隨機組合波動² = {best_var:.6f}（應接近且不顯著低於 SLSQP 的 {port_var(w):.6f}）")

# ========== 案例三：logistic 分類 ==========
m = 300
x1 = rng.normal(0, 1, m)
x2 = rng.normal(0, 1, m)
y = np.where(2*x1 - x2 + 0.5 + rng.normal(0, 1, m) > 0, 1, -1)
X = np.column_stack([np.ones(m), x1, x2])

def logistic_loss(w):
    z = y * (X @ w)
    return np.logaddexp(0, -z).mean()

def logistic_grad(w):
    z = y * (X @ w)
    sig = 1 / (1 + np.exp(-z))                     # σ(z)
    return -((1 - sig) * y) @ X / m

# 梯度檢查（10.2 節的必修儀式）
w0 = rng.normal(size=3)
eps = 1e-6
num = np.array([(logistic_loss(w0 + eps*np.eye(3)[i])
                 - logistic_loss(w0 - eps*np.eye(3)[i])) / (2*eps)
                for i in range(3)])
print("\n=== 案例三：logistic ===")
print("梯度檢查誤差 =", np.abs(num - logistic_grad(w0)).max())

w, eta = np.zeros(3), 0.5
for it in range(2000):
    g = logistic_grad(w)
    w -= eta * g
    if np.linalg.norm(g) < 1e-10:
        break
acc = np.mean(np.sign(X @ w) == y)
print(f"訓練 {it+1} 次迭代收斂，損失 {logistic_loss(w):.4f}，準確率 {acc:.2%}")
try:
    ref = minimize(logistic_loss, np.zeros(3), jac=logistic_grad,
                   method="BFGS")
    print(f"BFGS 對照損失 = {ref.fun:.4f}（應一致）")
except Exception as e:
    print("（BFGS 對照失敗：", e, "）")