"""
教育用再現コード: 2023年 統計データ分析コンペティション 審査員奨励賞 [高校生の部]
=================================================================
論文タイトル：日本人の英語能力の実態とその背景 ～諸外国と比較して～

【分析概要】
  データ：SSDSE-B-2026.csv（都道府県データ）、2022年度データ（47都道府県）

  ※ EF-EPI（国際英語能力指数）は公的統計データとして入手不可のため、
    都道府県別データで「英語力」に関連する指標を代理変数として分析。

  目的変数（英語力 proxy）：大学進学率
    = E4602 / E4601 × 100  （大学進学者数 / 卒業者数）
    ※ 大学入試には英語が必須であり、大学進学率は英語学力の代理指標

  説明変数：
    国際経験 proxy:  G5105 / A1101 × 10000  （国際旅券発行件数 per 1万人）
    経済格差:        L3221                   （消費支出, 千円）
                     C5401                   （住宅地価格, 千円/m²）
    教育投資:        L322108 / L3221 × 100   （教育費割合, %）
    英語環境:        G7102 / G7101 × 100     （外国人宿泊者割合, %）
    若年人口（男）:  A130201 / A1302 × 100   （15〜64歳男性人口比率, %）
    若年人口（女）:  A130202 / A1302 × 100   （15〜64歳女性人口比率, %）

  "国際比較" の代替：47都道府県を6地域に分けて地域間比較

  Step1. 散布図: 国際旅券発行率 vs 大学進学率
  Step2. 棒グラフ: 地域別（6グループ）大学進学率の比較
  Step3. 相関棒グラフ: 全変数 vs 大学進学率
  Step4. 散布図: 教育費割合 vs 大学進学率（回帰直線付き）

【データ出典】
  SSDSE-B-2026.csv: 社会・人口統計体系（都道府県データ）
  統計センターより公表の実データ

【データサイエンス学習ポイント】
  1. 記述統計（平均・標準偏差・最大最小）
  2. 相関分析（ピアソン相関係数と無相関検定）
  3. 地域間比較（一元配置分散分析的視点）
  4. 代理変数の概念（proxy variable）
=================================================================
"""

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


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

# ──────────────────────────────────────────────────────────────
# 共通設定
# ──────────────────────────────────────────────────────────────
plt.rcParams['font.family'] = 'Hiragino Sans'
plt.rcParams['axes.unicode_minus'] = False
plt.rcParams['figure.dpi'] = 150

import os
FIG_DIR = 'html/figures'
DATA_DIR = 'data/raw'
os.makedirs(FIG_DIR, exist_ok=True)

# ================================================================
# ■ 実データ読み込み（SSDSE-B-2026: 都道府県データ）
# ================================================================
df_raw = pd.read_csv(
    os.path.join(DATA_DIR, 'SSDSE-B-2026.csv'),
    header=0,        # 1行目がコード名（A1101 等）
    encoding='cp932'
)

# 都道府県コード（R01000〜R47000）の行のみ抽出
import re
mask_pref = df_raw['Code'].astype(str).str.match(r'^R\d{5}$')
df_pref = df_raw[mask_pref].copy()

# 年度列の名前を統一（先頭列 = 'SSDSE-B-2026' がそのまま年度値）
df_pref = df_pref.rename(columns={'SSDSE-B-2026': '年度'})
df_pref['年度'] = pd.to_numeric(df_pref['年度'], errors='coerce')

# 2022年度に絞る
df = df_pref[df_pref['年度'] == 2022].copy().reset_index(drop=True)

# 都道府県名の取得
pref_names = df['Prefecture'].astype(str).tolist()

print("=" * 60)
print(f"■ 対象: 2022年度 都道府県データ N={len(df)}")
print("=" * 60)

# ── 数値変換 ─────────────────────────────────────────────────
use_cols = ['E4601', 'E4602', 'G5105', 'A1101',
            'L3221', 'C5401', 'L322108',
            'G7101', 'G7102',
            'A130201', 'A130202', 'A1302']
for c in use_cols:
    df[c] = pd.to_numeric(df[c], errors='coerce')

# ── 指標の計算 ───────────────────────────────────────────────
# 目的変数: 大学進学率（英語力 proxy）
df['大学進学率']     = df['E4602'] / df['E4601'].clip(1) * 100

# 説明変数
df['国際旅券発行率'] = df['G5105'] / df['A1101'].clip(1) * 10000   # 件 / 万人
df['消費支出']       = df['L3221'] / 1000                           # 千円
df['住宅地価格']     = df['C5401']                                  # 千円/m²
df['教育費割合']     = df['L322108'] / df['L3221'].clip(1) * 100    # %
df['外国人宿泊割合'] = df['G7102'] / df['G7101'].clip(1) * 100      # %
df['若年人口比_男']  = df['A130201'] / df['A1302'].clip(1) * 100    # %
df['若年人口比_女']  = df['A130202'] / df['A1302'].clip(1) * 100    # %

TARGET = '大学進学率'
VARS = ['国際旅券発行率', '消費支出', '住宅地価格', '教育費割合',
        '外国人宿泊割合', '若年人口比_男', '若年人口比_女']

# 欠損確認・除外
needed = [TARGET] + VARS
df_clean = df.dropna(subset=needed).reset_index(drop=True)
pref_clean = [pref_names[i] for i in df_clean.index] if len(df_clean) == len(df) \
             else df_clean['Prefecture'].astype(str).tolist()
N = len(df_clean)
print(f"有効都道府県数（欠損除外後）: N={N}")

# ── 記述統計 ─────────────────────────────────────────────────
print("\n── 記述統計 ──")
print(df_clean[needed].describe().round(2))

# ── 地域グループ（6地域）─────────────────────────────────────
# 都道府県コード（R01000〜R47000）→ 番号で判定
df_clean['pref_no'] = df_clean['Code'].astype(str).str.extract(r'R(\d{2})').astype(int)

region_map = {
    '北海道・東北': [1, 2, 3, 4, 5, 6, 7],
    '関東':         [8, 9, 10, 11, 12, 13, 14],
    '中部':         [15, 16, 17, 18, 19, 20, 21, 22, 23],
    '近畿':         [24, 25, 26, 27, 28, 29, 30],
    '中国・四国':   [31, 32, 33, 34, 35, 36, 37, 38, 39],
    '九州・沖縄':   [40, 41, 42, 43, 44, 45, 46, 47],
}

def assign_region(pno):
    for region, nos in region_map.items():
        if pno in nos:
            return region
    return 'その他'

df_clean['地域'] = df_clean['pref_no'].apply(assign_region)

# ── 相関係数（全変数 vs 大学進学率）─────────────────────────
print("\n── 相関分析（vs 大学進学率）──")
print(f"  {'変数':<18} {'r':>8} {'p値':>10} {'有意':>6}")
print("  " + "-" * 46)
corrs = {}
pvals = {}
for v in VARS:
    x = df_clean[v].values
    y = df_clean[TARGET].values
    r, p = stats.pearsonr(x, y)
    corrs[v] = r
    pvals[v] = p
    sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else 'n.s.'
    print(f"  {v:<18} {r:>8.4f} {p:>10.4f} {sig:>6}")

# ================================================================
# ■ 図の生成（4枚）
# ================================================================

REGION_ORDER  = ['北海道・東北', '関東', '中部', '近畿', '中国・四国', '九州・沖縄']
REGION_COLORS = ['#1565C0', '#C62828', '#2E7D32', '#E65100', '#6A1B9A', '#00838F']

# ────────────────────────────────────────────────────────────────
# 図1: 散布図 — 国際旅券発行率 vs 大学進学率
# ────────────────────────────────────────────────────────────────
fig1, ax1 = plt.subplots(figsize=(9, 6))

x1 = df_clean['国際旅券発行率'].values
y1 = df_clean[TARGET].values

# 地域カラーで描画
for region, color in zip(REGION_ORDER, REGION_COLORS):
    mask_r = df_clean['地域'] == region
    ax1.scatter(x1[mask_r], y1[mask_r], color=color, s=60, alpha=0.8,
                label=region, zorder=3)

# 回帰直線
z1 = np.polyfit(x1, y1, 1)
xs1 = np.linspace(x1.min(), x1.max(), 200)
ax1.plot(xs1, np.poly1d(z1)(xs1), 'k--', linewidth=1.5, alpha=0.7, label='回帰直線')

r1, p1 = stats.pearsonr(x1, y1)
ax1.set_xlabel('国際旅券発行率（件 / 万人）\n＝国際経験の代理変数', fontsize=11)
ax1.set_ylabel('大学進学率（%）\n＝英語力の代理変数', fontsize=11)
ax1.set_title(f'国際経験と大学進学率の関係（都道府県別, 2022年）\nr = {r1:.3f}  p = {p1:.4f}',
              fontsize=12, fontweight='bold')
ax1.legend(fontsize=9, loc='lower right', ncol=2)
ax1.grid(True, alpha=0.3)

# 主要都道府県にラベル
highlight = {'東京都': 'top', '大阪府': 'bottom', '沖縄県': 'top', '鳥取県': 'bottom', '愛知県': 'top'}
for i, name in enumerate(df_clean['Prefecture'].tolist()):
    if name in highlight:
        va = highlight[name]
        ax1.annotate(name, (x1[i], y1[i]), textcoords='offset points',
                     xytext=(0, 6 if va == 'top' else -10),
                     ha='center', fontsize=8, color='#333')

plt.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2023_H5_1_fig1_scatter.png'), bbox_inches='tight', dpi=150)
plt.close(fig1)
print("\n図1保存: 2023_H5_1_fig1_scatter.png")

# ────────────────────────────────────────────────────────────────
# 図2: 棒グラフ — 地域別 大学進学率の比較（6地域グループ）
# ────────────────────────────────────────────────────────────────
fig2, ax2 = plt.subplots(figsize=(10, 6))

region_means = {}
region_stds  = {}
region_ns    = {}
for region in REGION_ORDER:
    vals = df_clean.loc[df_clean['地域'] == region, TARGET].values
    region_means[region] = vals.mean()
    region_stds[region]  = vals.std(ddof=1) if len(vals) > 1 else 0
    region_ns[region]    = len(vals)

means = [region_means[r] for r in REGION_ORDER]
stds  = [region_stds[r]  for r in REGION_ORDER]
ns    = [region_ns[r]    for r in REGION_ORDER]
x2pos = np.arange(len(REGION_ORDER))

bars = ax2.bar(x2pos, means, color=REGION_COLORS, alpha=0.80,
               edgecolor='white', linewidth=1.2, width=0.6)
ax2.errorbar(x2pos, means, yerr=stds, fmt='none',
             color='#333', capsize=5, linewidth=1.5)

# 全国平均
national_mean = df_clean[TARGET].mean()
ax2.axhline(national_mean, color='#555', linestyle='--', linewidth=1.5,
            label=f'全国平均 {national_mean:.1f}%')

for i, (m, s, n) in enumerate(zip(means, stds, ns)):
    ax2.text(i, m + s + 0.5, f'{m:.1f}%\n(N={n})', ha='center', va='bottom',
             fontsize=9, fontweight='bold')

ax2.set_xticks(x2pos)
ax2.set_xticklabels(REGION_ORDER, fontsize=10)
ax2.set_ylabel('大学進学率（%）', fontsize=11)
ax2.set_title('地域別 大学進学率の比較（2022年度 都道府県データ）\n（棒=平均, エラーバー=標準偏差）',
              fontsize=12, fontweight='bold')
ax2.legend(fontsize=10)
ax2.grid(axis='y', alpha=0.3)
ax2.set_ylim(0, max(means) + max(stds) + 8)

plt.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2023_H5_1_fig2_region.png'), bbox_inches='tight', dpi=150)
plt.close(fig2)
print("図2保存: 2023_H5_1_fig2_region.png")

# ────────────────────────────────────────────────────────────────
# 図3: 相関棒グラフ — 全変数 vs 大学進学率
# ────────────────────────────────────────────────────────────────
fig3, ax3 = plt.subplots(figsize=(9, 5))

var_labels = {
    '国際旅券発行率': '国際旅券発行率\n（国際経験 proxy）',
    '消費支出':       '消費支出\n（経済力）',
    '住宅地価格':     '住宅地価格\n（地価水準）',
    '教育費割合':     '教育費割合\n（教育投資）',
    '外国人宿泊割合': '外国人宿泊者割合\n（英語環境）',
    '若年人口比_男':  '若年人口比（男）\n（15〜64歳）',
    '若年人口比_女':  '若年人口比（女）\n（15〜64歳）',
}

r_vals  = [corrs[v] for v in VARS]
p_vals  = [pvals[v] for v in VARS]
labels3 = [var_labels[v] for v in VARS]
colors3 = []
for r, p in zip(r_vals, p_vals):
    if p < 0.05:
        colors3.append('#C62828' if r > 0 else '#1565C0')
    else:
        colors3.append('#9E9E9E')

y_pos3 = np.arange(len(VARS))
ax3.barh(y_pos3, r_vals, color=colors3, alpha=0.78, edgecolor='white', height=0.6)
ax3.axvline(0, color='gray', linestyle='--', linewidth=1.0)
ax3.set_yticks(y_pos3)
ax3.set_yticklabels(labels3, fontsize=9.5)
ax3.set_xlabel('ピアソン相関係数', fontsize=11)
ax3.set_title('各変数と大学進学率の相関係数（2022年度, N=47都道府県）\n'
              '（赤=正相関有意, 青=負相関有意, 灰=非有意 p≥0.05）',
              fontsize=11, fontweight='bold')
ax3.invert_yaxis()
ax3.grid(axis='x', alpha=0.3)

for i, (r, p) in enumerate(zip(r_vals, p_vals)):
    sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else ''
    offset = 0.015 if r >= 0 else -0.015
    ha = 'left' if r >= 0 else 'right'
    ax3.text(r + offset, i, f'{r:.3f}{sig}', va='center', ha=ha, fontsize=9)

from matplotlib.patches import Patch
legend_patches = [
    Patch(color='#C62828', alpha=0.78, label='正相関（p<0.05）'),
    Patch(color='#1565C0', alpha=0.78, label='負相関（p<0.05）'),
    Patch(color='#9E9E9E', alpha=0.78, label='非有意（p≥0.05）'),
]
ax3.legend(handles=legend_patches, fontsize=9, loc='lower right')

plt.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2023_H5_1_fig3_corr.png'), bbox_inches='tight', dpi=150)
plt.close(fig3)
print("図3保存: 2023_H5_1_fig3_corr.png")

# ────────────────────────────────────────────────────────────────
# 図4: 散布図 — 教育費割合 vs 大学進学率（回帰直線付き）
# ────────────────────────────────────────────────────────────────
fig4, ax4 = plt.subplots(figsize=(9, 6))

x4 = df_clean['教育費割合'].values
y4 = df_clean[TARGET].values

for region, color in zip(REGION_ORDER, REGION_COLORS):
    mask_r4 = df_clean['地域'] == region
    ax4.scatter(x4[mask_r4], y4[mask_r4], color=color, s=60, alpha=0.8,
                label=region, zorder=3)

# 回帰直線
z4 = np.polyfit(x4, y4, 1)
xs4 = np.linspace(x4.min(), x4.max(), 200)
ys4 = np.poly1d(z4)(xs4)
ax4.plot(xs4, ys4, 'k--', linewidth=1.8, alpha=0.8, label=f'回帰直線 y={z4[0]:.2f}x+{z4[1]:.1f}')

r4, p4 = stats.pearsonr(x4, y4)

# 信頼区間（bootstrap近似 → 簡易: SE from residuals）
n4 = len(x4)
y_pred4 = np.poly1d(z4)(x4)
resid4  = y4 - y_pred4
se4     = np.sqrt(np.sum(resid4**2) / (n4 - 2))
x_mean4 = x4.mean()
ss_x4   = np.sum((x4 - x_mean4)**2)
se_y4   = se4 * np.sqrt(1/n4 + (xs4 - x_mean4)**2 / ss_x4)
ax4.fill_between(xs4, ys4 - 1.96*se_y4, ys4 + 1.96*se_y4,
                 color='gray', alpha=0.15, label='95%信頼区間')

# 主要都道府県ラベル
highlight4 = {'東京都', '大阪府', '沖縄県', '秋田県', '京都府'}
for i, name in enumerate(df_clean['Prefecture'].tolist()):
    if name in highlight4:
        ax4.annotate(name, (x4[i], y4[i]), textcoords='offset points',
                     xytext=(4, 3), ha='left', fontsize=8, color='#333')

ax4.set_xlabel('教育費割合（消費支出に占める割合, %）\n＝教育投資の代理変数', fontsize=11)
ax4.set_ylabel('大学進学率（%）\n＝英語力の代理変数', fontsize=11)
ax4.set_title(f'教育投資と大学進学率の関係（都道府県別, 2022年）\n'
              f'r = {r4:.3f}  p = {p4:.4f}',
              fontsize=12, fontweight='bold')
ax4.legend(fontsize=9, loc='upper left')
ax4.grid(True, alpha=0.3)

plt.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, '2023_H5_1_fig4_educ.png'), bbox_inches='tight', dpi=150)
plt.close(fig4)
print("図4保存: 2023_H5_1_fig4_educ.png")

print("\n全図の生成完了（4枚）")
print(f"  fig1_scatter.png : 国際旅券発行率 vs 大学進学率（地域カラー散布図）")
print(f"  fig2_region.png  : 地域別 大学進学率の棒グラフ比較")
print(f"  fig3_corr.png    : 全変数 vs 大学進学率 相関棒グラフ")
print(f"  fig4_educ.png    : 教育費割合 vs 大学進学率（回帰直線・信頼区間）")