"""
2024_U4_katsuyo.py
金融資産購入経験の要因分析（代替: 大学進学率の決定要因）
統計活用奨励賞 [大学生・一般の部]
NGUYEN THI NGOC ANH、NGUYEN THI MINH QUY（青森中央学院大学）

実データ（SSDSE-B-2026, SSDSE-E-2026）による教育用再現コード
※個人レベルデータ非公開のため、都道府県別大学進学率（中央値超過=1）で代替
"""

# ============================================================
# 【データの準備】実行前に以下のデータファイルを用意してください
#
#   必要ファイル:
#     ・SSDSE-B-2026.csv
#       → data/raw/SSDSE-B-2026.csv に配置
#     ・SSDSE-E-2026.csv
#       → data/raw/SSDSE-E-2026.csv に配置
#
#   ダウンロード先:
#     https://www.nstac.go.jp/use/literacy/ssdse/
#     （SSDSE-B（社会・人口統計体系 都道府県データ） の CSV をダウンロード）
#     （SSDSE-E（社会・人口統計体系 都道府県の指標2） の CSV をダウンロード）
#
#   フォルダ配置（プロジェクトルートからの相対パス）:
#     code/                ← このスクリプトの場所
#     data/raw/            ← CSV ファイルをここに配置
#     html/figures/        ← 図の出力先（自動生成）
#
#   実行方法（ファイルを一切編集せず実行可能）:
#     python3 code/2024_U4_katsuyo.py
# ============================================================


import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import numpy as np
import pandas as pd
from scipy import stats
from numpy.linalg import lstsq

plt.rcParams['font.family'] = 'Hiragino Sans'
plt.rcParams['axes.unicode_minus'] = False

import os
FIGDIR = os.path.normpath('html/figures') + os.sep
DATA_B = 'data/raw/SSDSE-B-2026.csv'
DATA_E = 'data/raw/SSDSE-E-2026.csv'
os.makedirs(FIGDIR, exist_ok=True)

# ----------------------------------------------------------------
# データ読み込み: SSDSE-B 2022年断面
# ----------------------------------------------------------------
df_b = pd.read_csv(DATA_B, encoding='cp932', header=1)
mask_b = (df_b['地域コード'].str.match(r'^R\d{5}$', na=False) &
          (df_b['地域コード'] != 'R00000') &
          (df_b['年度'] == 2022))
df = df_b[mask_b].copy().reset_index(drop=True)

# SSDSE-E
df_e_raw = pd.read_csv(DATA_E, encoding='cp932', header=0)
df_e = df_e_raw.iloc[2:].copy()
df_e.columns = df_e_raw.iloc[1].values
df_e = df_e[df_e['都道府県'] != '全国'].reset_index(drop=True)

# ----------------------------------------------------------------
# 変数構築
# ----------------------------------------------------------------
df['高等学校卒業者数'] = pd.to_numeric(df['高等学校卒業者数'], errors='coerce')
df['高等学校卒業者のうち進学者数'] = pd.to_numeric(df['高等学校卒業者のうち進学者数'], errors='coerce')
df['進学率'] = df['高等学校卒業者のうち進学者数'] / df['高等学校卒業者数']

df['総人口'] = pd.to_numeric(df['総人口'], errors='coerce')
df['65歳以上人口'] = pd.to_numeric(df['65歳以上人口'], errors='coerce')
df['保育所等数'] = pd.to_numeric(df['保育所等数'], errors='coerce')
df['年平均気温'] = pd.to_numeric(df['年平均気温'], errors='coerce')

df['高齢化率'] = df['65歳以上人口'] / df['総人口']
df['保育所数千対'] = df['保育所等数'] / df['総人口'] * 10000   # 人口万人対

# SSDSE-E からの変数
def normalize_pref(s):
    return str(s).rstrip('県府都道').strip()

df['pref_short'] = df['都道府県'].apply(normalize_pref)
df_e['pref_short'] = df_e['都道府県'].apply(normalize_pref)

df_e_use = df_e.set_index('pref_short')
for col_e, col_new in [
    ('1人当たり県民所得（平成27年基準）', '1人当たり所得'),
    ('総面積（北方地域及び竹島を除く）', '総面積_ha'),
]:
    df[col_new] = df['pref_short'].map(
        pd.to_numeric(df_e_use[col_e], errors='coerce'))

df['人口密度'] = df['総人口'] / (df['総面積_ha'] / 100)   # 人/km²

# 目的変数: 進学率が全国中央値超 → 1
median_rate = df['進学率'].median()
df['進学率高（binary）'] = (df['進学率'] > median_rate).astype(int)

feature_names = ['1人当たり所得', '高齢化率', '人口密度', '保育所数千対', '年平均気温']
df_model = df[['都道府県', '進学率', '進学率高（binary）'] + feature_names].dropna()
df_model = df_model.reset_index(drop=True)

print(f"サンプル数: {len(df_model)}")
print(f"進学率 (全国): 中央値={median_rate:.3f}, "
      f"binary=1 が {df_model['進学率高（binary）'].sum()} 都道府県")

X_raw = df_model[feature_names].values.astype(float)
y = df_model['進学率高（binary）'].values.astype(float)
n = len(y)

# ================================================================
# ロジスティック回帰（勾配降下）
# ================================================================
def sigmoid(z):
    return 1 / (1 + np.exp(-np.clip(z, -500, 500)))

def logistic_fit(X, y_v, lr=0.1, n_iter=3000):
    n_s, k = X.shape
    b = np.zeros(k)
    for _ in range(n_iter):
        p = sigmoid(X @ b)
        grad = X.T @ (p - y_v) / n_s
        b -= lr * grad
    return b

X_m = X_raw.mean(axis=0)
X_s = X_raw.std(axis=0)
X_s[X_s == 0] = 1
X_std = (X_raw - X_m) / X_s
X_c = np.column_stack([np.ones(n), X_std])

coef = logistic_fit(X_c, y)

# 標準誤差（Fisher情報行列逆）
p_hat = sigmoid(X_c @ coef)
W = p_hat * (1 - p_hat)
W = np.maximum(W, 1e-6)
H = X_c.T @ (W[:, None] * X_c)
try:
    cov_mat = np.linalg.pinv(H)
except Exception:
    cov_mat = np.eye(len(coef)) * 0.01
se_coef = np.sqrt(np.abs(np.diag(cov_mat)))

z_stats = coef / (se_coef + 1e-12)
p_vals = 2 * (1 - stats.norm.cdf(np.abs(z_stats)))
or_vals = np.exp(coef)
ci_low = np.exp(coef - 1.96 * se_coef)
ci_high = np.exp(coef + 1.96 * se_coef)

# ================================================================
# 図1: 目的変数の分布（進学率の分布 + binary分け）
# ================================================================
fig, axes = plt.subplots(1, 2, figsize=(10, 5))

ax = axes[0]
counts = df_model['進学率高（binary）'].value_counts().sort_index()
labels_pie = ['進学率 低（0）', '進学率 高（1）']
colors_pie = ['#90A4AE', '#2E7D32']
wedges, texts, autotexts = ax.pie(
    counts.values, labels=labels_pie, colors=colors_pie,
    autopct='%1.1f%%', startangle=90,
    textprops={'fontsize': 11})
ax.set_title("大学進学率（中央値超=高）の割合", fontsize=12)

ax = axes[1]
groups_sorted = df_model.sort_values('進学率')
bars = ax.bar(range(len(groups_sorted)), groups_sorted['進学率'] * 100,
              color=['#C62828' if v == 1 else '#1565C0'
                     for v in groups_sorted['進学率高（binary）']],
              alpha=0.8)
ax.axhline(median_rate * 100, color='black', lw=2, linestyle='--',
           label=f'中央値 {median_rate*100:.1f}%')
ax.set_xlabel("都道府県（進学率昇順）")
ax.set_ylabel("大学進学率（%）")
ax.set_title("都道府県別大学進学率と中央値閾値", fontsize=12)
ax.legend(fontsize=10)
ax.grid(True, axis='y', alpha=0.3)

red_p = mpatches.Patch(color='#C62828', alpha=0.8, label='高（中央値超）')
blue_p = mpatches.Patch(color='#1565C0', alpha=0.8, label='低（中央値以下）')
axes[1].legend(handles=[red_p, blue_p, plt.Line2D([0], [0], color='black',
               lw=2, linestyle='--', label=f'中央値={median_rate*100:.1f}%')],
               fontsize=9)

fig.suptitle("図1: 目的変数（大学進学率 高/低）の分布", fontsize=13)
plt.tight_layout()
plt.savefig(FIGDIR + "2024_U4_fig1_dist.png", dpi=150)
plt.close()
print("fig1 saved")

# ================================================================
# 図2: VIF棒グラフ
# ================================================================
def calc_vif(X):
    n_s, k_s = X.shape
    vifs = []
    for j in range(k_s):
        Xj = X[:, j]
        Xi = np.delete(X, j, axis=1)
        Xi_c = np.column_stack([np.ones(n_s), Xi])
        b_j, _, _, _ = lstsq(Xi_c, Xj, rcond=None)
        pred_j = Xi_c @ b_j
        r2_j = 1 - np.var(Xj - pred_j) / max(np.var(Xj), 1e-10)
        vifs.append(1 / max(1 - r2_j, 1e-10))
    return vifs

vif_vals = calc_vif(X_raw)

fig, ax = plt.subplots(figsize=(8, 5))
colors_vif = ['#C62828' if v > 10 else ('#F9A825' if v > 5 else '#2E7D32') for v in vif_vals]
bars = ax.barh(feature_names, vif_vals, color=colors_vif, alpha=0.85)
ax.axvline(5, color='#F9A825', lw=1.5, linestyle='--', label='VIF=5（警戒水準）')
ax.axvline(10, color='#C62828', lw=1.5, linestyle='--', label='VIF=10（問題水準）')

for bar, v in zip(bars, vif_vals):
    ax.text(v + 0.05, bar.get_y() + bar.get_height() / 2,
            f'{v:.2f}', va='center', fontsize=10)

ax.set_xlabel("VIF（分散膨張係数）")
ax.set_title("図2: 説明変数のVIF（多重共線性チェック）", fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, axis='x', alpha=0.3)
plt.tight_layout()
plt.savefig(FIGDIR + "2024_U4_fig2_vif.png", dpi=150)
plt.close()
print("fig2 saved")

# ================================================================
# 図3: オッズ比 Forest Plot
# ================================================================
fig, ax = plt.subplots(figsize=(9, 6))

or_display = or_vals[1:]
ci_low_d = ci_low[1:]
ci_high_d = ci_high[1:]
p_display = p_vals[1:]

y_pos = np.arange(len(feature_names))
colors_or = ['#C62828' if or_v > 1 else '#1565C0' for or_v in or_display]

ax.axvline(1.0, color='black', lw=1.5, linestyle='-')

x_max = min(max(ci_high_d) * 1.1, 20)
for i, (or_v, cl, ch, p_v) in enumerate(zip(or_display, ci_low_d, ci_high_d, p_display)):
    cl_plot = max(cl, 0.05)
    ch_plot = min(ch, x_max - 0.5)
    ax.plot([cl_plot, ch_plot], [i, i], color=colors_or[i], lw=2.5, alpha=0.8)
    ax.scatter([or_v], [i], color=colors_or[i], s=80, zorder=5)
    sig = "***" if p_v < 0.001 else ("**" if p_v < 0.01 else ("*" if p_v < 0.05 else "n.s."))
    ax.text(min(ch_plot + 0.1, x_max - 0.2), i,
            f'OR={or_v:.2f} {sig}', va='center', fontsize=10)

ax.set_yticks(y_pos)
ax.set_yticklabels(feature_names, fontsize=11)
ax.set_xlabel("オッズ比（95% CI）")
ax.set_title("図3: ロジスティック回帰 オッズ比 Forest Plot\n（目的変数: 大学進学率 高/低）", fontsize=13)
ax.set_xlim(0.0, x_max)

red_patch = mpatches.Patch(color='#C62828', alpha=0.8, label='OR>1（進学率高を増加）')
blue_patch = mpatches.Patch(color='#1565C0', alpha=0.8, label='OR<1（進学率高を減少）')
ax.legend(handles=[red_patch, blue_patch], fontsize=10)
ax.grid(True, axis='x', alpha=0.3)
plt.tight_layout()
plt.savefig(FIGDIR + "2024_U4_fig3_odds.png", dpi=150)
plt.close()
print("fig3 saved")

# ================================================================
# 図4: 主要変数別の予測確率（線グラフ）
# ================================================================
fig, axes = plt.subplots(1, 3, figsize=(13, 5))

# 各変数の平均値（他変数固定用）
x_mean_std = np.zeros(len(feature_names))   # standardized mean = 0

def predict_prob(var_idx, var_range_raw, coef, X_m, X_s):
    """1変数を変化させたときの予測確率"""
    var_range_std = (var_range_raw - X_m[var_idx]) / X_s[var_idx]
    probs = []
    for v_std in var_range_std:
        x_pts = x_mean_std.copy()
        x_pts[var_idx] = v_std
        z = coef[0] + coef[1:] @ x_pts
        probs.append(sigmoid(z) * 100)
    return np.array(probs)

var_configs = [
    (0, '1人当たり所得', '1人当たり県民所得（万円）', None),
    (1, '高齢化率', '高齢化率', None),
    (3, '保育所数千対', '保育所数（人口万人あたり）', None),
]

for ax, (vi, vname, xlabel, _) in zip(axes, var_configs):
    col_data = X_raw[:, vi]
    var_range = np.linspace(col_data.min(), col_data.max(), 200)
    prob_line = predict_prob(vi, var_range, coef, X_m, X_s)

    ax.plot(var_range, prob_line, color='#1565C0', lw=2.5)
    ax.scatter(col_data, sigmoid(X_c @ coef) * 100, color='gray', alpha=0.5, s=20)
    ax.set_xlabel(xlabel)
    ax.set_ylabel("大学進学率高の予測確率（%）")
    ax.set_title(f"{vname}\nと進学率予測確率")
    ax.grid(True, alpha=0.3)
    ax.set_ylim(0, 100)

fig.suptitle("図4: 主要変数別 大学進学率高の予測確率", fontsize=13)
plt.tight_layout()
plt.savefig(FIGDIR + "2024_U4_fig4_prob.png", dpi=150)
plt.close()
print("fig4 saved")
print("All figures saved.")