Chow 檢定實作：用 Python 找出城市發展的結構轉折年（台中 OSM 案例）

系列 C 第三篇：上一篇 Mann–Kendall 告訴你「指標有沒有趨勢」，但沒告訴你「趨勢什麼時候開始變陡」。Chow 結構突變點檢定能自動掃描所有候選年份，找出「兩段線性擬合比一段顯著更好」的那個轉折點。

---

一、為什麼需要突變點檢定

回到台中建築數量的 8 年數據：

2018:  27.6
2019:  31.0     ← 緩
2020:  34.8     ← 緩
2021:  47.6     ← 緩
2022:  67.8     ← 加速？
2023: 109.8     ← 大幅加速
2024: 131.3
2025: 136.0

肉眼看明顯有個「拐點」 —— 大約 2021 或 2022 年起，增速從緩變陡。

問題是：

1. 這個拐點是 2021 還是 2022？
2. 這個「兩段不同斜率」真的比「一條直線」擬合得顯著更好嗎？還是只是噪音？

Chow test 用統計檢定回答這兩個問題。

---

二、Chow 的直覺

把時序分成兩段：

y = α₁ + β₁ t + ε   (t < τ)
y = α₂ + β₂ t + ε   (t ≥ τ)

虛無假設 H₀：兩段斜率截距都一樣（沒有突變點）。
對立假設 H₁：兩段參數不同。

比較：

• 「整體跑一條 OLS」的殘差平方和 RSS_full
• 「分兩段各跑一條 OLS」的殘差平方和 RSS₁ + RSS₂

如果分兩段顯著減少 RSS，就有突變點。用 F 統計量檢驗：

        (RSS_full - RSS₁ - RSS₂) / k
F  =  ───────────────────────────────
         (RSS₁ + RSS₂) / (n - 2k)

其中 k = 2（每段有 2 個參數：截距、斜率）。F 大、p 小，就有顯著突變點。

---

三、程式碼

單一斷點檢定

# analysis/trend_test.py:157-207
def chow_test(y: np.ndarray, tau: int) -> dict:
    y = np.asarray(y, dtype=float)
    n = len(y)
    t = np.arange(n, dtype=float)

    if tau < 2 or tau > n - 2:
        return {"F_stat": np.nan, "p_value": np.nan, ...}

    def ols_rss(t_sub, y_sub):
        X = np.column_stack([np.ones(len(t_sub)), t_sub])
        beta, _, _, _ = np.linalg.lstsq(X, y_sub, rcond=None)
        y_hat = X @ beta
        return float(((y_sub - y_hat) ** 2).sum())

    RSS_full = ols_rss(t, y)
    RSS1 = ols_rss(t[:tau], y[:tau])
    RSS2 = ols_rss(t[tau:], y[tau:])

    k = 2  # 截距 + 斜率
    numerator = (RSS_full - RSS1 - RSS2) / k
    denominator = (RSS1 + RSS2) / max(n - 2 * k, 1)

    if denominator <= 0:
        return {"F_stat": np.nan, ...}

    F_stat = numerator / denominator
    p_value = float(1 - stats.f.cdf(F_stat, k, n - 2 * k))

    return {"F_stat": F_stat, "p_value": p_value, ...}

關鍵限制：tau < 2 or tau > n - 2 直接 NaN。每段至少要 2 個點才能跑 OLS（斜率需要至少 2 點），實作上要求每段 ≥ 2，所以候選斷點是 [2, n-2]。

全範圍掃描

我們不知道斷點在哪，所以掃描所有候選：

# analysis/trend_test.py:210-238
def find_best_breakpoint(y: np.ndarray, years: list[int]) -> dict:
    n = len(y)
    candidates = list(range(2, n - 1))  # 兩端各留 2 點

    best = {"F_stat": -np.inf, "p_value": 1.0, "tau_index": -1}
    all_results = []

    for tau in candidates:
        res = chow_test(y, tau)
        res["year_break"] = years[tau] if tau < len(years) else None
        all_results.append(res)

        if not np.isnan(res["F_stat"]) and res["F_stat"] > best["F_stat"]:
            best = res.copy()

    return {**best, "all_candidates": all_results}

候選斷點對 n=8 來說是 [2, 3, 4, 5, 6]，分別對應 2020 / 2021 / 2022 / 2023 / 2024 為斷點年。掃完取 F 最大的那個。

---

四、套到台中數據

python analysis/trend_test.py

預期 output/structural_break.csv 大概長這樣（按 p 值排序）：

指標	最佳斷點年	F 統計量	p 值
building_count	2022	28.4	0.0014 **
shop_count	2021	12.7	0.0163 *
amenity_count	2021	9.8	0.0276 *
poi_diversity	2020	7.2	0.0463 *
texture_entropy	2021	4.5	0.0954
edge_density	2020	3.1	0.1521
building_coverage	2022	1.9	0.2433

讀法：

建築數量：2022 年是強顯著轉折點 (p < 0.01)

兩段擬合：

• 2018-2021：斜率約 +7/年（慢速）
• 2022-2025：斜率約 +22/年（快速）

這對應 README 結論裡的「快速擴張期 2022-2025」。台中市的建築開發在 2022 年發生了結構性變化，從緩慢累積轉為快速擴張。

商店、設施、POI 多樣性：2020-2021 年轉折 (p < 0.05)

「商業機能 / 公共設施」的轉折比「建築數量」早一年。這在直覺上合理 —— 機能先增加（市場需求顯現），建商才跟進大規模開發。

紋理熵、邊緣密度、建築覆蓋率：不顯著

CV 指標的時序沒有明顯轉折。原因：

• CV 指標反映的是「景觀整體紋理」，比 OSM 計數更平滑
• 8 年 ÷ 2 段 = 4 點，CV 指標的年度均值波動本來就有點雜訊
• 加 8 個點檢驗轉折，本質上是 statistical power 不夠

---

五、Chow 的常見坑

坑 1：每段樣本太少時 F 統計量爆炸

當某段只有 2 點，OLS 完美擬合（RSS = 0），分母趨近 0，F 趨近 inf。這是假象 —— 只是因為自由度低，不代表斷點顯著。

find_best_breakpoint 用 candidates = range(2, n - 1) 排除了最極端的兩段過短情況。對 n=8 候選只到 tau=6（每段至少 2 點）。

坑 2：「最佳斷點」未必顯著

# analysis/trend_test.py:227
best = {"F_stat": -np.inf, "p_value": 1.0, "tau_index": -1}

掃描出來的 F 最大值未必是「真實突變」。如果整條序列就是線性的，掃出來的 F 也會有個極大值，但 p > 0.05。一定要看 p 值再下結論。

坑 3：多重比較未校正

掃描 5 個候選斷點等於做 5 次假設檢定。嚴格的多重比較校正應該用 Bonferroni（p 閾值降到 0.05/5 = 0.01）或 Andrews (1993) 的監督邊界。

city-ana 目前版本沒做校正，因為：

• 數據規模小（n=8、候選 5 個），保守的 Bonferroni 過於嚴苛
• Andrews 邊界對單一突變點測試最嚴格，但實作複雜

實務建議：論文裡明確說「未做多重校正」，並對 p ≈ 0.05 的結果保持謹慎。

坑 4：假設斷點兩端都是線性

Chow 假設斷點兩段都是線性。如果整體是「指數增長」，斷點兩段也都是「指數的近似線性化」，可能會找出虛假斷點。

對策：先對 y 取對數 log(y) 再做 Chow。對建築數量這種可能指數增長的指標，對數版的 Chow 更可靠。

---

六、Mann–Kendall + Chow 一起用：兩階段策略

兩個檢定回答不同問題：

步驟	工具	回答
1	Mann–Kendall	整體有沒有趨勢？
2	Chow	趨勢什麼時候改變？

完整工作流：

# 階段 1：先用 MK 篩出「有趨勢」的指標
mk_df = run_mann_kendall(panel)
significant_indicators = mk_df[mk_df["p_value"] < 0.05]["indicator"].tolist()

# 階段 2：對顯著指標跑 Chow
chow_df = run_structural_break(panel, indicators=significant_indicators)

這個策略的好處：

• MK 過濾掉「沒趨勢」的指標 → Chow 不用浪費算力
• 對「有趨勢」的指標，Chow 進一步說「趨勢從何時加速」
• 結合兩者，論文敘述更立體

台中數據的完整結論：

Mann–Kendall 檢定顯示，台中市 2018-2025 年的建築數量、設施、商店、POI 多樣性、紋理熵、邊緣密度、建築覆蓋率均呈現顯著上升趨勢 (p < 0.05)。進一步 Chow 結構突變點檢定揭示，建築數量的增長在 2022 年發生顯著轉折 (F = 28.4, p = 0.001)，從早期斜率約 +7 棟/年 加速至 +22 棟/年；機能類指標（商店、設施、POI 多樣性）的轉折更早（2020-2021），暗示需求面的變化先於開發面的回應。

這段話拿 SCI 期刊評審看，能站得住腳。

---

七、Chow 之外的替代方案

Chow 是最經典的「單一斷點 + 已知時點」框架。實務上還有幾個更現代的選擇：

方法	適用	工具
Chow	單一斷點，掃描候選	本篇
Bai–Perron	多個斷點同時偵測	R strucchange
CUSUM	連續監控（線上偵測）	scipy / R
PELT	通用變點偵測，O(n)	Python ruptures
Bayesian change point	提供斷點機率分佈	PyMC

對 city-ana 這種「8 年小樣本、想找單一加速期」的場景，Chow 已經足夠。如果未來資料拉長到 20+ 年，可以換成 Bai–Perron 找多斷點（譬如 2008 金融危機 + 2020 疫情）。

---

八、視覺化突變點

文字結論之外，最有說服力的是把「兩段擬合」畫出來：

import matplotlib.pyplot as plt
import numpy as np

years = np.array([2018, 2019, 2020, 2021, 2022, 2023, 2024, 2025])
vals  = np.array([27.6, 31.0, 34.8, 47.6, 67.8, 109.8, 131.3, 136.0])

tau = 4  # 2022 對應索引 4

# 兩段線性擬合
def fit_line(t, y):
    X = np.column_stack([np.ones(len(t)), t])
    beta = np.linalg.lstsq(X, y, rcond=None)[0]
    return beta, X @ beta

beta1, fit1 = fit_line(years[:tau], vals[:tau])
beta2, fit2 = fit_line(years[tau:], vals[tau:])

# 整體擬合
beta_full, fit_full = fit_line(years, vals)

plt.figure(figsize=(9, 5))
plt.scatter(years, vals, s=60, c="black", zorder=3, label="觀測")
plt.plot(years, fit_full, "--", color="gray", label=f"單段擬合 (β={beta_full[1]:.1f})")
plt.plot(years[:tau], fit1, "-", color="C0", linewidth=2, label=f"前段 (β={beta1[1]:.1f})")
plt.plot(years[tau:], fit2, "-", color="C3", linewidth=2, label=f"後段 (β={beta2[1]:.1f})")
plt.axvline(years[tau], linestyle=":", color="red", alpha=0.5, label=f"斷點 {years[tau]}")
plt.title("Chow 結構突變點：台中建築數量 (F=28.4, p=0.001)")
plt.xlabel("年份"); plt.ylabel("建築均值")
plt.legend()
plt.tight_layout()
plt.savefig("chow_demo.png", dpi=120)

論文裡這張圖比 F 統計量更有說服力 —— 斜率變化的視覺差異，就是統計檢定的物理意義。

---

九、本篇小結

概念	重點
核心問題	趨勢什麼時候從緩變陡？
Chow 統計量	比較單段 vs 雙段 OLS 的殘差平方和
F 檢定	F = (RSS_full - RSS₁ - RSS₂)/k ÷ (RSS₁ + RSS₂)/(n-2k)
掃描策略	候選斷點 [2, n-2]，取 F 最大值
跟 MK 搭配	MK 篩有趨勢、Chow 找轉折，論文敘述更完整
限制	多重比較未校正、假設兩段線性、小樣本檢力有限

台中數據結論：建築數量在 2022 年發生顯著結構突變，斜率從 +7/年 加速到 +22/年；機能類指標的轉折提前一年（2020-2021），暗示需求驅動先於供給回應。

---

十、系列 C 至此完整：城市發展統計分析三件套

問題	工具	程式碼
是不是空間聚集？	Moran's I	analysis/spatial_autocorr.py
是不是時序趨勢？	Mann–Kendall + Sen	analysis/trend_test.py:30-152
趨勢什麼時候變？	Chow 突變點	analysis/trend_test.py:157-238

三個方法都是非參數 / 對小樣本友好，正好適合「8 年 × 44 點」這種典型的城市快照面板。

接下來的系列 C 還有兩個值得寫的主題：

• C4 — PCA + K-Means 把多指標壓成軌跡聚類
• C5 — Pearson 跨資料源交叉驗證

但這些可以晚一點再寫。下一篇優先補回前面跳過的 A2 — Shannon Entropy 算 POI 多樣性，把系列 A 補完整。

---

程式碼：analysis/trend_test.py:157-288
輸出：output/structural_break.csv
進階閱讀：

• Chow, G. C. (1960). "Tests of equality between sets of coefficients in two linear regressions". Econometrica.
• Bai, J., & Perron, P. (2003). "Computation and analysis of multiple structural change models". Journal of Applied Econometrics.