"""
教育用再現コード: 2024年 統計データ分析コンペティション 審査員奨励賞 [高校生の部]
=================================================================
論文タイトル：顎・足・枕が多様な観光客にどのような影響を及ぼすのか
著者：原田理矢（鳥取城北高等学校）

【分析概要】
  ※「顎足枕（あごあしまくら）」= 食事・交通・宿泊 の観光業用語
  データ：SSDSE-B-2026.csv（政府統計の総合窓口）、2022年度 都道府県データ
  目的変数1: 延べ宿泊者数 per capita（一人当たり宿泊需要 = 観光地としての総合的魅力の代理）
  目的変数2: 外国人延べ宿泊者数 / 延べ宿泊者数（外国人比率 = 国際的魅力の代理）

  説明変数:
    顎（食事）  : 食料費割合 = 食料費 / 消費支出
    足（交通）  : 交通・通信費割合 = 交通・通信費 / 消費支出
    枕（宿泊）  : 旅館数密度 = 旅館営業施設数 / 総人口 × 10000（10000人当たり）
    施設規模    : 客室数 / 施設数（1施設当たり客室数）
    気温        : 年平均気温
    高齢化率    : 65歳以上人口 / 総人口

  Step1. 6変数の相関行列ヒートマップ
  Step2. 旅館数密度 vs 一人当たり宿泊者数の散布図（都道府県ラベル付き）
  Step3. OLS標準化回帰係数の棒グラフ（2モデル比較）
  Step4. 一人当たり宿泊者数 Top10都道府県の棒グラフ

【データサイエンス学習ポイント】
  1. 実公的データ（SSDSE-B）の読み込みと前処理
  2. 派生変数（比率・密度）の作成
  3. 重回帰分析・標準化係数による変数重要度の比較
  4. 2つのアウトカム（宿泊需要 vs 外国人比率）での結果の違いを考察
=================================================================
"""

# ============================================================
# 【データの準備】実行前に以下のデータファイルを用意してください
#
#   必要ファイル:
#     ・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/2024_H5_6_shorei.py
# ============================================================


import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import numpy as np
from scipy import stats
import statsmodels.api as sm
from sklearn.preprocessing import StandardScaler
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'
os.makedirs(FIG_DIR, exist_ok=True)

DATA_PATH = 'data/raw/SSDSE-B-2026.csv'

# ================================================================
# ■ データ読み込み（SSDSE-B-2026）
# ================================================================
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("=" * 60)
print(f"■ SSDSE-B-2026 読み込み完了: {len(df)}都道府県（2022年度）")
print("=" * 60)

# ================================================================
# ■ 派生変数の作成
# ================================================================
# 元列をnumericに変換
num_cols = [
    '総人口', '65歳以上人口',
    '延べ宿泊者数', '外国人延べ宿泊者数',
    '旅館営業施設数（ホテルを含む）', '旅館営業施設客室数（ホテルを含む）',
    '食料費（二人以上の世帯）', '消費支出（二人以上の世帯）',
    '交通・通信費（二人以上の世帯）',
    '年平均気温',
]
for c in num_cols:
    df[c] = pd.to_numeric(df[c], errors='coerce')

# 目的変数
df['宿泊者数_per_capita'] = df['延べ宿泊者数'] / df['総人口']          # 一人当たり宿泊需要
df['外国人比率']          = df['外国人延べ宿泊者数'] / df['延べ宿泊者数']  # 外国人シェア

# 説明変数
df['食料費割合']    = df['食料費（二人以上の世帯）'] / df['消費支出（二人以上の世帯）']   # 顎（食）
df['交通費割合']    = df['交通・通信費（二人以上の世帯）'] / df['消費支出（二人以上の世帯）']  # 足（交通）
df['旅館数密度']    = df['旅館営業施設数（ホテルを含む）'] / df['総人口'] * 10000          # 枕（宿泊インフラ）
df['客室数_per_施設'] = df['旅館営業施設客室数（ホテルを含む）'] / df['旅館営業施設数（ホテルを含む）']  # 施設規模
df['気温']          = df['年平均気温']
df['高齢化率']      = df['65歳以上人口'] / df['総人口']

# 分析列を揃えてNaN行を削除
FEATURES = ['食料費割合', '交通費割合', '旅館数密度', '客室数_per_施設', '気温', '高齢化率']
OUTCOMES = ['宿泊者数_per_capita', '外国人比率']
df_model = df[['都道府県'] + FEATURES + OUTCOMES].dropna().reset_index(drop=True)

print(f"■ 欠損除去後: {len(df_model)}都道府県\n")
print(df_model[FEATURES + OUTCOMES].describe().round(4))

# ================================================================
# ■ 重回帰分析（statsmodels OLS）
# ================================================================
def run_ols(y_col, df_m, feature_cols, label):
    X = sm.add_constant(df_m[feature_cols].values)
    y = df_m[y_col].values
    model = sm.OLS(y, X).fit()
    print(f"\n{'='*60}")
    print(f"■ OLS: 目的変数 = {label} ({y_col})")
    print(model.summary(xname=['const'] + feature_cols))
    return model

model1 = run_ols('宿泊者数_per_capita', df_model, FEATURES, '一人当たり延べ宿泊者数')
model2 = run_ols('外国人比率',          df_model, FEATURES, '外国人延べ宿泊者数比率')

# 標準化係数
scaler = StandardScaler()
X_std = scaler.fit_transform(df_model[FEATURES].values)

def std_coefs(y_col, X_s, df_m):
    y = df_m[y_col].values
    y_s = (y - y.mean()) / y.std()
    m = sm.OLS(y_s, sm.add_constant(X_s)).fit()
    return m.params[1:]

beta1 = std_coefs('宿泊者数_per_capita', X_std, df_model)
beta2 = std_coefs('外国人比率',          X_std, df_model)

FEATURE_LABELS = ['食料費割合\n（顎）', '交通費割合\n（足）', '旅館数密度\n（枕）',
                  '客室数/施設', '気温', '高齢化率']

print("\n■ 標準化回帰係数（一人当たり宿泊者数モデル）:", beta1.round(4))
print("■ 標準化回帰係数（外国人比率モデル）:",        beta2.round(4))

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

# ────────────────────────────────────────────────────────────────
# 図1: 相関行列ヒートマップ（6変数）
# ────────────────────────────────────────────────────────────────
ALL_VARS = FEATURES + OUTCOMES
ALL_LABELS = FEATURE_LABELS + ['宿泊者数\nper capita', '外国人比率']
corr_mat = df_model[ALL_VARS].corr().values
n_vars = len(ALL_VARS)

fig1, ax1 = plt.subplots(figsize=(10, 8))
im = ax1.imshow(corr_mat, cmap='RdBu_r', vmin=-1, vmax=1)
ax1.set_xticks(range(n_vars))
ax1.set_yticks(range(n_vars))
ax1.set_xticklabels(ALL_LABELS, rotation=30, ha='right', fontsize=9)
ax1.set_yticklabels(ALL_LABELS, fontsize=9)
ax1.set_title('顎（食料費）・足（交通費）・枕（旅館密度）と観光需要の相関行列\n（SSDSE-B 2022年度 都道府県別）',
              fontsize=11, fontweight='bold')
for i in range(n_vars):
    for j in range(n_vars):
        r = corr_mat[i, j]
        txt_col = 'white' if abs(r) > 0.55 else 'black'
        ax1.text(j, i, f'{r:.2f}', ha='center', va='center',
                 fontsize=9, color=txt_col, fontweight='bold')
plt.colorbar(im, ax=ax1, shrink=0.8, label='相関係数')
plt.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2024_H5_6_fig1_corr.png'), bbox_inches='tight', dpi=150)
plt.close(fig1)
print("\n図1保存: 2024_H5_6_fig1_corr.png")

# ────────────────────────────────────────────────────────────────
# 図2: 旅館数密度 vs 一人当たり宿泊者数（散布図 + 回帰直線）
# ────────────────────────────────────────────────────────────────
x2 = df_model['旅館数密度'].values
y2 = df_model['宿泊者数_per_capita'].values
prefs = df_model['都道府県'].tolist()

r2_val, p2_val = stats.pearsonr(x2, y2)

fig2, ax2 = plt.subplots(figsize=(10, 7))
ax2.scatter(x2, y2, color='#E65100', s=50, alpha=0.7, zorder=3)

# 都道府県ラベル（全都道府県）
HIGHLIGHT = {'山梨県': '#C62828', '長野県': '#C62828', '静岡県': '#C62828',
             '北海道': '#1565C0', '沖縄県': '#1565C0', '鳥取県': '#2E7D32'}
for i, pref in enumerate(prefs):
    color = HIGHLIGHT.get(pref, '#555555')
    fontsize = 8.5 if pref in HIGHLIGHT else 7
    fw = 'bold' if pref in HIGHLIGHT else 'normal'
    ax2.annotate(pref, (x2[i], y2[i]),
                 textcoords='offset points', xytext=(4, 3),
                 fontsize=fontsize, color=color, fontweight=fw)

# ハイライト点
for pref, col in HIGHLIGHT.items():
    idx = prefs.index(pref) if pref in prefs else None
    if idx is not None:
        ax2.scatter(x2[idx], y2[idx], color=col, s=100, zorder=5)

# 回帰直線
z2 = np.polyfit(x2, y2, 1)
xs2 = np.linspace(x2.min() - 0.1, x2.max() + 0.1, 200)
ax2.plot(xs2, np.poly1d(z2)(xs2), 'r-', linewidth=2, alpha=0.8,
         label=f'回帰直線  r={r2_val:.2f}, p={p2_val:.4f}')

ax2.set_xlabel('旅館数密度（10,000人当たり旅館数）【枕の代理指標】', fontsize=11)
ax2.set_ylabel('一人当たり延べ宿泊者数（泊/人）', fontsize=11)
ax2.set_title('旅館インフラ密度（枕）と観光宿泊需要の関係\n（SSDSE-B 2022年度 都道府県別）',
              fontsize=12, fontweight='bold')
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)
plt.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2024_H5_6_fig2_scatter.png'), bbox_inches='tight', dpi=150)
plt.close(fig2)
print("図2保存: 2024_H5_6_fig2_scatter.png")

# ────────────────────────────────────────────────────────────────
# 図3: 標準化回帰係数の棒グラフ（2モデル比較）
# ────────────────────────────────────────────────────────────────
fig3, ax3 = plt.subplots(figsize=(10, 5.5))

x_pos3 = np.arange(len(FEATURES))
width = 0.38
bars1 = ax3.bar(x_pos3 - width/2, beta1, width,
                color='#E65100', alpha=0.78, edgecolor='white', label='一人当たり延べ宿泊者数')
bars2 = ax3.bar(x_pos3 + width/2, beta2, width,
                color='#1565C0', alpha=0.78, edgecolor='white', label='外国人延べ宿泊者数比率')

for bars, betas in [(bars1, beta1), (bars2, beta2)]:
    for bar, b in zip(bars, betas):
        va = 'bottom' if b >= 0 else 'top'
        offset = 0.01 if b >= 0 else -0.01
        ax3.text(bar.get_x() + bar.get_width() / 2,
                 bar.get_height() + offset,
                 f'β={b:.2f}', ha='center', va=va, fontsize=8.5, fontweight='bold')

ax3.axhline(0, color='gray', linestyle='--', linewidth=0.8)
ax3.set_xticks(x_pos3)
ax3.set_xticklabels(FEATURE_LABELS, fontsize=10)
ax3.set_ylabel('標準化回帰係数 β', fontsize=11)
ax3.set_title('顎・足・枕の標準化回帰係数\n一人当たり宿泊者数モデル vs 外国人比率モデル',
              fontsize=12, fontweight='bold')
ax3.legend(fontsize=10)
ax3.grid(axis='y', alpha=0.3)

ymax = max(abs(beta1).max(), abs(beta2).max())
ax3.set_ylim(-ymax * 1.35, ymax * 1.45)

# R²テキスト
ax3.text(0.02, 0.97,
         f'R²（宿泊者数モデル）= {model1.rsquared:.3f}\nR²（外国人比率モデル）= {model2.rsquared:.3f}',
         transform=ax3.transAxes, fontsize=9, va='top',
         bbox=dict(facecolor='#F5F5F5', edgecolor='#BDBDBD', boxstyle='round,pad=0.4'))

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

# ────────────────────────────────────────────────────────────────
# 図4: 一人当たり延べ宿泊者数 Top10都道府県の棒グラフ
# ────────────────────────────────────────────────────────────────
top10 = df_model.nlargest(10, '宿泊者数_per_capita').reset_index(drop=True)
top10_prefs  = top10['都道府県'].tolist()
top10_vals   = top10['宿泊者数_per_capita'].values
top10_fratio = top10['外国人比率'].values

fig4, ax4a = plt.subplots(figsize=(10, 6))
ax4b = ax4a.twinx()

colors4 = ['#C62828' if p == '鳥取県' else '#E65100' for p in top10_prefs]
bars4 = ax4a.bar(range(10), top10_vals, color=colors4, alpha=0.78, edgecolor='white',
                 label='一人当たり延べ宿泊者数（左軸）')
ax4b.plot(range(10), top10_fratio * 100, 'o-',
          color='#1565C0', linewidth=2, markersize=7, label='外国人比率 %（右軸）')

for i, (v, f) in enumerate(zip(top10_vals, top10_fratio)):
    ax4a.text(i, v + 0.02, f'{v:.2f}', ha='center', va='bottom', fontsize=8.5, fontweight='bold')

ax4a.set_xticks(range(10))
ax4a.set_xticklabels(top10_prefs, rotation=25, ha='right', fontsize=10)
ax4a.set_ylabel('一人当たり延べ宿泊者数（泊/人）', fontsize=11)
ax4b.set_ylabel('外国人延べ宿泊者数比率（%）', fontsize=11, color='#1565C0')
ax4b.tick_params(axis='y', colors='#1565C0')
ax4a.set_title('一人当たり延べ宿泊者数 上位10都道府県\n（SSDSE-B 2022年度／外国人比率も併記）',
               fontsize=12, fontweight='bold')

# 凡例を統合
lines1, labs1 = ax4a.get_legend_handles_labels()
lines2, labs2 = ax4b.get_legend_handles_labels()
ax4a.legend(lines1 + lines2, labs1 + labs2, fontsize=9, loc='upper right')
ax4a.grid(axis='y', alpha=0.3)

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

print("\n全図の生成完了（4枚）")