自治体ごとの
ふるさと納税の必要性を定義する

# ============================================================
# 図4: 再現可能性トリアージ（この教材で扱えるデータの整理）
#   どの変数が SSDSE で再現でき、どれが原論文の報告値かを一覧表として図示。
# ============================================================
print("\n=== 図4: 再現可能性トリアージ ===")
rows = [
    ('ふるさと納税 受け入れ金額', 'ふるさと納税ガイド', '×', '報告値'),
    ('市区町村・都道府県の面積', '国土地理院 面積調2020', '×', '報告値'),
    ('市区町村別 総人口', '住民基本台帳 2017', '×', '報告値'),
    ('都道府県別 歳入額', '地方財政統計年報 H29', '×', '報告値'),
    ('転入者数 − 転出者数', 'SSDSE-2020A/B（現B）', '△', '都道府県は収録'),
    ('人口密度（＝人口／面積）', '面積が SSDSE 未収録', '×', '報告値'),
    ('歳入割合（納税／歳入）', '市区町村レベル未収録', '×', '報告値'),
]
print(f"  {'変数':<22}{'SSDSE収録':<10}{'扱い'}")
for name, src, ok, how in rows:
    print(f"  {name:<22}{ok:<10}{how}（{src}）")

fig, ax = plt.subplots(figsize=(9.5, 5))
ax.axis('off')
col_x = [0.02, 0.46, 0.70, 0.86]
headers = ['分析に用いた変数', '出典', 'SSDSE収録', 'この教材での扱い']
ax.text(0.5, 0.97, '再現可能性トリアージ：この論文の変数と SSDSE',
        ha='center', fontsize=13, fontweight='bold', transform=ax.transAxes)
for cx, h in zip(col_x, headers):
    ax.text(cx, 0.88, h, fontsize=10.5, fontweight='bold', transform=ax.transAxes, color='#1565C0')
ax.axhline(0.855, color='#1565C0', linewidth=1.2)
for i, (name, src, ok, how) in enumerate(rows):
    yy = 0.80 - i * 0.105
    face = '#E8F5E9' if ok != '×' else '#FFF5F5'
    ax.add_patch(plt.Rectangle((0.0, yy - 0.045), 1.0, 0.092, transform=ax.transAxes,
                               facecolor=face, edgecolor='none', zorder=0))
    ax.text(col_x[0], yy, name, fontsize=9.8, transform=ax.transAxes, va='center')
    ax.text(col_x[1], yy, src, fontsize=8.8, transform=ax.transAxes, va='center', color='#555')
    ok_color = '#2E7D32' if ok != '×' else '#C62828'
    ax.text(col_x[2] + 0.05, yy, ok, fontsize=12, transform=ax.transAxes, va='center',
            ha='center', color=ok_color, fontweight='bold')
    ax.text(col_x[3], yy, how, fontsize=9.2, transform=ax.transAxes, va='center')
ax.text(0.99, -0.02, '×＝SSDSE未収録のため原論文の報告値を可視化／△＝一部（都道府県）のみ収録',
        transform=ax.transAxes, ha='right', fontsize=8, color='#B71C1C', style='italic')
out4 = os.path.join(FIG_DIR, '2020_H4_fig4.png')
plt.savefig(out4, dpi=150, bbox_inches='tight')
plt.close()
print(f"  保存: {out4}")

import os
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from matplotlib.patches import Patch

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

FIG_DIR = 'html/figures'
os.makedirs(FIG_DIR, exist_ok=True)

# ============================================================
# 図1: 人口密度との相関係数（原論文の報告値）
#   原論文 5.1「仮説の検証」より。すべて原論文に印字された報告値。
# ============================================================
print("=== 図1: 相関係数（原論文の報告値） ===")
corr_labels = [
    '都道府県別\n歳入 × 人口密度',
    '都道府県別\n(転入−転出) × 人口密度',
    '市区町村別*\n歳入 × 人口密度',
    '市区町村別*\n(転入−転出) × 人口密度',
]
corr_values = [0.8468049, 0.8249227, 0.6286418, 0.5282505]  # 原論文 図3/図4/図7/図8
corr_colors = ['#1565C0', '#1565C0', '#E65100', '#E65100']
for lab, v in zip([l.replace('\n', ' ') for l in corr_labels], corr_values):
    print(f"  r = {v:.7f}  |  {lab}")

fig, ax = plt.subplots(figsize=(9, 5))
ypos = np.arange(len(corr_values))[::-1]
bars = ax.barh(ypos, corr_values, color=corr_colors, edgecolor='white', height=0.62)
for y, v in zip(ypos, corr_values):
    ax.text(v + 0.01, y, f'{v:.4f}', va='center', ha='left', fontsize=12, fontweight='bold')
ax.axvline(0.7, color='#888', linestyle='--', linewidth=1)
ax.text(0.705, len(corr_values) - 0.5, '強い相関の目安 (0.7)', fontsize=9, color='#666')
ax.set_yticks(ypos)
ax.set_yticklabels(corr_labels, fontsize=10.5)
ax.set_xlim(0, 1.0)
ax.set_xlabel('Pearson 相関係数 r（原論文の報告値）', fontsize=11)
ax.set_title('人口密度と各指標の相関（2017年・原論文の報告値）\n'
             '*市区町村別は政令指定都市を外れ値として除外', fontsize=12.5, fontweight='bold', pad=10)
legend = [Patch(facecolor='#1565C0', label='都道府県別 (n=47)'),
          Patch(facecolor='#E65100', label='市区町村別（政令市除く）')]
ax.legend(handles=legend, loc='lower right', fontsize=9.5)
ax.grid(axis='x', linestyle='--', alpha=0.4)
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.text(0.01, 0.02, '原論文の報告値の可視化（再計算ではない）',
        transform=ax.transAxes, ha='left', fontsize=8.5, color='#B71C1C', style='italic')
plt.tight_layout()
out1 = os.path.join(FIG_DIR, '2020_H4_fig1.png')
plt.savefig(out1, dpi=150, bbox_inches='tight')
plt.close()
print(f"  保存: {out1}")

# ============================================================
# 図2: 人口密度で説明のつかない4自治体（原論文 図9の報告値）
#   原論文 図9 の緑字（各自治体の人口密度[人/km^2]）より。
# ============================================================
print("\n=== 図2: 過剰にふるさと納税を集めている4自治体（原論文の報告値） ===")
outlier_names = ['上峰町', 'みやき町', '湯浅町', '泉佐野市']
outlier_density = [303.0, 490.3, 750.4, 1802.5]  # 原論文 図9 緑字
# 参考：人口密度から必要性が高い（適正）とされた低密度自治体の例（原論文 図9 緑字）
ref_names = ['北山村', '小谷村', '奈半利町', '都農町', '根室市']
ref_density = [9.4, 11.0, 16.1, 39.1, 51.0]
for n, d in zip(outlier_names, outlier_density):
    print(f"  外れ値: {n:<6} 人口密度 = {d:>7.1f} 人/km^2")
for n, d in zip(ref_names, ref_density):
    print(f"  適正例: {n:<6} 人口密度 = {d:>7.1f} 人/km^2")

fig, ax = plt.subplots(figsize=(9, 5))
ypos = np.arange(len(outlier_names))
bars = ax.barh(ypos, outlier_density, color='#C62828', edgecolor='white', height=0.6)
for y, (n, d) in enumerate(zip(outlier_names, outlier_density)):
    ax.text(d + 20, y, f'{d:,.1f}', va='center', ha='left', fontsize=12, fontweight='bold')
ax.set_yticks(ypos)
ax.set_yticklabels(outlier_names, fontsize=13)
ax.set_xlim(0, 2100)
ax.set_xlabel('人口密度（人/km²・原論文の報告値）', fontsize=11)
ax.set_title('人口密度では説明できなかった4自治体\n（歳入割合が高いのに人口密度も高い＝過剰収集の疑い）',
             fontsize=12.5, fontweight='bold', pad=10)
ax.grid(axis='x', linestyle='--', alpha=0.4)
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ref_txt = '参考（人口密度が低く必要性が高いとされた自治体）：\n' + \
          '　'.join(f'{n} {d}' for n, d in zip(ref_names, ref_density)) + ' 人/km²'
ax.text(0.98, 0.06, ref_txt, transform=ax.transAxes, ha='right', va='bottom',
        fontsize=8.5, color='#333',
        bbox=dict(boxstyle='round,pad=0.4', facecolor='#FFF9C4', edgecolor='#F9A825'))
ax.text(0.99, 0.02, '原論文 図9 の報告値（緑字）の可視化（再計算ではない）',
        transform=ax.transAxes, ha='right', fontsize=8, color='#B71C1C', style='italic')
plt.tight_layout()
out2 = os.path.join(FIG_DIR, '2020_H4_fig2.png')
plt.savefig(out2, dpi=150, bbox_inches='tight')
plt.close()
print(f"  保存: {out2}")

# ============================================================
# 図3: 必要性モデルの概念図（原論文 図10・図11 の考え方の図解）
#   ※実データではない。原論文の主張（人口密度と歳入割合は反比例が望ましい）を図解。
# ============================================================
print("\n=== 図3: 必要性モデルの概念図（実データではない・考え方の図解） ===")
fig, ax = plt.subplots(figsize=(9, 5.5))
x = np.linspace(20, 2000, 300)
y_curve = 40.0 / x  # 反比例の「成功曲線」を模した概念曲線（任意スケール）
ax.plot(x, y_curve, color='#C62828', linewidth=2.4, label='必要性の目安（反比例＝成功例の曲線）')
# 曲線上（うまく活用＝成功例）の概念点
ax.scatter([60, 120, 300], 40.0 / np.array([60, 120, 300]), s=70, color='#2E7D32',
           zorder=5, edgecolors='white', label='曲線上：制度をうまく活用（成功例）')
# 曲線から外れた4自治体（概念的に上方＝密度が高いのに歳入割合も高い）
off_x = [303, 490, 750, 1802]
off_y = [0.26, 0.30, 0.34, 0.40]
off_lab = ['上峰町', 'みやき町', '湯浅町', '泉佐野市']
ax.scatter(off_x, off_y, s=90, color='#E65100', zorder=6, marker='D',
           edgecolors='white', label='曲線から外れる：過剰収集の疑い')
for xx, yy, ll in zip(off_x, off_y, off_lab):
    ax.annotate(ll, (xx, yy), fontsize=10, xytext=(6, 4), textcoords='offset points')
ax.set_xlabel('人口密度（人/km²）→ 大きいほど「必要性が低い」', fontsize=11)
ax.set_ylabel('歳入割合（ふるさと納税 ÷ 歳入）', fontsize=11)
ax.set_title('人口密度と歳入割合：ふるさと納税の「必要性」モデル（概念図）',
             fontsize=12.5, fontweight='bold', pad=10)
ax.set_xlim(0, 2000)
ax.set_ylim(0, 0.45)
ax.legend(fontsize=9, loc='upper center')
ax.grid(linestyle='--', alpha=0.35)
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.text(0.99, 0.02, '概念図（実データではない）／グラフの実データは原論文 図9〜11参照',
        transform=ax.transAxes, ha='right', fontsize=8, color='#B71C1C', style='italic')
plt.tight_layout()
out3 = os.path.join(FIG_DIR, '2020_H4_fig3.png')
plt.savefig(out3, dpi=150, bbox_inches='tight')
plt.close()
print(f"  保存: {out3}")

# ── 要約 ─────────────────────────────────────────────────
print("\n=== 要約（すべて原論文の報告値） ===")
print("都道府県別  歳入×人口密度       r = 0.8468049（強い相関）")
print("都道府県別  (転入−転出)×人口密度 r = 0.8249227（強い相関）")
print("市区町村別  歳入×人口密度       r = 0.6286418（政令市除く）")
print("市区町村別  (転入−転出)×人口密度 r = 0.5282505（政令市除く）")
print("人口密度で説明できない4自治体：湯浅町(750.4) 上峰町(303.0) みやき町(490.3) 泉佐野市(1802.5) 人/km²")

print("\n=== 全図生成完了 ===")
for f in [out1, out2, out3, out4]:
    print(f"  {'OK' if os.path.exists(f) else 'NG'}: {os.path.basename(f)}")