市区町村別でみる
合計特殊出生率推移の特徴分析

# ===== ステップ2: 図1 都道府県別 TFR 推移（実データ・沖縄強調） =====
pv = df.pivot_table(index='年度', columns='都道府県', values='TFR')
fig, ax = plt.subplots(figsize=(9.5, 5.8))
for pref in pv.columns:
    if pref not in ('沖縄県', '東京都'):
        ax.plot(pv.index, pv[pref], color='#B0BEC5', lw=0.8, alpha=0.7, zorder=1)
ax.plot(pv.index, pv['沖縄県'], color='#E65100', lw=3.0, zorder=3,
        marker='o', ms=4, label='沖縄県（一貫して最高水準）')
ax.plot(pv.index, pv['東京都'], color='#1565C0', lw=2.2, zorder=2,
        marker='s', ms=3.5, label='東京都（最低水準で推移）')
ax.set_xlabel('年')
ax.set_ylabel('合計特殊出生率（TFR）')
ax.set_title('都道府県別 合計特殊出生率の推移（SSDSE-B-2026 実データ、2012〜2023年）\n'
             '原論文 図1・図2（SSDSE-2019B、2005〜2016年）の再表現：沖縄の高さが際立つ',
             fontsize=11)
ax.legend(loc='upper right', fontsize=10)
ax.grid(alpha=0.3)
fig.tight_layout()
fig.savefig(f'{FIG_DIR}/2019_U4_fig1.png', bbox_inches='tight')
plt.close(fig)

tfr15 = pv.loc[2015].sort_values(ascending=False)
tfr23 = pv.loc[2023].sort_values(ascending=False)
print('図1 保存: 2019_U4_fig1.png')
print(f"2015年 TFR 上位: {'、'.join(f'{k} {v:.2f}' for k, v in tfr15.head(3).items())}")
print(f"2015年 TFR 下位: {'、'.join(f'{k} {v:.2f}' for k, v in tfr15.tail(3).items())}")
print(f"都道府県単位の TFR 標準偏差（実データ計算）: "
      f"2015年 {tfr15.std(ddof=1):.4f} ／ 2023年 {tfr23.std(ddof=1):.4f}")
print('  ※ 原論文の報告値 0.1327242 は「2010年の都道府県単位 TFR」の標準偏差。')
print('     原論文はこれが市区町村単位の分布幅の約半分だと指摘した（図8・図9）')

import os
import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

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

FIG_DIR = 'html/figures'
DATA_B = 'data/raw/SSDSE-B-2026.csv'
os.makedirs(FIG_DIR, exist_ok=True)

# ===== ステップ1: SSDSE-B（都道府県・時系列）の読み込み =====
df = pd.read_csv(DATA_B, encoding='cp932', header=1)
df = df[df['地域コード'].astype(str).str.match(r'^R\d{5}', na=False)].copy()
df['年度'] = df['年度'].astype(int)
df['TFR'] = df['合計特殊出生率'].astype(float)
years = sorted(df['年度'].unique())
print(f"【SSDSE-B-2026】{df['都道府県'].nunique()}都道府県 × "
      f"{len(years)}年（{years[0]}〜{years[-1]}年）")
print('  ※ 原論文が使った SSDSE-2019B は 2005〜2016 年収録。現行版とは期間が異なる')
print('  ※ 市区町村別 TFR（ベイズ推定）・年齢別女性人口は SSDSE 未収録')
print('     → 原論文の子ども女性比による近似計算は再計算不可（報告値で提示）')

# ===== ステップ3: 図2 原論文 表3 の報告値（年齢範囲別 R²・R）を可視化 =====
# 出典: 原論文 表3（2010年・全1734市区町村、TFR（ベイズ推定）×子ども女性比）
t3 = pd.DataFrame({
    '年齢範囲': ['15〜44歳', '15〜49歳', '20〜39歳', '20〜44歳', '20〜49歳', '25〜39歳'],
    '決定係数R2': [0.605, 0.576, 0.645, 0.657, 0.625, 0.623],
    '相関係数R':  [0.777, 0.759, 0.803, 0.810, 0.791, 0.789],
})
best = t3.loc[t3['決定係数R2'].idxmax()]
fig, ax = plt.subplots(figsize=(8.6, 5.2))
colors = ['#90A4AE' if a != best['年齢範囲'] else '#E65100' for a in t3['年齢範囲']]
bars = ax.bar(t3['年齢範囲'], t3['決定係数R2'], color=colors, width=0.6)
for b, r2, r in zip(bars, t3['決定係数R2'], t3['相関係数R']):
    ax.text(b.get_x() + b.get_width() / 2, b.get_height() + 0.008,
            f'R²={r2:.3f}\nR={r:.3f}', ha='center', fontsize=9)
ax.set_ylim(0.5, 0.72)
ax.set_ylabel('決定係数 R²')
ax.set_xlabel('子ども女性比の分母に用いる「若年女性」の年齢範囲')
ax.set_title('子ども女性比と TFR（ベイズ推定）の当てはまり（女性年齢範囲別）\n'
             '原論文 表3・図3 の報告値の可視化（2010年・全1734市区町村。再計算ではない）',
             fontsize=11)
ax.grid(axis='y', alpha=0.3)
fig.tight_layout()
fig.savefig(f'{FIG_DIR}/2019_U4_fig2.png', bbox_inches='tight')
plt.close(fig)
print('図2 保存: 2019_U4_fig2.png')
print(t3.to_string(index=False))
print(f"→ 最適な年齢範囲は {best['年齢範囲']}（R²={best['決定係数R2']}）")
print('→ さらに人口1万人以上の1248市区町村に限定すると R²=0.8424（原論文 図4・表5）')
print('   回帰式: TFR ＝ 4.8949 × 子ども女性比 ＋ 0.1704（原論文 表5 の報告値）')

# ===== ステップ4: 図3 原論文 表7 の報告値（平均・変動係数の年次推移）を可視化 =====
# 出典: 原論文 表7・図10（TFR（子ども女性比推定）、市区町村単位）
t7_year = [1995, 2000, 2005, 2010, 2015]
t7_mean = [1.690327, 1.663421, 1.451268, 1.50186, 1.632668]
t7_cv   = [0.155119, 0.143068, 0.156253, 0.151681, 0.166862]
fig, ax1 = plt.subplots(figsize=(8.6, 5.2))
ax1.plot(t7_year, t7_mean, color='#1565C0', marker='o', lw=2.5, label='平均')
ax1.set_xlabel('年')
ax1.set_ylabel('TFR（子ども女性比推定）の平均', color='#1565C0')
ax1.tick_params(axis='y', labelcolor='#1565C0')
ax1.set_xticks(t7_year)
ax1.set_ylim(1.40, 1.75)
ax2 = ax1.twinx()
ax2.plot(t7_year, t7_cv, color='#2E7D32', marker='s', lw=2.5, ls='--', label='変動係数')
ax2.set_ylabel('変動係数（標準偏差÷平均）', color='#2E7D32')
ax2.tick_params(axis='y', labelcolor='#2E7D32')
ax2.set_ylim(0.10, 0.20)
for y, m in zip(t7_year, t7_mean):
    ax1.annotate(f'{m:.3f}', (y, m), textcoords='offset points',
                 xytext=(0, 9), ha='center', fontsize=9, color='#1565C0')
ax1.set_title('TFR（子ども女性比推定）の平均と変動係数の年次推移\n'
              '原論文 表7・図10 の報告値の可視化（再計算ではない）', fontsize=11)
ax1.grid(alpha=0.3)
h1, l1 = ax1.get_legend_handles_labels()
h2, l2 = ax2.get_legend_handles_labels()
ax1.legend(h1 + h2, l1 + l2, loc='lower right', fontsize=10)
fig.tight_layout()
fig.savefig(f'{FIG_DIR}/2019_U4_fig3.png', bbox_inches='tight')
plt.close(fig)
print('図3 保存: 2019_U4_fig3.png')
print(pd.DataFrame({'年': t7_year, '平均': t7_mean, '変動係数': t7_cv}).to_string(index=False))
print('→ 平均は 2005 年に最小（1.451）となり 2010 年には約 0.2 上昇。')
print('   変動係数は平均が増減してもあまり変化しない（原論文の指摘）')

# ===== ステップ5: 図4 原論文 表8 の報告値（9領域の分割と市区町村数）を再表現 =====
# 出典: 原論文 表8（人口2000人以上の1424市区町村）
#   A = 2015年 TFR（子ども女性比推定）、B = その差（2015年−2005年）
#   列: A<1.3（領域3-x）／1.3≦A<1.8（領域2-x）／A≧1.8（領域1-x）
#   行: B>0.5（x-1）／0.1<B≦0.5（x-2）／B<0.1（x-3）
counts = np.array([[0, 9, 35],      # B > 0.5
                   [35, 725, 217],  # 0.1 < B ≦ 0.5
                   [56, 329, 9]])   # B < 0.1
region = [['領域3-1', '領域2-1', '領域1-1'],
          ['領域3-2', '領域2-2', '領域1-2'],
          ['領域3-3', '領域2-3', '領域1-3']]
note = [['', '東京3区（千代田・\n中央・港）を含む\nTFRが2005年から倍増', 'V字回復・高水準\n九州・沖縄に多い'],
        ['', '', ''],
        ['低水準・回復なし\n関東以北に多い', '', '']]
fig, ax = plt.subplots(figsize=(9.2, 6.4))
im = ax.imshow(np.log1p(counts), cmap='Blues', vmin=0, vmax=np.log1p(counts).max() * 1.15)
for i in range(3):
    for j in range(3):
        c = 'white' if counts[i, j] > 300 else '#0D2B45'
        ax.text(j, i - 0.24, region[i][j], ha='center', fontsize=12,
                fontweight='bold', color=c)
        ax.text(j, i, f'{counts[i, j]} 自治体', ha='center', fontsize=12, color=c)
        if note[i][j]:
            ax.text(j, i + 0.28, note[i][j], ha='center', fontsize=8.2, color=c)
ax.set_xticks([0, 1, 2])
ax.set_xticklabels(['A < 1.3', '1.3 ≦ A < 1.8', 'A ≧ 1.8'])
ax.set_yticks([0, 1, 2])
ax.set_yticklabels(['B > 0.5', '0.1 < B ≦ 0.5', 'B < 0.1'])
ax.set_xlabel('A ＝ 2015年 TFR（子ども女性比推定）')
ax.set_ylabel('B ＝ TFR（子ども女性比推定）の差\n（2015年 − 2005年）')
ax.set_title('散布図平面の9領域分割と各領域の市区町村数（人口2000人以上・計1424自治体）\n'
             '原論文 表8 の報告値の再表現（散布図そのもの＝原論文 図11 は原論文参照）',
             fontsize=11)
fig.tight_layout()
fig.savefig(f'{FIG_DIR}/2019_U4_fig4.png', bbox_inches='tight')
plt.close(fig)
print('図4 保存: 2019_U4_fig4.png')
tot = counts.sum()
print(f'9領域の市区町村数の単純和: {tot}')
print('  ※ 原論文 表8 の「合計」欄は 1424。各領域の値の単純和 1415 とは一致しない')
print('    （原論文内の不整合。本ページは各セルの値をそのまま転記している）')
print('着目された3領域（原論文 第4章）:')
print('  領域1-1（A≧1.8 かつ B>0.5, 35自治体）: 高水準かつ伸びも大。九州・沖縄に多い')
print('  領域2-1（1.3≦A<1.8 かつ B>0.5, 9自治体）: 伸びが特に大。東京3区は2005年から倍増')
print('  領域3-3（A<1.3 かつ B<0.1, 56自治体）: 低水準で回復なし。関東以北に多い')
print('done')