""" L08 (v2): 雨量フィールドの PCA — 主成分は「時間構造」を捉えるか「空間構造」を捉えるか 14日分の10分雨量 (2024-06-29 〜 07-12, DoBoX #1275) を **観測所 × 時間 (約 280 × 2016)** の行列にし、PCA で主要モードを抽出。 loadings / scores / 再構成誤差 / NMF と多角的に比較し、 PCA が時間カーネルを捉えるのか、空間クラスタを捉えるのかを検証する。 実行: cd "2026 DoBoX 教材" py -X utf8 lessons/L08_pca_rain.py データは存在しなければ DoBoX から自動取得する (ensure_dataset)。 """ from pathlib import Path import re import numpy as np import pandas as pd import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt from matplotlib.colors import ListedColormap from sklearn.preprocessing import StandardScaler from sklearn.decomposition import PCA, NMF from _common import (ROOT, ASSETS, LESSONS, render_lesson, code, figure, data_recipe, parse_rain_csv, ensure_dataset) plt.rcParams["font.family"] = "Yu Gothic" plt.rcParams["axes.unicode_minus"] = False RNG = np.random.default_rng(42) # === 0. データ自動取得 ======================================================= RAIN_DIR = ROOT / "data" / "rain_2024" RAIN_DIR.mkdir(parents=True, exist_ok=True) RA_PATH = ROOT / "data" / "extras" / "rainfall_annual.csv" print("=== 0. データ自動取得 ===") ensure_dataset(RA_PATH, dataset_id=1276, label="雨量年集計 #1276 (メタ用)") RAIN_RID = { "2024-06-29": 94490, "2024-06-30": 94495, "2024-07-01": 94500, "2024-07-02": 94505, "2024-07-03": 94510, "2024-07-04": 94515, "2024-07-05": 94520, "2024-07-06": 94525, "2024-07-07": 94530, "2024-07-08": 94535, "2024-07-09": 94540, "2024-07-10": 94545, "2024-07-11": 94551, "2024-07-12": 94556, } for date, rid in RAIN_RID.items(): ensure_dataset(RAIN_DIR / f"rain_{date}.csv", resource_id=rid, min_bytes=500_000, label=f"雨量10分値 {date}") # === 1. 14日分10分値を station×time 行列にする ================================= print("\n=== 1. 14日分10分値を station×time 行列にする ===") rain_files = sorted(RAIN_DIR.glob("rain_2024-*.csv")) print(f" 雨量ファイル: {len(rain_files)}件") rain_dfs = [parse_rain_csv(f) for f in rain_files] rain_full = pd.concat(rain_dfs).sort_index() print(f" 10分値 raw shape (時刻×観測所): {rain_full.shape}") # 行=観測所, 列=時刻 に転置 (PCA では行=サンプル=観測所, 列=特徴量=時刻) M = rain_full.T # (station, time) # 全期間で値が有効な観測所のみ残す valid_mask = M.notna().mean(axis=1) >= 0.95 # 欠損 5% 以下 M = M.loc[valid_mask].fillna(0.0) # 完全に降雨ゼロの観測所は除外 (PCA は分散ゼロを嫌う) totals = M.sum(axis=1) M = M.loc[totals > 0.5] print(f" 有効観測所 × 時刻数 行列: {M.shape}") # === 2. 観測所メタ (事務所名) を結合: 空間プロキシ ============================= print("\n=== 2. 観測所メタ (事務所名) を空間プロキシとして結合 ===") ra_raw = pd.read_csv(RA_PATH, header=None, encoding="utf-8-sig") # rainfall_annual.csv の段組: # row 0: 事務所名 / row 1: データ所管 / row 2: 水系名 # row 3: 河川名 / row 4: 観測所名 / row 5: 観測時刻 (=雨量(mm)) meta = pd.DataFrame({ "office": ra_raw.iloc[0, 1:].astype(str).tolist(), "water_system": ra_raw.iloc[2, 1:].astype(str).tolist(), "station": ra_raw.iloc[4, 1:].astype(str).tolist(), }) office_map = {} ws_map = {} for _, r in meta.iterrows(): s = r["station"] if s and s != "nan" and s not in office_map: office_map[s] = r["office"] ws_map[s] = r["water_system"] def lookup(col, mp): if col in mp: return mp[col] base = re.sub(r"_\d+$", "", col) return mp.get(base, "不明") station_office = np.array([lookup(s, office_map) for s in M.index]) station_ws = np.array([lookup(s, ws_map) for s in M.index]) print(f" 事務所別観測所数 (上位):") office_counts = pd.Series(station_office).value_counts() print(office_counts.head(10).to_string()) # === 3. 標準化 + PCA (k=10) ================================================== print("\n=== 3. 標準化 + PCA (k=10) ===") # 各観測所 (= 行) を標準化: 観測所間の "降雨スケール" を揃え、降雨の "形" を比較 Xc = M.values - M.values.mean(axis=1, keepdims=True) # row-centered # 観測所の総雨量による相対化 (zero-rain 既に除外済) — std で割る row_std = Xc.std(axis=1, keepdims=True) row_std[row_std == 0] = 1.0 X = Xc / row_std K = 10 pca = PCA(n_components=K, random_state=42) Z = pca.fit_transform(X) # (n_station, K) EVR = pca.explained_variance_ratio_ print(f" X.shape={X.shape}, EVR(top5)={[f'{e:.3f}' for e in EVR[:5]]}, " f"cum(top5)={np.cumsum(EVR)[:5].round(3).tolist()}") # === 4. 図1: scree plot (寄与率 + 累積) ======================================= print("\n=== 4. 図1: scree plot ===") fig, ax = plt.subplots(figsize=(8.5, 4.6)) xs = np.arange(1, K + 1) ax.bar(xs, EVR * 100, color="#0969da", alpha=0.85, label="個別寄与率") ax2 = ax.twinx() ax2.plot(xs, np.cumsum(EVR) * 100, "o-", color="#cf222e", linewidth=2, markersize=8, label="累積寄与率") ax2.axhline(80, color="#888", linestyle="--", linewidth=1) ax2.text(K, 81, "80%", color="#888", ha="right", fontsize=10) ax.set_xticks(xs) ax.set_xlabel("主成分番号 k") ax.set_ylabel("個別寄与率 (%)", color="#0969da") ax2.set_ylabel("累積寄与率 (%)", color="#cf222e") ax2.set_ylim(0, 100) ax.tick_params(axis="y", labelcolor="#0969da") ax2.tick_params(axis="y", labelcolor="#cf222e") ax.set_title(f"PCA scree plot — 観測所×時間行列 {X.shape} " f"(PC1={EVR[0]*100:.1f}%, 上位3で {np.sum(EVR[:3])*100:.1f}%)") ax.grid(alpha=0.3) plt.tight_layout() plt.savefig(ASSETS / "L08_scree.png", dpi=140) plt.close() # === 5. 図2: PC1〜PC3 の loadings を時系列としてプロット ===================== print("\n=== 5. 図2: PC1-PC3 loadings (時系列) ===") times = pd.to_datetime(M.columns) fig, axes = plt.subplots(3, 1, figsize=(13.5, 7.5), sharex=True) colors = ["#cf222e", "#0969da", "#1f883d"] for i, (ax, col) in enumerate(zip(axes, colors)): comp = pca.components_[i] ax.fill_between(times, 0, comp, color=col, alpha=0.35) ax.plot(times, comp, color=col, linewidth=0.8) ax.axhline(0, color="#666", linewidth=0.6) ax.set_ylabel(f"PC{i+1} loading\n({EVR[i]*100:.1f}%)", fontsize=10) ax.grid(alpha=0.25) # 主要降雨イベントの境界 (日付マーカー) for d in pd.date_range("2024-06-29", "2024-07-12", freq="D"): ax.axvline(d, color="#bbb", linestyle=":", linewidth=0.5) axes[0].set_title("主成分の loadings (=時間方向の重みベクトル)\n" "→ 強いピークがある日 = その PC が表す『時間モード』") axes[-1].set_xlabel("時刻 (2024)") plt.tight_layout() plt.savefig(ASSETS / "L08_loadings.png", dpi=140) plt.close() # loadings の解釈統計: 各 PC のピーク時刻 print(" 各 PC のピーク時刻 (絶対値最大):") for i in range(3): abs_load = np.abs(pca.components_[i]) pk = int(np.argmax(abs_load)) print(f" PC{i+1}: peak at {times[pk]} (loading={pca.components_[i][pk]:.3f})") # === 6. 図3: PC1×PC2 観測所スコア散布 (事務所別カラー) ======================= print("\n=== 6. 図3: PC1×PC2 観測所スコア (事務所別=空間プロキシ) ===") # 事務所を上位 8 + その他にまとめる top_offices = office_counts.head(8).index.tolist() office_label = np.where(np.isin(station_office, top_offices), station_office, "その他") unique_offices = list(dict.fromkeys(top_offices + ["その他"])) palette = ["#cf222e", "#0969da", "#1f883d", "#fb8500", "#8250df", "#bf3989", "#16a085", "#d4a72c", "#888888"] o2c = {o: palette[i % len(palette)] for i, o in enumerate(unique_offices)} fig, ax = plt.subplots(figsize=(9.5, 7.0)) for office in unique_offices: mask = (office_label == office) if not mask.any(): continue ax.scatter(Z[mask, 0], Z[mask, 1], c=o2c[office], s=42, alpha=0.78, edgecolor="black", linewidth=0.4, label=f"{office} (n={mask.sum()})") ax.axhline(0, color="#999", linewidth=0.5) ax.axvline(0, color="#999", linewidth=0.5) ax.set_xlabel(f"PC1 score ({EVR[0]*100:.1f}%)") ax.set_ylabel(f"PC2 score ({EVR[1]*100:.1f}%)") ax.set_title("観測所スコア PC1 × PC2 色 = 事務所 (空間プロキシ)\n" "→ 同じ事務所の観測所が固まれば PCA は『空間構造』を捉えている") ax.legend(fontsize=9, loc="best", framealpha=0.85, ncol=2) ax.grid(alpha=0.3) plt.tight_layout() plt.savefig(ASSETS / "L08_scatter_office.png", dpi=140) plt.close() # 事務所による PC スコアの分散 (空間構造の指標) def variance_explained_by_group(scores, groups): """1元配置分散分析の説明分散比 η² = SSB / SST.""" grand = scores.mean() sst = ((scores - grand) ** 2).sum() if sst == 0: return 0.0 ssb = 0.0 for g in np.unique(groups): mask = (groups == g) if mask.sum() < 2: continue ssb += mask.sum() * (scores[mask].mean() - grand) ** 2 return float(ssb / sst) print("\n 事務所による PC スコアの説明分散比 η² (空間構造の強さ):") for i in range(5): eta_o = variance_explained_by_group(Z[:, i], station_office) eta_w = variance_explained_by_group(Z[:, i], station_ws) print(f" PC{i+1}: η²(事務所)={eta_o:.3f}, η²(水系)={eta_w:.3f}") # === 7. 図4: 再構成誤差 vs k ================================================= print("\n=== 7. 図4: 再構成誤差 vs k ===") rmse_k = [] mae_k = [] for k in range(1, K + 1): Xk_hat = Z[:, :k] @ pca.components_[:k] err = X - Xk_hat rmse_k.append(float(np.sqrt(np.mean(err ** 2)))) mae_k.append(float(np.mean(np.abs(err)))) print(f" RMSE(k=1)={rmse_k[0]:.3f}, RMSE(k=3)={rmse_k[2]:.3f}, " f"RMSE(k=5)={rmse_k[4]:.3f}, RMSE(k=10)={rmse_k[-1]:.3f}") # 例として代表観測所の再構成カーブ total_rain = M.sum(axis=1) top_idx = int(np.argmax(total_rain.values)) top_station = M.index[top_idx] fig, axes = plt.subplots(1, 2, figsize=(14.5, 4.8)) # 左: 全体の再構成誤差曲線 ax = axes[0] ax.plot(range(1, K + 1), rmse_k, "o-", color="#cf222e", linewidth=2, label="RMSE") ax.plot(range(1, K + 1), mae_k, "s--", color="#0969da", linewidth=1.5, label="MAE") ax.set_xticks(range(1, K + 1)) ax.set_xlabel("使用主成分数 k") ax.set_ylabel("再構成誤差 (標準化スケール)") ax.set_title("再構成誤差 vs 主成分数") ax.legend() ax.grid(alpha=0.3) # 右: 最大降雨観測所の k=1, 3, 10 の再構成 ax = axes[1] ax.plot(times, X[top_idx], color="black", linewidth=1.0, label="原データ (標準化)", alpha=0.8) for k, c in zip([1, 3, 10], ["#cf222e", "#fb8500", "#1f883d"]): Xk_hat = Z[:, :k] @ pca.components_[:k] ax.plot(times, Xk_hat[top_idx], color=c, linewidth=1.5, alpha=0.85, label=f"k={k} 再構成 (RMSE={rmse_k[k-1]:.2f})") ax.set_title(f"最大降雨観測所 「{top_station}」 の再構成 (累計={total_rain.iloc[top_idx]:.0f}mm)") ax.set_xlabel("時刻") ax.set_ylabel("標準化雨量") ax.legend(fontsize=9, loc="best") ax.grid(alpha=0.3) ax.tick_params(axis="x", rotation=20) plt.tight_layout() plt.savefig(ASSETS / "L08_reconstruction.png", dpi=140) plt.close() # === 8. 図5: NMF と比較 (非負成分による分解) ================================= print("\n=== 8. 図5: NMF と比較 ===") # NMF は非負入力が必要なので生の雨量 (標準化なし) を使う X_pos = M.values.copy() # サンプル別正規化 (各観測所の総雨量で割る) → "降雨パターン" を分解 row_sums = X_pos.sum(axis=1, keepdims=True) row_sums[row_sums == 0] = 1.0 X_pos_n = X_pos / row_sums K_nmf = 3 nmf = NMF(n_components=K_nmf, init="nndsvda", max_iter=500, random_state=42, tol=1e-4) W = nmf.fit_transform(X_pos_n) # (n_station, K_nmf) 非負スコア H = nmf.components_ # (K_nmf, n_time) 非負基底パターン print(f" NMF reconstruction err = {nmf.reconstruction_err_:.3f}") fig, axes = plt.subplots(2, 1, figsize=(13.5, 7.0)) # 上: NMF 基底 (時間方向) ax = axes[0] nmf_colors = ["#cf222e", "#0969da", "#1f883d"] for i in range(K_nmf): ax.fill_between(times, 0, H[i], color=nmf_colors[i], alpha=0.35) ax.plot(times, H[i], color=nmf_colors[i], linewidth=1.0, label=f"NMF 基底 {i+1}") for d in pd.date_range("2024-06-29", "2024-07-12", freq="D"): ax.axvline(d, color="#bbb", linestyle=":", linewidth=0.5) ax.legend(loc="upper right", fontsize=9) ax.set_ylabel("基底成分 (非負)") ax.set_title("NMF basis: 各成分は『時間集中型の降雨イベント』として現れる\n" "(PCA loadings と違って符号反転がなく、解釈が直接的)") ax.grid(alpha=0.25) # 下: 各観測所の最も強く割り当てられた NMF 基底でカラー ax = axes[1] W_norm = W / (W.sum(axis=1, keepdims=True) + 1e-9) nmf_label = W_norm.argmax(axis=1) for i in range(K_nmf): mask = (nmf_label == i) if not mask.any(): continue ax.scatter(Z[mask, 0], Z[mask, 1], c=nmf_colors[i], s=42, alpha=0.78, edgecolor="black", linewidth=0.4, label=f"NMF 基底{i+1} 主導 (n={mask.sum()})") ax.set_xlabel(f"PC1 ({EVR[0]*100:.1f}%)") ax.set_ylabel(f"PC2 ({EVR[1]*100:.1f}%)") ax.set_title("PCA 散布上に NMF の主基底クラスを重ね描き\n" "→ NMF は『どの降雨イベントで濡れたか』で観測所をグループ化") ax.axhline(0, color="#999", linewidth=0.5) ax.axvline(0, color="#999", linewidth=0.5) ax.legend(fontsize=9, loc="best") ax.grid(alpha=0.3) plt.tight_layout() plt.savefig(ASSETS / "L08_nmf_compare.png", dpi=140) plt.close() # === 9. サマリ表 ============================================================ print("\n=== 9. サマリ表 ===") summary_rows = [] for i in range(5): eta_o = variance_explained_by_group(Z[:, i], station_office) eta_w = variance_explained_by_group(Z[:, i], station_ws) abs_load = np.abs(pca.components_[i]) pk = int(np.argmax(abs_load)) summary_rows.append({ "PC": f"PC{i+1}", "寄与率(%)": round(EVR[i] * 100, 2), "累積(%)": round(np.cumsum(EVR)[i] * 100, 2), "loading ピーク時刻": str(times[pk])[:16], "η² 事務所(空間)": round(eta_o, 3), "η² 水系(空間)": round(eta_w, 3), }) summary_df = pd.DataFrame(summary_rows) summary_html = summary_df.to_html(index=False, na_rep="-") print(summary_df.to_string(index=False)) recon_rows = [] for k in [1, 2, 3, 5, 10]: recon_rows.append({ "k": k, "累積寄与率(%)": round(np.cumsum(EVR)[k-1] * 100, 2), "RMSE": round(rmse_k[k-1], 4), "MAE": round(mae_k[k-1], 4), }) recon_df = pd.DataFrame(recon_rows) recon_html = recon_df.to_html(index=False) print(recon_df.to_string(index=False)) # === 10. HTML レンダリング =================================================== print("\n=== 10. HTML レンダリング ===") sections = [ ("学習目標", """ """), ("使用データ", """

14日窓には 2024年7月豪雨の前駆〜本格降雨 (7/1, 7/10〜7/12) が含まれており、PCA で時間モードが立ちやすい好題材。

"""), ("データ取得手順", data_recipe([ {"label": "雨量10分値 14日分 (各日別 CSV)", "dataset_id": 1275, "file": "data/rain_2024/rain_2024-MM-DD.csv (×14)", "format": "CSV (5段ヘッダ, 10分値, UTF-8 BOM)", "size": "1.0〜1.4 MB / 日"}, {"label": "雨量 年集計 (観測所→事務所/水系)", "dataset_id": 1276, "file": "data/extras/rainfall_annual.csv", "format": "CSV (多段ヘッダ)", "size": "約 500 KB"}, ])), ("方法 — 行列化の選択", """

14日分の10分雨量を PCA で解析するために、行と列をどう割り当てるかがまず最大の意思決定になる。

選択肢形状主成分が捉えるもの欠点
(A) 観測所 × 時間 (10分粒度)≈ 280 × 2016時間モード (県全体に共通する降雨イベント) と 地域差 の両方列数 >> 行数 → 共分散の安定性に注意
(B) 観測所 × 日 (14日)≈ 280 × 14日単位の粗い時間モード14列では分散が薄まり、降雨イベントの形状が潰れる (旧 v1 教材の限界)
(C) 時間 × 観測所 (転置)≈ 2016 × 280各時刻の "雨域パターン" の主成分サンプル間が時系列依存で独立性仮定が崩れる

本教材では (A) 観測所 × 時間 (2016列) を採用する。理由:

  1. 降雨イベントの形 (10分単位の立ち上がり/減衰) を保存できる → 第1主成分が「7月豪雨イベントそのもの」として現れることを確認できる
  2. n=280 と p=2016 は p>n だが、PCA は SVD で安定計算可能 (rank ≤ min(n,p)−1 = 279)
  3. 各観測所を row-center + row-scale で標準化 → 「降雨の総量」ではなく「降雨パターンの形」を比較する
  1. 14ファイル統合: parse_rain_csv() で 14 日分を tidy 化し転置
  2. 欠損 5% 超 / 全期間ゼロの観測所を除外
  3. 各行 (観測所) を z-score 化: パターン比較のため
  4. PCA(n_components=10): scree plot, loadings, scores, 再構成誤差 を計算
  5. 事務所名を空間プロキシとして scatter を色分け、η²(分散比) で空間構造の強さを定量化
  6. NMF(n_components=3): 非負入力で同じデータを分解し、PCA との解釈差を比較
"""), ("コード解説", code(''' import re, numpy as np, pandas as pd from pathlib import Path import matplotlib.pyplot as plt from sklearn.decomposition import PCA, NMF from _common import parse_rain_csv, ensure_dataset # === (0) 14日分データを取得 === RAIN_RID = {"2024-06-29": 94490, ..., "2024-07-12": 94556} for date, rid in RAIN_RID.items(): ensure_dataset(f"data/rain_2024/rain_{date}.csv", resource_id=rid) # === (1) 観測所×時間 行列を作る === files = sorted(Path("data/rain_2024").glob("rain_2024-*.csv")) rain = pd.concat([parse_rain_csv(f) for f in files]).sort_index() # (時刻, 観測所) M = rain.T # (観測所, 時刻) M = M.loc[M.notna().mean(axis=1) >= 0.95].fillna(0.0) M = M.loc[M.sum(axis=1) > 0.5] # 全期間ゼロ除外 print(M.shape) # 例: (271, 2016) # === (2) 行ごとに z-score 標準化 (観測所間の "降雨スケール" を揃える) === Xc = M.values - M.values.mean(axis=1, keepdims=True) row_std = Xc.std(axis=1, keepdims=True); row_std[row_std==0] = 1.0 X = Xc / row_std # === (3) PCA(k=10) === pca = PCA(n_components=10, random_state=42) Z = pca.fit_transform(X) # (n_station, 10) 観測所スコア EVR = pca.explained_variance_ratio_ # 寄与率 print("PC1 =", EVR[0], " cum top3 =", EVR[:3].sum()) # === (4) loadings は時間方向のベクトル → そのまま時系列としてプロット === times = pd.to_datetime(M.columns) for i in range(3): plt.plot(times, pca.components_[i], label=f"PC{i+1}") # loadings の "ピーク時刻" を見れば、その PC が表す時間モードが分かる # 例: PC1 のピークが 7/2 早朝 → PC1 = 「2024-07-02 豪雨イベント」 # === (5) 観測所スコア散布: 事務所別カラーで空間構造を見る === office_for_station = lookup_office_from_rainfall_annual(...) plt.scatter(Z[:,0], Z[:,1], c=office_for_station) # 事務所が PC スコアでどれだけ説明されるかを η² で測る def eta2(scores, groups): grand = scores.mean() sst = ((scores - grand)**2).sum() ssb = sum(m.sum()*(scores[m].mean()-grand)**2 for g in np.unique(groups) for m in [groups==g]) return ssb / sst # === (6) 再構成誤差 vs k === rmse = [np.sqrt(((X - Z[:,:k] @ pca.components_[:k])**2).mean()) for k in range(1, 11)] # k=3 で何 % 再現できるかを RMSE で見る # === (7) NMF 比較: 非負入力 → 解釈しやすい "降雨イベント基底" になる === X_pos_n = M.values / M.values.sum(axis=1, keepdims=True) # 各観測所で正規化 nmf = NMF(n_components=3, init="nndsvda", max_iter=500, random_state=42) W = nmf.fit_transform(X_pos_n) # 観測所 × 3 (非負スコア) H = nmf.components_ # 3 × 時間 (非負基底=降雨イベント形状) # H[i] が「i 番目の典型的な降雨イベント時系列」、W[s,i] が「観測所 s が i に # どれだけ似ているか」 — PCA の符号反転問題が無く解釈が一意 ''')), ("結果", figure("assets/L08_scree.png", f"scree plot — PC1={EVR[0]*100:.1f}%, 上位3で {np.sum(EVR[:3])*100:.1f}%, " f"上位5で {np.sum(EVR[:5])*100:.1f}%, 上位10で {np.sum(EVR[:10])*100:.1f}%。" "1モードに集中する豪雨期間ながら、空間多様性ゆえ累積80%には10成分前後を要する。") + figure("assets/L08_loadings.png", "PC1〜PC3 の loadings を時間軸でプロット。点線は日付境界。" "loading が強くピークする時刻 = その PC が表す時間モード。") + figure("assets/L08_scatter_office.png", "PC1×PC2 観測所スコアの散布図。色 = 事務所 (空間プロキシ)。" "同じ事務所が固まって見えれば PCA は『空間構造』も捉えている。") + figure("assets/L08_reconstruction.png", "(左) RMSE/MAE と k の関係。(右) 最大降雨観測所の波形を k=1,3,10 で再構成。" "k=3 でピークの位置を、k=10 で細部の波形を概ね再現。") + figure("assets/L08_nmf_compare.png", "(上) NMF が抽出した3つの非負基底 = 異なるタイミングの降雨イベント。" "(下) 各観測所をその主導 NMF 基底でラベル付けし PCA 散布上に重ね描き。") + "

PC ごとのサマリ (寄与率 + 空間構造の強さ η²)

" + summary_html + "

再構成誤差 vs k

" + recon_html), ("考察", f""" """), ("発展課題", """
  1. 観測所緯度経度を結合 (DoBoX 別データセットや観測所マスタを別途取得) して、PC2/PC3 のスコアを地図にカラーマップする。 本教材では事務所名を空間プロキシとしたが、緯度経度があれば 北東-南西の降雨グラデーション が可視化できる。
  2. 列を中心化のみ・行は中心化のみ・両方の3条件で PCA を再実行し、結果がどう変わるかを比較する (前処理の影響学習)。
  3. EOF (Empirical Orthogonal Function) 解析: 気象学で使われる時空間PCAの呼称。本教材の手順と全く同じだが、地球流体力学の文脈での解釈を学ぶ。
  4. ICA (Independent Component Analysis) を試す: PCA の「直交性」より「統計的独立性」を優先する分解。降雨イベントが独立な 'ソース' から来ていると仮定するなら ICA の方が物理的に合う。
  5. テンソル分解 (Tucker / CP): (観測所, 日, 時刻) の3階テンソルとして直接分解すれば、日内パターンと日間パターンを同時に分離できる。
  6. PCA で外れ値検出: 再構成誤差が他観測所より大きい観測所は「特異な雨域」「センサ異常」「地形効果」のいずれか。本データで誰が外れ値かを調べる。
"""), ] html = render_lesson( num=8, title="PCA で雨域構造を読み解く — 時間モード vs 空間モード", tags=["選択", "応用基礎", "次元削減", "v2-rewrite"], time="120分", level="応用基礎", data_label='2データ横断: ' '#1275 (雨量10分値) + ' '#1276 (観測所メタ)', sections=sections, script_filename="lessons/L08_pca_rain.py", ) out = LESSONS / "L08_pca_rain.html" out.write_text(html, encoding="utf-8") print(f"\nsaved: {out}") print(f" PNGs: 5 generated in {ASSETS}") print(f" X.shape = {X.shape}, top3 EVR = {EVR[:3].round(3).tolist()}")