人口規模によって異なる
保育所数・保育所在所児数・定員充足率の関係

# 原論文（小野島 2018）の本文・表から転記した報告値。いずれも市区町村（1668自治体）
# を単位とする 2015 年データの集計値であり、このブロックで数値を作り出してはいない。
# （A）人口規模別の「保育所等1つあたり在所児数」要約統計量（原論文 表2）
POP_STATS = pd.DataFrame([
    # 区分,          市区町村数, 最小値, 中央値, 平均値, 最大値, 標準偏差
    ['0〜5万人',    1125,  5.0,  74.67,  79.22, 309.00, 37.13],
    ['5〜20万人',    414, 46.0,  96.78,  99.83, 196.00, 23.16],
    ['20〜50万人',    94, 74.49, 100.62, 101.46, 154.09, 15.10],
    ['50万人以上',    35, 82.72, 101.91, 104.62, 156.54, 15.75],
], columns=['人口規模', '市区町村数', '最小値', '中央値', '平均値', '最大値', '標準偏差'])

# （B）保育所等数×在所児数の定数項なし回帰の報告値（原論文 4.2・4.3、市区町村1668）
REPORT_REG = {
    '実数':    {'beta': 98.344, 'resid_sd': 577.0, 'r2': 0.9701},
    '平方根':  {'beta': 9.6536, 'resid_sd': 4.71,  'r2': 0.9839},
}

# （C）平方根変換モデルで残差が正で大きい市区町村（原論文 4.2 の表）
RESID_TOP = pd.DataFrame([
    ['福岡県', '福岡市', 40.01, '指定都市'],
    ['広島県', '広島市', 25.68, '指定都市'],
    ['愛知県', '豊橋市', 20.28, '中核市'],
    ['大阪府', '大阪市', 17.81, '指定都市'],
    ['大阪府', '堺市',   17.31, '指定都市'],
    ['岡山県', '岡山市', 16.21, '指定都市'],
    ['京都府', '宇治市', 14.35, '市'],
    ['愛知県', '岡崎市', 14.09, '中核市'],
    ['富山県', '富山市', 13.80, '中核市'],
], columns=['都道府県', '市区町村', '残差', '都市区分'])

# （D）人口規模別「定員充足率×回帰残差」の相関係数（原論文 4.3）
CORR_BY_POP = pd.DataFrame([
    ['0〜5万人',   0.45],
    ['5〜20万人',  0.23],
    ['20〜50万人', 0.00],
    ['50万人以上', 0.11],
], columns=['人口規模', '相関係数'])

# 原論文 表2 の報告値（平均値と標準偏差）を人口規模区分ごとに棒グラフ化する。
# 元データは市区町村版 SSDSE（現行 SSDSE-B には市区町村の粒度が無い）ため、これは
# 再計算ではなく報告値の可視化。エラーバー＝標準偏差で「ばらつきの大きさ」も示す。
labels = POP_STATS['人口規模'].tolist()
means = POP_STATS['平均値'].values
sds = POP_STATS['標準偏差'].values
xpos = np.arange(len(labels))

fig1, ax1 = plt.subplots(figsize=(9, 6))
bars = ax1.bar(xpos, means, yerr=sds, capsize=8, width=0.6,
               color=['#E65100', '#1565C0', '#2E7D32', '#6A1B9A'],
               alpha=0.85, edgecolor='white', error_kw={'elinewidth': 1.8, 'ecolor': '#444'})
for xi, m, s in zip(xpos, means, sds):
    ax1.text(xi, m + s + 3, f'平均 {m:.1f}\n(SD {s:.1f})', ha='center', fontsize=10)
ax1.set_xticks(xpos)
ax1.set_xticklabels(labels, fontsize=11)
ax1.set_ylabel('保育所等1つあたりの在所児数（人）', fontsize=11)
ax1.set_ylim(0, 160)
ax1.set_title('図1 人口規模別 保育所等1つあたりの在所児数（平均±標準偏差）\n'
              '【原論文 表2 の報告値の可視化（再計算ではない・市区町村1668）】',
              fontsize=12, fontweight='bold', pad=12)
ax1.grid(axis='y', alpha=0.3)
plt.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2018_U2_fig1.png'), bbox_inches='tight')
plt.close(fig1)

# 原論文 4.2 の中心的な関係「在所児数 = β × 保育所等数（定数項なし）」を、実データで
# 再現する。原論文は市区町村1668だが、現行 SSDSE-B は都道府県47。ここでは原論文と
# 同じ 2015 年の都道府県データを使い、粒度差を明示したうえで関係の再現を試みる。
# データ仕様：cp932・先頭列 'SSDSE-B-2026' が年度・2行目(単位行)は skiprows=[1]。
df = pd.read_csv(DATA_B, encoding='cp932', skiprows=[1])
d15 = df[df['SSDSE-B-2026'] == 2015].copy()          # 原論文と同じ 2015 年
x = d15['J2503'].astype(float).values                # 保育所等数
y = d15['J2506'].astype(float).values                # 保育所等在所児数

def fit_no_intercept(xv, yv):
    """定数項なし最小二乗 β=Σxy/Σx²、分散説明率 R²=1−Σ残差²/Σy²。"""
    beta = (xv * yv).sum() / (xv * xv).sum()
    resid = yv - beta * xv
    r2 = 1.0 - (resid ** 2).sum() / (yv ** 2).sum()
    return beta, resid, r2

b_raw, res_raw, r2_raw = fit_no_intercept(x, y)                       # 実数モデル
b_sqrt, res_sqrt, r2_sqrt = fit_no_intercept(np.sqrt(x), np.sqrt(y))  # 平方根変換

fig2, (axL, axR) = plt.subplots(1, 2, figsize=(13, 6))
# 左：実数
axL.scatter(x, y, s=45, color='#1565C0', alpha=0.75, edgecolor='white', zorder=3)
xl = np.linspace(0, x.max(), 100)
axL.plot(xl, b_raw * xl, color='crimson', lw=2.2, zorder=4,
         label=f'実再現 β={b_raw:.3f}（原論文報告 98.344）')
axL.set_title('左：実数モデル', fontsize=12)
axL.set_xlabel('保育所等数 J2503（所）', fontsize=11)
axL.set_ylabel('保育所等在所児数 J2506（人）', fontsize=11)
axL.legend(loc='upper left', fontsize=10)
axL.grid(alpha=0.3)
# 右：平方根変換
axR.scatter(np.sqrt(x), np.sqrt(y), s=45, color='#2E7D32', alpha=0.75, edgecolor='white', zorder=3)
xs = np.linspace(0, np.sqrt(x).max(), 100)
axR.plot(xs, b_sqrt * xs, color='crimson', lw=2.2, zorder=4,
         label=f'実再現 β={b_sqrt:.4f}（原論文報告 9.6536）')
axR.set_title('右：平方根変換モデル', fontsize=12)
axR.set_xlabel('√保育所等数', fontsize=11)
axR.set_ylabel('√保育所等在所児数', fontsize=11)
axR.legend(loc='upper left', fontsize=10)
axR.grid(alpha=0.3)
fig2.suptitle('図2 保育所等数 × 保育所等在所児数（定数項なし回帰）\n'
              '【実SSDSE-B-2026・2015年・都道府県47の実再現：市区町村1668とは粒度が異なる】',
              fontsize=12.5, fontweight='bold')
plt.tight_layout(rect=[0, 0, 1, 0.94])
fig2.savefig(os.path.join(FIG_DIR, '2018_U2_fig2.png'), bbox_inches='tight')
plt.close(fig2)

# 原論文 4.2 の「平方根変換モデルで残差が正に大きい市区町村」の報告値を横棒で並べる。
# 残差＝保育所等数では説明できなかった在所児数のずれ。市区町村版データが現行SSDSEに
# 無いため再計算はできず、原論文の報告値をそのまま可視化する。都市区分で色分け。
resid_df = RESID_TOP.sort_values('残差')
citycolor = {'指定都市': '#C62828', '中核市': '#E65100', '市': '#1565C0'}
colors = [citycolor[c] for c in resid_df['都市区分']]
ylab = [f"{r['市区町村']}（{r['都道府県']}）" for _, r in resid_df.iterrows()]

fig3, ax3 = plt.subplots(figsize=(9.5, 6.5))
ax3.barh(range(len(resid_df)), resid_df['残差'], color=colors, alpha=0.88, edgecolor='white')
ax3.set_yticks(range(len(resid_df)))
ax3.set_yticklabels(ylab, fontsize=10)
for i, v in enumerate(resid_df['残差']):
    ax3.text(v + 0.4, i, f'{v:.2f}', va='center', fontsize=9)
handles = [plt.Rectangle((0, 0), 1, 1, color=c) for c in citycolor.values()]
ax3.legend(handles, citycolor.keys(), title='都市区分', loc='lower right', fontsize=10)
ax3.set_xlabel('平方根変換モデルの正の残差（原論文 報告値）', fontsize=11)
ax3.set_xlim(0, 46)
ax3.set_title('図3 保育所等数で説明されない在所児数（正の残差）が大きい市区町村\n'
              '【原論文 4.2 の報告値の可視化（再計算ではない）】',
              fontsize=12, fontweight='bold', pad=12)
ax3.grid(axis='x', alpha=0.3)
plt.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2018_U2_fig3.png'), bbox_inches='tight')
plt.close(fig3)

# 原論文 4.3 の結論を支える数値。人口規模ごとに「定員充足率」と「回帰残差」の相関係数
# を報告値のまま棒グラフ化する。定員充足率は社会福祉施設等調査由来で SSDSE に列が
# 無く再現不能のため、報告値の可視化にとどめる。人口規模が小さいほど相関が強い。
cpos = np.arange(len(CORR_BY_POP))
cvals = CORR_BY_POP['相関係数'].values
fig4, ax4 = plt.subplots(figsize=(9, 6))
bars = ax4.bar(cpos, cvals, width=0.6,
               color=['#E65100', '#F9A825', '#90A4AE', '#78909C'],
               alpha=0.9, edgecolor='white')
for xi, v in zip(cpos, cvals):
    ax4.text(xi, v + 0.012, f'r = {v:.2f}', ha='center', fontsize=12, fontweight='bold')
ax4.set_xticks(cpos)
ax4.set_xticklabels(CORR_BY_POP['人口規模'], fontsize=11)
ax4.set_ylabel('定員充足率 × 回帰残差 の相関係数', fontsize=11)
ax4.set_ylim(0, 0.55)
ax4.axhline(0, color='#333', lw=0.8)
ax4.set_title('図4 人口規模別「定員充足率 × 回帰残差」の相関係数\n'
              '【原論文 4.3 の報告値の可視化（再計算ではない）】',
              fontsize=12, fontweight='bold', pad=12)
ax4.grid(axis='y', alpha=0.3)
plt.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, '2018_U2_fig4.png'), bbox_inches='tight')
plt.close(fig4)