#!/usr/bin/env python3
# coding: utf-8
# ============================================================
# 【データの準備】実行前に以下のデータファイルを用意してください
#
#   必要ファイル:
#     ・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/2022_U1_daijin.py
# ============================================================

"""
==========================================================================
2022年度 統計データ分析コンペティション 総務大臣賞（大学生・一般の部）
「パンデミックは人流をどう変えたか ―地域の特性別に―」

【このスクリプトの目的】
論文の手法をトレースすることで、次の分析手法を学ぶ。
  - 時系列分析（延べ宿泊者数の2012〜2023年推移）
  - 地域別（6区分）比較
  - 重回帰分析（標準化係数）
  - 変化率と説明変数の散布図

【データ】
  data/raw/SSDSE-B-2026.csv（実データ）
  出典: 情報・システム研究機構「教育用標準データセット SSDSE-B」
  目的変数プロキシ: 延べ宿泊者数（G7101）
  比較期間:
    コロナ前 (2018-2019) vs コロナ禍 (2020-2021) vs 回復期 (2022-2023)

【必要ライブラリ】
  pip install pandas numpy matplotlib scipy statsmodels scikit-learn
==========================================================================
"""

import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
from matplotlib.gridspec import GridSpec
import statsmodels.api as sm
from scipy import stats
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import StandardScaler
import warnings
import os

warnings.filterwarnings('ignore')

# 日本語フォント設定（macOS）
plt.rcParams['font.family'] = 'Hiragino Sans'
plt.rcParams['axes.unicode_minus'] = False
plt.rcParams['font.size'] = 10

# 出力先ディレクトリ
OUT_DIR  = 'html/figures'
DATA_B   = 'data/raw/SSDSE-B-2026.csv'
os.makedirs(OUT_DIR, exist_ok=True)

def save_fig(name):
    path = os.path.join(OUT_DIR, f'2022_U1_{name}.png')
    plt.savefig(path, dpi=150, bbox_inches='tight', facecolor='white')
    plt.close()
    print(f'  → 保存: {path}')

# ==========================================================================
# 1. データ読み込み・前処理
# ==========================================================================
print("=" * 60)
print("【Step 1】SSDSE-B-2026 読み込み・前処理")
print("  出典: SSDSE-B-2026.csv（実データ）")
print("=" * 60)

df_b = pd.read_csv(DATA_B, encoding='cp932', header=1)
df_b = df_b[df_b['地域コード'].str.match(r'^R\d{5}$', na=False)].copy()

# 数値変換（エラーは NaN）
numeric_cols = [
    '延べ宿泊者数', '総人口', '年平均気温', '合計特殊出生率',
    '標準価格（平均価格）（住宅地）', '消費支出（二人以上の世帯）',
    '65歳以上人口', '月間有効求人数（一般）', '月間有効求職者数（一般）',
]
for col in numeric_cols:
    if col in df_b.columns:
        df_b[col] = pd.to_numeric(df_b[col], errors='coerce')

# 65歳以上人口 → 高齢化率
df_b['高齢化率'] = df_b['65歳以上人口'] / df_b['総人口'] * 100

# 有効求人倍率
df_b['有効求人倍率'] = df_b['月間有効求人数（一般）'] / df_b['月間有効求職者数（一般）']

# 宿泊者数の人口あたり（千人あたり）
df_b['宿泊者数_千人対'] = df_b['延べ宿泊者数'] / df_b['総人口'] * 1000

df_b['年度'] = pd.to_numeric(df_b['年度'], errors='coerce').astype(int)

print(f"  データ形状: {df_b.shape}  ( {df_b['年度'].min()}〜{df_b['年度'].max()}年 × 47都道府県 )")
print(f"  延べ宿泊者数: 非欠損={df_b['延べ宿泊者数'].notna().sum()}, "
      f"最小={df_b['延べ宿泊者数'].min():.0f}, 最大={df_b['延べ宿泊者数'].max():.0f}")
print()

# ==========================================================================
# 2. 変化率の計算
# ==========================================================================
print("=" * 60)
print("【Step 2】期間別 延べ宿泊者数の変化率を計算")
print("  コロナ前平均 (2018-2019) → コロナ禍 (2020-2021) → 回復期 (2022-2023)")
print("=" * 60)

def period_mean(df, years, col):
    return df[df['年度'].isin(years)].groupby('都道府県')[col].mean()

pre_mean  = period_mean(df_b, [2018, 2019], '延べ宿泊者数')
covid_mean = period_mean(df_b, [2020, 2021], '延べ宿泊者数')
rec_mean  = period_mean(df_b, [2022, 2023], '延べ宿泊者数')

# COVID期の変化率（対前期比%）
covid_chg = ((covid_mean - pre_mean) / pre_mean * 100).rename('covid_chg')
rec_chg   = ((rec_mean   - pre_mean) / pre_mean * 100).rename('rec_chg')

print(f"  全国平均 コロナ禍変化率: {covid_chg.mean():.1f}%")
print(f"  全国平均 回復期変化率:   {rec_chg.mean():.1f}%")
print()

# ==========================================================================
# 3. 地域区分の設定（8ブロック→6区分）
# ==========================================================================
REGION_MAP = {
    '北海道': '北海道・東北', '青森県': '北海道・東北', '岩手県': '北海道・東北',
    '宮城県': '北海道・東北', '秋田県': '北海道・東北', '山形県': '北海道・東北',
    '福島県': '北海道・東北',
    '茨城県': '関東', '栃木県': '関東', '群馬県': '関東', '埼玉県': '関東',
    '千葉県': '関東', '東京都': '関東', '神奈川県': '関東',
    '新潟県': '中部', '富山県': '中部', '石川県': '中部', '福井県': '中部',
    '山梨県': '中部', '長野県': '中部', '岐阜県': '中部', '静岡県': '中部',
    '愛知県': '中部',
    '三重県': '近畿', '滋賀県': '近畿', '京都府': '近畿', '大阪府': '近畿',
    '兵庫県': '近畿', '奈良県': '近畿', '和歌山県': '近畿',
    '鳥取県': '中国・四国', '島根県': '中国・四国', '岡山県': '中国・四国',
    '広島県': '中国・四国', '山口県': '中国・四国', '徳島県': '中国・四国',
    '香川県': '中国・四国', '愛媛県': '中国・四国', '高知県': '中国・四国',
    '福岡県': '九州・沖縄', '佐賀県': '九州・沖縄', '長崎県': '九州・沖縄',
    '熊本県': '九州・沖縄', '大分県': '九州・沖縄', '宮崎県': '九州・沖縄',
    '鹿児島県': '九州・沖縄', '沖縄県': '九州・沖縄',
}
REGION_ORDER = ['北海道・東北', '関東', '中部', '近畿', '中国・四国', '九州・沖縄']
REGION_COLORS = {
    '北海道・東北': '#1565C0',
    '関東':         '#E53935',
    '中部':         '#F57F17',
    '近畿':         '#2E7D32',
    '中国・四国':   '#6A1B9A',
    '九州・沖縄':   '#00838F',
}

df_b['地域区分'] = df_b['都道府県'].map(REGION_MAP)

# ==========================================================================
# 4. 図1：全国平均 延べ宿泊者数の時系列（2012〜2023年）
# ==========================================================================
print("=" * 60)
print("【Step 3】図の生成")
print("=" * 60)
print("  図1：全国平均 延べ宿泊者数の時系列（2012〜2023年）")

national_ts = df_b.groupby('年度')['延べ宿泊者数'].mean() / 1e6  # 百万人単位

fig, ax = plt.subplots(figsize=(11, 5.5))

# 期間背景色
ax.axvspan(2017.5, 2019.5, color='#E8F5E9', alpha=0.7, label='コロナ前 (2018-2019)')
ax.axvspan(2019.5, 2021.5, color='#FFEBEE', alpha=0.7, label='コロナ禍 (2020-2021)')
ax.axvspan(2021.5, 2023.5, color='#E3F2FD', alpha=0.7, label='回復期 (2022-2023)')

# 都道府県別の薄い折れ線
for pref in df_b['都道府県'].unique():
    pref_ts = df_b[df_b['都道府県'] == pref].set_index('年度')['延べ宿泊者数'] / 1e6
    pref_ts = pref_ts.sort_index()
    region = REGION_MAP.get(pref, '関東')
    ax.plot(pref_ts.index, pref_ts.values,
            color=REGION_COLORS.get(region, '#999'), alpha=0.15, linewidth=0.8)

# 全国平均（太線）
ax.plot(national_ts.index, national_ts.values,
        color='black', linewidth=2.5, marker='o', markersize=5,
        label='全国平均', zorder=5)

# 注釈
ax.annotate('COVID-19\n感染拡大',
            xy=(2020, national_ts[2020]), xytext=(2020.3, national_ts[2020] * 1.25),
            fontsize=9, color='#C62828',
            arrowprops=dict(arrowstyle='->', color='#C62828', lw=1.2))

ax.set_xlabel('年度', fontsize=12)
ax.set_ylabel('延べ宿泊者数（百万人）', fontsize=12)
ax.set_title('図1：全国平均 延べ宿泊者数の時系列推移（2012〜2023年）\n'
             '出典: SSDSE-B-2026（実データ）  各細線=都道府県別、太線=全国平均', fontsize=11)
ax.legend(loc='upper left', fontsize=9)
ax.set_xlim(2011.5, 2023.8)
ax.grid(True, alpha=0.3)
plt.tight_layout()
save_fig('fig1_timeseries')

# ==========================================================================
# 5. 図2：地域別（6区分）の変化率比較
# ==========================================================================
print("  図2：地域別（6区分）宿泊者数の変化率比較")

# COVID変化率のDataFrame
chg_df = pd.DataFrame({'covid_chg': covid_chg, 'rec_chg': rec_chg}).reset_index()
chg_df['地域区分'] = chg_df['都道府県'].map(REGION_MAP)
chg_df = chg_df.dropna(subset=['地域区分', 'covid_chg', 'rec_chg'])

region_covid = chg_df.groupby('地域区分')['covid_chg'].mean()
region_rec   = chg_df.groupby('地域区分')['rec_chg'].mean()

fig, axes = plt.subplots(1, 2, figsize=(13, 5.5), sharey=True)

for ax, (col_data, col_name, period_label, bg_color) in zip(
        axes,
        [(region_covid, 'covid_chg', 'コロナ禍 (2020-2021) vs 前期', '#FFEBEE'),
         (region_rec,   'rec_chg',   '回復期 (2022-2023) vs 前期',   '#E3F2FD')]):

    y_pos = np.arange(len(REGION_ORDER))
    vals  = [col_data.get(r, np.nan) for r in REGION_ORDER]
    colors = [REGION_COLORS[r] for r in REGION_ORDER]
    bars = ax.barh(y_pos, vals, color=colors, alpha=0.75, height=0.6)

    # 値ラベル
    for i, v in enumerate(vals):
        if not np.isnan(v):
            ax.text(v + (1 if v >= 0 else -1), i,
                    f'{v:.1f}%', va='center', fontsize=9,
                    ha='left' if v >= 0 else 'right')

    ax.axvline(0, color='black', linewidth=0.9, linestyle='--')
    ax.set_yticks(y_pos)
    ax.set_yticklabels(REGION_ORDER, fontsize=10)
    ax.set_xlabel('変化率（%）', fontsize=10)
    ax.set_title(period_label, fontsize=10)
    ax.set_facecolor(bg_color)
    ax.grid(True, axis='x', alpha=0.3)
    ax.invert_yaxis()

axes[0].set_title('コロナ禍 (2020-2021)\nvs コロナ前 (2018-2019)', fontsize=10)
axes[1].set_title('回復期 (2022-2023)\nvs コロナ前 (2018-2019)', fontsize=10)

plt.suptitle('図2：地域別（6区分）延べ宿泊者数の変化率\n'
             '出典: SSDSE-B-2026（実データ）  各区分は都道府県平均', fontsize=11, y=1.01)
plt.tight_layout()
save_fig('fig2_region')

# ==========================================================================
# 6. 重回帰分析：COVID変化率の地域間差を説明変数で解析
# ==========================================================================
print("  図3：重回帰 標準化係数プロット（コロナ禍変化率の規定要因）")

# 説明変数：コロナ前（2018-2019）平均
def pre_mean_var(df, years, col):
    return df[df['年度'].isin(years)].groupby('都道府県')[col].mean()

x_vars_raw = {
    '年平均気温':     pre_mean_var(df_b, [2018, 2019], '年平均気温'),
    '合計特殊出生率': pre_mean_var(df_b, [2018, 2019], '合計特殊出生率'),
    '住宅地価':       pre_mean_var(df_b, [2018, 2019], '標準価格（平均価格）（住宅地）'),
    '消費支出':       pre_mean_var(df_b, [2018, 2019], '消費支出（二人以上の世帯）'),
    '高齢化率':       pre_mean_var(df_b, [2018, 2019], '高齢化率'),
    '有効求人倍率':   pre_mean_var(df_b, [2018, 2019], '有効求人倍率'),
}

# 旅館施設数（観光地ダミーの代理変数）
x_vars_raw['旅館施設数'] = pre_mean_var(df_b, [2018, 2019], '旅館営業施設数（ホテルを含む）') \
    if '旅館営業施設数（ホテルを含む）' in df_b.columns else \
    pre_mean_var(df_b, [2018, 2019], '旅館営業施設客室数（ホテルを含む）')

x_df = pd.DataFrame(x_vars_raw)
reg_df = x_df.join(covid_chg).dropna()

print(f"  回帰サンプル数: N={len(reg_df)}")

X_raw = reg_df[list(x_vars_raw.keys())]
y_reg = reg_df['covid_chg']

# 標準化
scaler = StandardScaler()
X_std = scaler.fit_transform(X_raw)
X_std_df = pd.DataFrame(X_std, columns=X_raw.columns, index=X_raw.index)

# OLS（statsmodels）
X_sm = sm.add_constant(X_std_df)
ols = sm.OLS(y_reg, X_sm).fit()

print()
print(f"  R² = {ols.rsquared:.3f}  adj-R² = {ols.rsquared_adj:.3f}")
print(f"  {'変数名':<15} {'係数':>8} {'p値':>8} {'有意'}")
print("  " + "─" * 42)
for var in X_raw.columns:
    b = ols.params[var]; p = ols.pvalues[var]
    sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else ''
    print(f"  {var:<15} {b:>8.4f} {p:>8.4f} {sig}")
print()

# 係数プロット
var_labels = {
    '年平均気温':     '年平均気温（℃）',
    '合計特殊出生率': '合計特殊出生率',
    '住宅地価':       '住宅地標準価格',
    '消費支出':       '消費支出（世帯）',
    '高齢化率':       '高齢化率（%）',
    '有効求人倍率':   '有効求人倍率',
    '旅館施設数':     '旅館・ホテル施設数',
}

coefs  = [ols.params[v] for v in X_raw.columns]
pvals  = [ols.pvalues[v] for v in X_raw.columns]
ci_lo  = [ols.conf_int().loc[v, 0] for v in X_raw.columns]
ci_hi  = [ols.conf_int().loc[v, 1] for v in X_raw.columns]
labs   = [var_labels.get(v, v) for v in X_raw.columns]

fig, ax = plt.subplots(figsize=(10, 6))
y_pos = np.arange(len(coefs))
bar_colors = ['#C62828' if b < 0 else '#1565C0' for b, p in zip(coefs, pvals)]
bar_alpha  = [0.85 if p < 0.05 else 0.35 for p in pvals]

for i, (b, cl, ch, bc, ba) in enumerate(zip(coefs, ci_lo, ci_hi, bar_colors, bar_alpha)):
    ax.barh(i, b, color=bc, alpha=ba, height=0.6)
    ax.plot([cl, ch], [i, i], color='black', linewidth=1.5)
    ax.plot([cl, cl], [i - 0.15, i + 0.15], color='black', linewidth=1.5)
    ax.plot([ch, ch], [i - 0.15, i + 0.15], color='black', linewidth=1.5)

ax.axvline(0, color='black', linewidth=1.0, linestyle='--')
ax.set_yticks(y_pos)
ax.set_yticklabels(labs, fontsize=10)
ax.set_xlabel('標準化偏回帰係数（95% CI）\n※ 濃色: p<0.05 有意、薄色: 非有意', fontsize=10)
ax.set_title(f'図3：コロナ禍 延べ宿泊者数変化率の規定要因（重回帰 標準化係数）\n'
             f'N=47都道府県, R²={ols.rsquared:.3f}  出典: SSDSE-B-2026（実データ）', fontsize=11)
ax.grid(True, axis='x', alpha=0.3)

# 有意性マーク
for i, (b, ch, p) in enumerate(zip(coefs, ci_hi, pvals)):
    sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else ''
    if sig:
        ax.text(ch + 0.02, i, sig, va='center', fontsize=10,
                color='#C62828', fontweight='bold')

sig_patch   = mpatches.Patch(color='#1565C0', alpha=0.85, label='正の効果（p<0.05）')
nosig_patch = mpatches.Patch(color='#1565C0', alpha=0.35, label='p≥0.05（非有意）')
neg_patch   = mpatches.Patch(color='#C62828', alpha=0.85, label='負の効果（p<0.05）')
ax.legend(handles=[sig_patch, neg_patch, nosig_patch], fontsize=9, loc='lower right')

plt.tight_layout()
save_fig('fig3_coef')

# ==========================================================================
# 7. 図4：変化率 vs 主要説明変数の散布図
# ==========================================================================
print("  図4：COVID変化率 vs 主要説明変数 散布図")

# p値の低い順で上位4変数を選択
pval_sorted = sorted(zip(pvals, X_raw.columns), key=lambda x: x[0])
top4 = [v for _, v in pval_sorted[:4]]

fig, axes = plt.subplots(2, 2, figsize=(12, 10))
axes = axes.flatten()

scatter_colors = [REGION_COLORS.get(REGION_MAP.get(p, '関東'), '#666')
                  for p in reg_df.index]

for ax, var in zip(axes, top4):
    x_vals = reg_df[var].values
    y_vals = reg_df['covid_chg'].values
    N_pts  = len(x_vals)

    # jitter: np.linspaceベースの擬似ジッター（重なり軽減）
    jitter = np.linspace(-0.15, 0.15, N_pts)

    ax.scatter(x_vals + jitter * (x_vals.max() - x_vals.min()) * 0.03,
               y_vals, c=scatter_colors, alpha=0.75, s=60, edgecolors='white', linewidths=0.5)

    # 回帰直線
    slope, intercept, r_val, p_val, _ = stats.linregress(x_vals, y_vals)
    x_line = np.linspace(x_vals.min(), x_vals.max(), 100)
    ax.plot(x_line, intercept + slope * x_line,
            color='black', linewidth=1.5, linestyle='--',
            label=f'r={r_val:.3f}, p={p_val:.3f}')

    # 都道府県名ラベル（上位・下位3点）
    sort_idx = np.argsort(y_vals)
    label_idx = np.concatenate([sort_idx[:3], sort_idx[-3:]])
    for idx in label_idx:
        pref_name = reg_df.index[idx]
        ax.annotate(pref_name, (x_vals[idx], y_vals[idx]),
                    xytext=(4, 2), textcoords='offset points', fontsize=7, alpha=0.8)

    ax.axhline(0, color='gray', linewidth=0.8, linestyle=':')
    ax.set_xlabel(var_labels.get(var, var), fontsize=10)
    ax.set_ylabel('コロナ禍変化率（%）', fontsize=9)
    ax.set_title(f'{var_labels.get(var, var)} vs 変化率', fontsize=10)
    ax.legend(fontsize=9)
    ax.grid(True, alpha=0.25)

# 地域カラー凡例
legend_patches = [mpatches.Patch(color=c, label=r, alpha=0.75)
                  for r, c in REGION_COLORS.items()]
fig.legend(handles=legend_patches, loc='lower center',
           ncol=3, fontsize=9, bbox_to_anchor=(0.5, -0.02))

plt.suptitle('図4：コロナ禍 延べ宿泊者数変化率 vs 主要説明変数（上位4変数）\n'
             '出典: SSDSE-B-2026（実データ）  色=地域区分', fontsize=11, y=1.01)
plt.tight_layout()
save_fig('fig4_scatter')

# ==========================================================================
# 8. 記述統計サマリー
# ==========================================================================
print()
print("=" * 60)
print("【Step 4】記述統計サマリー")
print("=" * 60)

desc_cols = ['covid_chg', 'rec_chg']
desc_df = pd.DataFrame({'covid_chg': covid_chg, 'rec_chg': rec_chg})
print(desc_df.describe().round(2))
print()

# 地域別サマリー
region_summary = chg_df.groupby('地域区分')[['covid_chg', 'rec_chg']].mean().round(1)
region_summary.columns = ['コロナ禍変化率(%)', '回復期変化率(%)']
print("  ─ 地域別平均変化率 ─")
for r in REGION_ORDER:
    if r in region_summary.index:
        row = region_summary.loc[r]
        print(f"  {r:<12}  コロナ禍: {row.iloc[0]:>7.1f}%  回復期: {row.iloc[1]:>7.1f}%")
print()

# ==========================================================================
# 学習ポイント
# ==========================================================================
print("=" * 60)
print("【学習ポイント】この論文から学べるデータサイエンスの技術")
print("=" * 60)
print("""
1. 時系列分析の基本
   ・単純な集計（年度別平均）から変化のパターンを読み取る
   ・「外れ値」（コロナ禍の急落）の視覚的把握が重要

2. 期間比較と変化率
   ・「コロナ前 → コロナ禍 → 回復期」を定量化
   ・変化率 = (後期 - 前期) / 前期 × 100

3. 地域分類による集計
   ・47都道府県を6区分に分類し、地域差を浮き彫りにする
   ・観光地（旅館施設数が多い地域）ほど落ち込みが大きい傾向

4. 重回帰分析による規定要因の特定
   ・目的変数: COVID変化率（宿泊者数の落ち込み度）
   ・説明変数: 地域の人口・経済・観光・気候特性
   ・標準化係数で変数間の効果の大きさを比較

5. 散布図による個別確認
   ・各都道府県のポジションを確認（外れ値・傾向の確認）
   ・回帰直線と相関係数rでの変数関係の要約

6. 実データの活用
   ・SSDSE-B：47都道府県×12年の多変数データ（全て実測値）
   ・宿泊者数は人流（mobility）の代理変数として活用
""")

print("=" * 60)
print("すべての図を保存しました。")
print("  2022_U1_fig1_timeseries.png : 全国平均の時系列（2012-2023）")
print("  2022_U1_fig2_region.png     : 地域別6区分の変化率比較")
print("  2022_U1_fig3_coef.png       : 重回帰 標準化係数プロット")
print("  2022_U1_fig4_scatter.png    : 変化率 vs 主要説明変数の散布図")
print("=" * 60)