機械学習による
15歳未満人口の推定

# --- 原論文 表1「スコア（決定係数）」5回分をそのまま転記【報告値】 ---
#     ※ これは原論文がランダムフォレストで学習した結果であり、
#        本スクリプトが計算した値ではない（報告値の転記）。
scores = [0.9886389795692969, 0.9908011704351442, 0.9854055839021446,
          0.9814707632552758, 0.9905348525568049]
score_mean = 0.987370269943733           # 原論文 表1 の平均

# --- 原論文 図2・4.2節の報告統計量【報告値】 ---
trend_slope     = 0.9635                  # 図2 トレンドライン y = 0.9635x + 272.21
trend_intercept = 272.21
corr_report     = 0.9948                  # 5回平均と実際の15歳未満人口の相関係数
r2_report       = 0.9896                  # 同・決定係数
mean_rel_err    = 0.038                   # 図3 相対誤差の平均

print("原論文 表1 スコア（決定係数, 5回）:")
for i, s in enumerate(scores, 1):
    print(f"  {i}回目: {s}")
print(f"  平均  : {score_mean}")
print(f"報告トレンドライン: y = {trend_slope}x + {trend_intercept}")
print(f"報告 相関係数 r = {corr_report} / 決定係数 R^2 = {r2_report}")
print(f"報告 相対誤差の平均 = {mean_rel_err}")

# --- 原論文 表1 の5回スコアを棒グラフにする（報告値の可視化・再計算ではない） ---
fig1, ax1 = plt.subplots(figsize=(9, 5.5))
labels = [f'{i}回目' for i in range(1, 6)]
bars = ax1.bar(labels, scores, color='#4e9af1', edgecolor='white',
               linewidth=0.6, alpha=0.9, width=0.62)
ax1.axhline(score_mean, color='#c0392b', linewidth=1.8, linestyle='--',
            label=f'5回平均: {score_mean:.4f}')
for b, s in zip(bars, scores):
    ax1.text(b.get_x() + b.get_width() / 2, s + 0.0004, f'{s:.4f}',
             ha='center', va='bottom', fontsize=9.5)
ax1.set_ylim(0.975, 0.995)
ax1.set_ylabel('決定係数（スコア）', fontsize=12)
ax1.set_title('図1：ランダムフォレストの決定係数（原論文 表1・5回）\n'
              '【報告値の可視化（再計算ではない）】',
              fontsize=13, fontweight='bold', pad=12)
ax1.grid(axis='y', alpha=0.3, linewidth=0.8)
ax1.legend(fontsize=10, loc='lower right')
ax1.spines['top'].set_visible(False)
ax1.spines['right'].set_visible(False)
plt.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2018_H4_fig1.png'), bbox_inches='tight')
plt.close(fig1)

# --- 原論文 図2 の報告回帰直線 y=0.9635x+272.21 を描画（報告値の可視化） ---
#     ※ 個々のテスト推定値（散布点）は原論文に数値が無いため描かない。
#        報告された直線と統計量だけを転記して可視化する。
xr = np.linspace(0, 300000, 200)
y_trend = trend_slope * xr + trend_intercept

fig2, ax2 = plt.subplots(figsize=(8.5, 7))
ax2.plot(xr, y_trend, color='#c0392b', linewidth=2.2,
         label=f'報告トレンドライン: y = {trend_slope}x + {trend_intercept}')
ax2.plot(xr, xr, color='#888', linewidth=1.2, linestyle=':',
         label='参照線 y = x（完全一致）')
ax2.set_xlim(0, 300000)
ax2.set_ylim(0, 300000)
ax2.set_xlabel('実際の15歳未満人口（人）', fontsize=12)
ax2.set_ylabel('5回の推定値の平均（人）', fontsize=12)
ax2.set_title('図2：実際の15歳未満人口 vs 推定値の報告関係\n'
              f'【報告値の可視化】 r = {corr_report}, R² = {r2_report}',
              fontsize=13, fontweight='bold', pad=12)
ax2.text(0.04, 0.90,
         f'相関係数 r = {corr_report}\n決定係数 R² = {r2_report}\n'
         f'相対誤差の平均 = {mean_rel_err}',
         transform=ax2.transAxes, fontsize=10.5,
         bbox=dict(boxstyle='round', facecolor='#fff6e5', edgecolor='#f0a500'))
ax2.legend(fontsize=9.5, loc='lower right')
ax2.grid(alpha=0.25, linewidth=0.8)
ax2.spines['top'].set_visible(False)
ax2.spines['right'].set_visible(False)
plt.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2018_H4_fig2.png'), bbox_inches='tight')
plt.close(fig2)

# --- SSDSE-B-2026 を cp932 で読み込み、都道府県行のみ・2023年に絞る ---
df_b = pd.read_csv(DATA_B, encoding='cp932', header=1)
df_b = df_b[df_b['地域コード'].astype(str).str.match(r'^R\d{5}', na=False)].copy()
df_b['年度'] = df_b['年度'].astype(int)

df23 = df_b[df_b['年度'] == 2023].dropna(subset=['総人口', '15歳未満人口']).copy()
xp = df23['総人口'].astype(float).values / 10000        # 万人
yp = df23['15歳未満人口'].astype(float).values / 10000  # 万人

# 実データの相関と回帰直線（実python実行＝実再現値）
r_real, p_real = stats.pearsonr(xp, yp)
slope_real, intercept_real = np.polyfit(xp, yp, 1)

fig3, ax3 = plt.subplots(figsize=(9, 7))
ax3.scatter(xp, yp, s=60, color='#5cb85c', alpha=0.8,
            edgecolors='white', linewidths=0.7, zorder=3)
xr3 = np.linspace(xp.min(), xp.max(), 100)
ax3.plot(xr3, slope_real * xr3 + intercept_real, 'k-', linewidth=1.8, alpha=0.75,
         label=f'回帰直線 y = {slope_real:.3f}x + {intercept_real:.2f}')
# 突出する県にラベル
for name, xv, yv in zip(df23['都道府県'], xp, yp):
    if xv > 500 or name in ('沖縄県', '鳥取県'):
        ax3.annotate(name.replace('県', '').replace('都', '').replace('府', ''),
                     (xv, yv), xytext=(4, 3), textcoords='offset points',
                     fontsize=8, color='#333')
ax3.set_xlabel('総人口（万人・SSDSE-B 2023）', fontsize=12)
ax3.set_ylabel('15歳未満人口（万人・SSDSE-B 2023）', fontsize=12)
ax3.set_title('図3：総人口 vs 15歳未満人口（47都道府県・2023年）\n'
              f'【実SSDSE-B・実再現】 実測 r = {r_real:.4f}, R² = {r_real**2:.4f}',
              fontsize=13, fontweight='bold', pad=12)
ax3.legend(fontsize=10, loc='upper left')
ax3.grid(alpha=0.25, linewidth=0.8)
ax3.spines['top'].set_visible(False)
ax3.spines['right'].set_visible(False)
plt.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2018_H4_fig3.png'), bbox_inches='tight')
plt.close(fig3)

# --- 原論文と同じ RandomForestRegressor を、SSDSEの総人口→15歳未満人口で実学習 ---
#     全都道府県×全年（2012-2023）を使い、原論文と同じ test_size=0.2, random_state=70。
#     説明変数は SSDSE に実在する「総人口」のみ（総面積・可住地面積は現行SSDSEに無い）。
df_ml = df_b.dropna(subset=['総人口', '15歳未満人口']).copy()
X = df_ml[['総人口']].astype(float).values
y = df_ml['15歳未満人口'].astype(float).values

X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.2, random_state=70)
scaler = StandardScaler().fit(X_tr)                 # 原論文どおり標準化
model = RandomForestRegressor(random_state=70).fit(scaler.transform(X_tr), y_tr)

pred = model.predict(scaler.transform(X_te))
r2_rf = model.score(scaler.transform(X_te), y_te)   # 実測の決定係数
r_rf, _ = stats.pearsonr(y_te, pred)                # 実測の相関係数

fig4, ax4 = plt.subplots(figsize=(8.5, 7))
ax4.scatter(y_te / 10000, pred / 10000, s=45, color='#9b59b6', alpha=0.7,
            edgecolors='white', linewidths=0.6, zorder=3)
lim = max(y_te.max(), pred.max()) / 10000 * 1.05
ax4.plot([0, lim], [0, lim], color='#888', linestyle=':', linewidth=1.3,
         label='y = x（完全一致）')
ax4.set_xlim(0, lim); ax4.set_ylim(0, lim)
ax4.set_xlabel('実際の15歳未満人口（万人）', fontsize=12)
ax4.set_ylabel('ランダムフォレストの推定値（万人）', fontsize=12)
ax4.set_title('図4：ランダムフォレストによる推定（総人口のみ・都道府県×年）\n'
              f'【実SSDSE-B・実学習・実再現】 実測 R² = {r2_rf:.4f}, r = {r_rf:.4f}',
              fontsize=12.5, fontweight='bold', pad=12)
ax4.legend(fontsize=10, loc='upper left')
ax4.grid(alpha=0.25, linewidth=0.8)
ax4.spines['top'].set_visible(False)
ax4.spines['right'].set_visible(False)
plt.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, '2018_H4_fig4.png'), bbox_inches='tight')
plt.close(fig4)