地方創生を推進する多角的複合指標の提案
―SVMに基づく主観的でない変数選択と重み付け手法の検討―

print("\n" + "=" * 60)
print("■ Step3. 複合指標（SVM重み付き）の構築")
print("=" * 60)

weights_norm   = svm_weights / svm_weights.sum()
composite_index = X_scaled @ weights_norm   # (47,)

rank_df = pd.DataFrame({
    '都道府県':     PREFS,
    '複合指標（SVM）': composite_index,
}).sort_values('複合指標（SVM）', ascending=False).reset_index(drop=True)
rank_df['順位'] = range(1, len(PREFS) + 1)

print("【複合指標ランキング（上位10）】")
print(rank_df.head(10).to_string(index=False))

print("\n図1: SVM重み棒グラフを作成中...")

n_vars = len(VAR_NAMES)
colors_bar = ['#1565C0' if i < 3 else '#90CAF9' for i in range(n_vars)]

fig1, ax1 = plt.subplots(figsize=(10, 6))
# weight_df は降順 → 横棒グラフは上から重み大の順になるよう逆転
bars = ax1.barh(
    weight_df['変数'][::-1],
    weight_df['SVM重み（絶対値）'][::-1],
    color=colors_bar[::-1],
    edgecolor='white',
    alpha=0.88,
)
mean_w = weight_df['SVM重み（絶対値）'].mean()
ax1.axvline(mean_w, color='gray', linestyle='--', linewidth=1.2,
            label=f'平均 = {mean_w:.3f}')
ax1.set_xlabel('LinearSVR係数の絶対値（変数重要度）', fontsize=12)
ax1.set_title(
    'SVM（線形カーネル）による変数重要度\n地方創生関連変数のランキング（2022年度, 47都道府県）',
    fontsize=13, fontweight='bold',
)
ax1.legend(fontsize=10)
ax1.grid(axis='x', alpha=0.3)
for bar, val in zip(bars, weight_df['SVM重み（絶対値）'][::-1]):
    ax1.text(val + mean_w * 0.02, bar.get_y() + bar.get_height() / 2,
             f'{val:.3f}', va='center', fontsize=9, color='#333')

plt.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2024_U5_2_fig1_svm_weights.png'),
             bbox_inches='tight', dpi=150)
plt.close(fig1)
print("  → 2024_U5_2_fig1_svm_weights.png 保存完了")

print("図2: 複合指標ランキングを作成中...")

top15 = rank_df.head(15).copy()
bot15 = rank_df.tail(15).iloc[::-1].copy()

# 上位15の最高スコア都道府県を強調
top1_name = rank_df.iloc[0]['都道府県']
last_name  = rank_df.iloc[-1]['都道府県']

fig2, axes2 = plt.subplots(1, 2, figsize=(14, 7))
fig2.suptitle('地方創生度 複合指標ランキング（LinearSVR重み付き, 2022年度）',
              fontsize=13, fontweight='bold')

for ax, sub_df, title, clr_hi, clr_lo in zip(
    axes2,
    [top15, bot15],
    ['上位15都道府県', '下位15都道府県'],
    ['#1565C0', '#C62828'],
    ['#90CAF9',  '#EF9A9A'],
):
    bar_colors = [
        clr_hi if row['都道府県'] == top1_name or row['都道府県'] == last_name
        else clr_lo
        for _, row in sub_df.iterrows()
    ]
    ax.barh(sub_df['都道府県'], sub_df['複合指標（SVM）'],
            color=bar_colors, edgecolor='white', alpha=0.88)
    ax.axvline(0, color='black', linewidth=0.8)
    ax.set_xlabel('複合指標スコア（標準化加重和）', fontsize=11)
    ax.set_title(title, fontsize=12, fontweight='bold')
    ax.grid(axis='x', alpha=0.3)

plt.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2024_U5_2_fig2_index_rank.png'),
             bbox_inches='tight', dpi=150)
plt.close(fig2)
print("  → 2024_U5_2_fig2_index_rank.png 保存完了")

import os
import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings('ignore')

from sklearn.svm import LinearSVR
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from scipy import stats as scipy_stats

# ── パス設定 ─────────────────────────────────────────────────────────────────
BASE_DIR = os.path.join(_script_dir, '..')
FIG_DIR  = os.path.join(BASE_DIR, 'html', 'figures')
DATA_PATH = os.path.join(BASE_DIR, 'data', 'raw', 'SSDSE-B-2026.csv')
os.makedirs(FIG_DIR, exist_ok=True)

plt.rcParams.update({
    'font.family':        'Hiragino Sans',
    'axes.unicode_minus': False,
    'figure.dpi':         150,
    'axes.spines.top':    False,
    'axes.spines.right':  False,
})

print("=" * 60)
print("■ Step1. SSDSE-B-2026 実データ読み込み（2022年度）")
print("=" * 60)

df_b = pd.read_csv(DATA_PATH, encoding='cp932', header=1)
# 都道府県行のみ（市区町村を除外）
df_b = df_b[df_b['地域コード'].str.match(r'^R\d{5}$', na=False)]
df = df_b[df_b['年度'] == 2022].copy().reset_index(drop=True)
print(f"  読み込み完了: {len(df)} 都道府県")

# 都道府県名を短縮（「県」「都」「府」「道」を保持しつつ表示用に整理）
PREFS = df['都道府県'].tolist()

# ── 変数算出 ──────────────────────────────────────────────────────────────────
pop   = df['総人口']                                      # A1101
tfr   = df['合計特殊出生率']                             # A4103  ← 目的変数
inflow    = df['転入者数（日本人移動者）']               # A5101
outflow   = df['転出者数（日本人移動者）']               # A5102
nursery_n = df['保育所等数']                             # J2503
nursery_c = df['保育所等定員数']                         # J2505
univ_grads= df['高等学校卒業者のうち進学者数']          # E4602
hs_grads  = df['高等学校卒業者数']                       # E4601
consumption= df['消費支出（二人以上の世帯）']           # L3221
land_price = df['標準価格（平均価格）（住宅地）']       # C5401
job_open   = df['月間有効求人数（一般）']                # F3103
job_seek   = df['月間有効求職者数（一般）']              # F3102
elderly    = df['65歳以上人口']                          # A1303
hospitals  = df['一般病院数']                            # I510120
tourists   = df['延べ宿泊者数']                          # G7101

# 派生指標算出
in_rate      = inflow    / pop * 100         # 転入率 (%)
out_rate     = outflow   / pop * 100         # 転出率 (%)
nursery_dens = nursery_n / pop * 10000       # 保育所密度 (10000人当たり)
nursery_cap  = nursery_c / pop * 10000       # 保育定員密度 (10000人当たり)
univ_rate    = univ_grads / hs_grads * 100   # 高校→大学進学率 (%)
job_ratio    = job_open   / job_seek         # 有効求人倍率
aging_rate   = elderly    / pop * 100        # 高齢化率 (%)
hosp_dens    = hospitals  / pop * 10000      # 一般病院密度 (10000人当たり)
tourist_pc   = tourists   / pop             # 宿泊者数per capita

# 説明変数・変数名を整理
VAR_NAMES = [
    '転入率(%)',
    '転出率(%)',
    '保育所密度',
    '保育定員密度',
    '大学進学率(%)',
    '消費支出(万円)',
    '住宅地標準価格',
    '有効求人倍率',
    '高齢化率(%)',
    '病院密度',
    '宿泊者数/人',
]

X_raw = np.column_stack([
    in_rate.values,
    out_rate.values,
    nursery_dens.values,
    nursery_cap.values,
    univ_rate.values,
    consumption.values / 10000,   # 万円単位にスケール調整
    land_price.values,
    job_ratio.values,
    aging_rate.values,
    hosp_dens.values,
    tourist_pc.values,
])

y = tfr.values   # 目的変数：合計特殊出生率

print(f"  説明変数 ({len(VAR_NAMES)}個): {VAR_NAMES}")
print(f"  目的変数: 合計特殊出生率（TFR）, 範囲 {y.min():.2f} – {y.max():.2f}")

print("\n" + "=" * 60)
print("■ Step2. LinearSVR変数重み算出")
print("=" * 60)

scaler   = StandardScaler()
X_scaled = scaler.fit_transform(X_raw)

svr = LinearSVR(C=1.0, max_iter=10000)
svr.fit(X_scaled, y)

svm_weights = np.abs(svr.coef_)   # shape: (n_features,)
weight_df = pd.DataFrame({
    '変数':          VAR_NAMES,
    'SVM重み（絶対値）': svm_weights,
}).sort_values('SVM重み（絶対値）', ascending=False).reset_index(drop=True)

print("【SVM変数重み（降順）】")
print(weight_df.round(4).to_string(index=False))

print("\n" + "=" * 60)
print("■ Step4. PCA（主成分分析）との比較")
print("=" * 60)

pca = PCA(n_components=2)
X_pca = pca.fit_transform(X_scaled)
pca_index = X_pca[:, 0]   # 第1主成分をPCA指標として使用

print(f"  PCA説明分散比: PC1={pca.explained_variance_ratio_[0]:.3f}, "
      f"PC2={pca.explained_variance_ratio_[1]:.3f}")

print("図3: PCAバイプロットを作成中...")

# 都道府県ラベル表示: 複合指標上位10・下位10
top10_prefs = set(rank_df.head(10)['都道府県'].tolist())
bot10_prefs = set(rank_df.tail(10)['都道府県'].tolist())
label_prefs = top10_prefs | bot10_prefs

fig3, ax3 = plt.subplots(figsize=(11, 8))
sc = ax3.scatter(
    X_pca[:, 0], X_pca[:, 1],
    c=composite_index, cmap='RdYlBu',
    s=65, alpha=0.82, zorder=3,
)
plt.colorbar(sc, ax=ax3, label='SVM複合指標スコア')

for i, pref in enumerate(PREFS):
    if pref in label_prefs:
        ax3.annotate(pref, (X_pca[i, 0], X_pca[i, 1]),
                     textcoords='offset points', xytext=(6, 3),
                     fontsize=7.5, color='#333')

# ローディングベクトル（矢印）
loadings = pca.components_.T   # shape: (n_features, 2)
scale = 2.2
for j, var in enumerate(VAR_NAMES):
    ax3.annotate(
        '', xy=(loadings[j, 0] * scale, loadings[j, 1] * scale), xytext=(0, 0),
        arrowprops=dict(arrowstyle='->', color='#E65100', lw=1.5),
    )
    ax3.text(
        loadings[j, 0] * scale * 1.12,
        loadings[j, 1] * scale * 1.12,
        var, color='#E65100', fontsize=8,
    )

ax3.axhline(0, color='gray', linewidth=0.5)
ax3.axvline(0, color='gray', linewidth=0.5)
ax3.set_xlabel(f'PC1（説明率 {pca.explained_variance_ratio_[0]*100:.1f}%）', fontsize=12)
ax3.set_ylabel(f'PC2（説明率 {pca.explained_variance_ratio_[1]*100:.1f}%）', fontsize=12)
ax3.set_title(
    'PCAバイプロット（都道府県 × 地方創生変数, 2022年度）\n色：SVM複合指標スコア',
    fontsize=13, fontweight='bold',
)
ax3.grid(True, alpha=0.2)

plt.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2024_U5_2_fig3_pca.png'),
             bbox_inches='tight', dpi=150)
plt.close(fig3)
print("  → 2024_U5_2_fig3_pca.png 保存完了")

print("図4: SVM指標 vs PCA指標散布図を作成中...")

# 両指標を標準化（方向を揃えるため）
svm_std = (composite_index - composite_index.mean()) / composite_index.std()
pca_std = (pca_index - pca_index.mean()) / pca_index.std()

# PCAとSVM指標の方向が逆の場合は符号を揃える（正の相関になるよう）
if scipy_stats.pearsonr(svm_std, pca_std)[0] < 0:
    pca_std = -pca_std

r_val, p_val = scipy_stats.pearsonr(svm_std, pca_std)

# 上位・下位5都道府県をラベル表示
top5_prefs  = set(rank_df.head(5)['都道府県'].tolist())
bot5_prefs  = set(rank_df.tail(5)['都道府県'].tolist())
label5_prefs = top5_prefs | bot5_prefs

scatter_c = [
    '#E53935' if p in top5_prefs else
    '#1E88E5' if p in bot5_prefs else
    '#BDBDBD'
    for p in PREFS
]

fig4, ax4 = plt.subplots(figsize=(9, 7))
ax4.scatter(svm_std, pca_std, c=scatter_c, s=70, alpha=0.85, zorder=3)

for i, pref in enumerate(PREFS):
    if pref in label5_prefs:
        ax4.annotate(pref, (svm_std[i], pca_std[i]),
                     textcoords='offset points', xytext=(7, 3),
                     fontsize=9.5, fontweight='bold',
                     color='#C62828' if pref in top5_prefs else '#0D47A1')

# 回帰直線
x_fit = np.linspace(svm_std.min(), svm_std.max(), 200)
coef  = np.polyfit(svm_std, pca_std, 1)
ax4.plot(x_fit, np.polyval(coef, x_fit),
         color='#444', linewidth=1.8, linestyle='--')

ax4.axhline(0, color='gray', linewidth=0.5)
ax4.axvline(0, color='gray', linewidth=0.5)
ax4.set_xlabel('SVM複合指標（標準化）', fontsize=12)
ax4.set_ylabel('PCA第1主成分スコア（標準化）', fontsize=12)

p_str = '<0.001' if p_val < 0.001 else f'={p_val:.3f}'
ax4.set_title(
    f'SVM複合指標 vs PCA第1主成分（2022年度, 47都道府県）\nr = {r_val:.3f} (p{p_str})',
    fontsize=13, fontweight='bold',
)
ax4.grid(True, alpha=0.2)
ax4.text(
    0.05, 0.95, f'相関係数 r = {r_val:.3f}',
    transform=ax4.transAxes, fontsize=11, va='top',
    bbox=dict(boxstyle='round', facecolor='#E3F2FD', alpha=0.8),
)

# 凡例
from matplotlib.patches import Patch
legend_elements = [
    Patch(facecolor='#E53935', label='SVM複合指標 上位5'),
    Patch(facecolor='#1E88E5', label='SVM複合指標 下位5'),
    Patch(facecolor='#BDBDBD', label='その他'),
]
ax4.legend(handles=legend_elements, fontsize=9, loc='lower right')

plt.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, '2024_U5_2_fig4_scatter.png'),
             bbox_inches='tight', dpi=150)
plt.close(fig4)
print("  → 2024_U5_2_fig4_scatter.png 保存完了")

print("\n" + "=" * 60)
print("全図の生成完了（4枚）")
print("=" * 60)
print("\n【主要知見】")
print(f"  SVM変数重み上位3変数: {weight_df['変数'].head(3).tolist()}")
print(f"  複合指標1位:  {rank_df.iloc[0]['都道府県']}（スコア = {rank_df.iloc[0]['複合指標（SVM）']:.3f}）")
print(f"  複合指標最下位: {rank_df.iloc[-1]['都道府県']}（スコア = {rank_df.iloc[-1]['複合指標（SVM）']:.3f}）")
print(f"  SVM指標 vs PCA指標 相関: r = {r_val:.3f}")
print(f"  PCA説明率: PC1={pca.explained_variance_ratio_[0]*100:.1f}%, "
      f"PC2={pca.explained_variance_ratio_[1]*100:.1f}%")

import numpy as np
import pandas as pd
from sklearn.preprocessing import StandardScaler
from scipy.cluster.hierarchy import linkage, fcluster

# SSDSE-A（市区町村データ、原論文と同じ粒度）で複合指標を作る
a = pd.read_csv('data/raw/SSDSE-A-2025.csv', encoding='cp932', header=2)
a = a[a['地域コード'].str.match(r'^R\d+', na=False)].copy()
a = a[pd.to_numeric(a['総人口'], errors='coerce') > 0]
num = lambda c: pd.to_numeric(a[c], errors='coerce')

# サブ指標（原論文の「主観に頼らない変数選択 → 重み付け」に対応する素材）
X = pd.DataFrame(index=(a['都道府県'] + a['市区町村']).values)
X['社会増減率'] = ((num('転入者数（日本人移動者）') - num('転出者数（日本人移動者）')) / num('総人口') * 100).values
X['自然増減率'] = ((num('出生数') - num('死亡数')) / num('総人口') * 100).values
X['若年度'] = (num('15歳未満人口') / num('総人口') * 100).values
X['財政力'] = (num('地方税（市町村財政）') / num('歳入決算総額（市町村財政）') * 100).values
X['雇用集積'] = (num('従業者数（民営）') / num('総人口') * 100).values
X['医療環境'] = (num('医師数') / num('総人口') * 10000).values
X = X.dropna()
print(f'市区町村数: {len(X)}')

# 階層クラスター分析（Ward法）でサブ指標を「要因」に分解（原論文の中核工程）
# → 変数側をクラスタリングして相関の強い指標群をまとめる
Xz = pd.DataFrame(StandardScaler().fit_transform(X), columns=X.columns, index=X.index)
corr_dist = 1 - Xz.corr().abs()          # 変数間の距離 = 1 − |相関|
Z = linkage(corr_dist.values[np.triu_indices(len(X.columns), 1)], method='ward')
groups = fcluster(Z, t=3, criterion='maxclust')
print('\n=== 変数の階層クラスタリング（要因分解） ===')
for g in sorted(set(groups)):
    print(f'  要因{g}: ' + '、'.join(X.columns[groups == g]))

# 要因ごとに平均して複合指標化 → 総合ランキング（原論文の複合指標の考え方）
comp = pd.DataFrame(index=X.index)
for g in sorted(set(groups)):
    comp[f'要因{g}'] = Xz.loc[:, X.columns[groups == g]].mean(axis=1)
comp['総合指標'] = comp.mean(axis=1)
print('\n=== 複合指標 総合ランキング ===')
print('上位5:', '、'.join(comp['総合指標'].nlargest(5).index))
print('下位5:', '、'.join(comp['総合指標'].nsmallest(5).index))