import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
df = df[df['SSDSE-B-2026'] == 2023]
# 基本散布図
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
ax = axes[0, 0]
ax.scatter(df['C5401'], df['L3221'], alpha=0.7, s=60)
ax.set_xlabel('標準地価（円/㎡）'); ax.set_ylabel('消費支出（円）')
ax.set_title('① 基本散布図')

# 第3変数で色分け
ax = axes[0, 1]
sc = ax.scatter(df['C5401'], df['L3221'],
                c=df['A1101'], cmap='viridis', s=60, alpha=0.8)
plt.colorbar(sc, ax=ax, label='総人口')
ax.set_title('② 総人口で色分け')

# 回帰線つき
ax = axes[1, 0]
sns.regplot(x='C5401', y='L3221', data=df,
            ax=ax, scatter_kws={'alpha':0.6}, line_kws={'color':'red'})
ax.set_title('③ 回帰線（95% CI 帯）')

# 六角ビン（大きいデータ向け）
ax = axes[1, 1]
hb = ax.hexbin(df['C5401'], df['L3221'], gridsize=10, cmap='Blues')
plt.colorbar(hb, ax=ax, label='密度')
ax.set_title('④ 六角ビン')

plt.tight_layout(); plt.savefig('scatter_4panel.png', dpi=110)

import numpy as np
x = np.array([2,4,7,5,3])
y = np.array([3,7,12,9,4])
print(f"r: {np.corrcoef(x, y)[0,1]:.3f}")

import matplotlib.pyplot as plt

plt.figure(figsize=(8, 5))
plt.scatter(x, y, s=60, alpha=0.6, color='steelblue', edgecolors='white')
plt.xlabel('X変数')
plt.ylabel('Y変数')
plt.title('散布図のタイトル')
plt.grid(alpha=0.3)
plt.show()

# 回帰直線を追加
import numpy as np
z = np.polyfit(x, y, 1)
plt.plot(x, z[0]*x + z[1], 'r-', label=f'y={z[0]:.2f}x+{z[1]:.2f}')
plt.legend()

# ── この抜粋で使うデータを用意します（SSDSE-B の 47 都道府県・最新年度）──
# 英字の項目コード（C5401 など）も後のブロックで使うので skiprows=[1] で読み、
# 本文で使う短い呼び名は、実在する列の別名として足す。
import pandas as pd

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
df = df[df['SSDSE-B-2026'] == df['SSDSE-B-2026'].max()].copy()   # 最新年度の 47 行
df['総人口'] = df['A1101']
df['食料費'] = df['L322101']      # 食料費（二人以上の世帯）
df['教育費'] = df['L322108']      # 教育費（二人以上の世帯）
df['住居費'] = df['L322102']      # 住居費（二人以上の世帯）
# 「地方」は SSDSE に無い列なので、都道府県名から 8 地方区分を作る
_BLOCKS = {
    '北海道': ['北海道'],
    '東北': ['青森県', '岩手県', '宮城県', '秋田県', '山形県', '福島県'],
    '関東': ['茨城県', '栃木県', '群馬県', '埼玉県', '千葉県', '東京都', '神奈川県'],
    '中部': ['新潟県', '富山県', '石川県', '福井県', '山梨県', '長野県',
             '岐阜県', '静岡県', '愛知県'],
    '近畿': ['三重県', '滋賀県', '京都府', '大阪府', '兵庫県', '奈良県', '和歌山県'],
    '中国': ['鳥取県', '島根県', '岡山県', '広島県', '山口県'],
    '四国': ['徳島県', '香川県', '愛媛県', '高知県'],
    '九州・沖縄': ['福岡県', '佐賀県', '長崎県', '熊本県', '大分県', '宮崎県',
                   '鹿児島県', '沖縄県'],
}
_R = {p: g for g, ps in _BLOCKS.items() for p in ps}
df['地域'] = df['Prefecture'].map(_R)

import seaborn as sns

# 基本
sns.scatterplot(x='食料費', y='教育費', data=df)

# 色・サイズ分け
sns.scatterplot(x='食料費', y='教育費', hue='地域', size='総人口', data=df)

# 回帰直線 + 信頼区間
sns.regplot(x='食料費', y='教育費', data=df)

# 散布図行列
sns.pairplot(df[['食料費', '教育費', '住居費']], diag_kind='kde')

# joint plot
sns.jointplot(x='食料費', y='教育費', data=df, kind='reg')

import pandas as pd
import plotly.express as px

# このページは英字の項目コード（C5401 など）で通しているので、同じ流儀で読む
df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
df = df[df['SSDSE-B-2026'] == df['SSDSE-B-2026'].max()].copy()
for _c in ('A1101', 'L3221', 'L322108', 'C5401'):
    df[_c] = pd.to_numeric(df[_c], errors='coerce')

fig = px.scatter(df, x='L3221', y='L322108',      # 食料費 と 教育費
                 size='A1101',                     # 総人口
                 hover_data=['Prefecture'], trendline='ols')
fig.show()
# マウスを載せると都道府県名が出る

from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures
from sklearn.pipeline import Pipeline
import numpy as np

# 散布図に非線形回帰曲線を重ねる（多項式）
# このブロックは英字の項目コードを使うので、ここで読み直す
import pandas as pd
import matplotlib.pyplot as plt
d = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
d = d[d['SSDSE-B-2026'] == d['SSDSE-B-2026'].max()]
X = d[['C5401']].astype(float).values   # 住宅地の標準価格
y = d['L3221'].astype(float).values     # 消費支出

fig, ax = plt.subplots(figsize=(8, 5))
ax.scatter(X, y, alpha=0.6, s=50)

x_grid = np.linspace(X.min(), X.max(), 200).reshape(-1, 1)
for degree, color in [(1, 'red'), (2, 'green'), (3, 'blue')]:
    model = Pipeline([('poly', PolynomialFeatures(degree)), ('lr', LinearRegression())])
    model.fit(X, y)
    ax.plot(x_grid, model.predict(x_grid), color=color, label=f'degree={degree}')

ax.legend(); ax.set_xlabel('標準地価'); ax.set_ylabel('消費支出')
plt.savefig('scatter_poly.png', dpi=110)

from scipy import stats
from scipy.stats import gaussian_kde
import numpy as np

# 2D カーネル密度推定で散布図を補強
x = df['C5401'].values
y = df['L3221'].values
kde = gaussian_kde(np.vstack([x, y]))

xx, yy = np.mgrid[x.min():x.max():100j, y.min():y.max():100j]
positions = np.vstack([xx.ravel(), yy.ravel()])
zz = kde(positions).reshape(xx.shape)

fig, ax = plt.subplots(figsize=(8, 6))
ax.contourf(xx, yy, zz, levels=10, cmap='Blues', alpha=0.7)
ax.scatter(x, y, c='red', s=20, alpha=0.7, edgecolor='white')

# 相関と回帰
r, p = stats.pearsonr(x, y)
slope, intercept, _, _, _ = stats.linregress(x, y)
xs = np.array([x.min(), x.max()])
ax.plot(xs, intercept + slope * xs, 'k--', label=f'r={r:.2f}, p={p:.3f}')
ax.legend(); plt.savefig('scatter_kde.png', dpi=110)

# LOWESS の bandwidth 最適化（散布図に滑らかな曲線を重ねる）
import statsmodels.api as sm
import optuna

def objective(trial):
    frac = trial.suggest_float('frac', 0.1, 0.9)
    smoothed = sm.nonparametric.lowess(y, x, frac=frac, return_sorted=False)
    return ((y - smoothed)**2).mean()  # MSE 最小化

study = optuna.create_study(direction='minimize')
study.optimize(objective, n_trials=30, show_progress_bar=False)
print('Best frac:', study.best_params)

# 最適 frac で散布図にプロット
best_frac = study.best_params['frac']
smoothed = sm.nonparametric.lowess(y, x, frac=best_frac)
fig, ax = plt.subplots(figsize=(8, 5))
ax.scatter(x, y, alpha=0.6)
ax.plot(smoothed[:, 0], smoothed[:, 1], 'r-', label=f'LOWESS frac={best_frac:.2f}')
ax.legend(); plt.savefig('scatter_lowess.png', dpi=110)

# 散布図 を SSDSE-B-2026 で実行する最小コード
import pandas as pd
df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
df = df[df['SSDSE-B-2026'] == 2023]  # 2023 年のみ抽出
print(df.shape)  # (47, 112)
print(df[['Prefecture','A1101','A1303','L3221']].head())

import matplotlib.pyplot as plt
fig, axes = plt.subplots(1,2, figsize=(12,5))
axes[0].scatter(df['A1101'], df['A1303'], alpha=0.7)
axes[0].set_xlabel('A1101 総人口'); axes[0].set_ylabel('A1303 65歳以上')
axes[0].set_title('総人口 vs 高齢人口')
axes[1].scatter(df['A1101'], df['A1303']/df['A1101'], alpha=0.7, color='red')
axes[1].set_xlabel('A1101 総人口'); axes[1].set_ylabel('高齢化率')
axes[1].set_title('総人口 vs 高齢化率')
plt.tight_layout(); plt.savefig('scatter_demo.png', dpi=100)
print('r1:', df[['A1101','A1303']].corr().iloc[0,1])

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

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
d = df[df['SSDSE-B-2026'] == 2023].copy()
d['aging'] = d['A1303'] / d['A1101'] * 100        # 高齢化率 (%)

east = ['北海道','青森県','岩手県','宮城県','秋田県','山形県','福島県','茨城県','栃木県',
        '群馬県','埼玉県','千葉県','東京都','神奈川県','新潟県','富山県','石川県','福井県',
        '山梨県','長野県','岐阜県','静岡県','愛知県']
d['area'] = d['Prefecture'].apply(lambda p: '東日本' if p in east else '西日本')

fig, ax = plt.subplots(figsize=(7, 5))
for area, mk in [('東日本', 'o'), ('西日本', '^')]:
    g = d[d['area'] == area]
    ax.scatter(g['A1101'] / 1e4, g['A4101'] / 1e3,
               marker=mk, s=(g['aging'] - 20) * 12,   # 面積で高齢化率を表す
               alpha=0.7, edgecolor='white', label=f'{area} (n={len(g)})')
ax.set_xlabel('総人口（万人）'); ax.set_ylabel('出生数（千人）')
ax.set_title('マーカーの形＝地域、 面積＝高齢化率')
ax.legend(); plt.tight_layout(); plt.savefig('scatter_marker.png', dpi=120)

print(f'東日本 {len(d[d["area"]=="東日本"])} 県 / 西日本 {len(d[d["area"]=="西日本"])} 県')
print(f'高齢化率の範囲: {d["aging"].min():.1f}% 〜 {d["aging"].max():.1f}%')
print(f'マーカー面積の比 = {(d["aging"].max()-20) / (d["aging"].min()-20):.1f} 倍')
print(f'東日本の平均高齢化率 {d[d["area"]=="東日本"]["aging"].mean():.2f}% / '
      f'西日本 {d[d["area"]=="西日本"]["aging"].mean():.2f}%')

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

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
d = df[df['SSDSE-B-2026'] == 2023].copy()

# 合計特殊出生率を小数第 1 位に丸めると、 値が 8 通りしか無い離散変数になる
d['tfr_round'] = d['A4103'].round(1)
n_unique = d['tfr_round'].nunique()
dup = len(d) - n_unique
print(f'丸めた出生率の異なる値: {n_unique} 通り（47 県中 {dup} 県が他県と同じ値）')
print('最も混み合う値:', d['tfr_round'].value_counts().idxmax(),
      '→', int(d['tfr_round'].value_counts().max()), '県が完全に重なる')

rng = np.random.default_rng(20260614)      # 揺らし方を固定して毎回同じ図にする
jitter = rng.uniform(-0.02, 0.02, len(d))

fig, axes = plt.subplots(1, 2, figsize=(11, 4.5), sharey=True)
axes[0].scatter(d['tfr_round'], d['A1303'] / d['A1101'] * 100, alpha=0.8)
axes[0].set_title('ジッターなし（点が重なる）')
axes[1].scatter(d['tfr_round'] + jitter, d['A1303'] / d['A1101'] * 100, alpha=0.8)
axes[1].set_title('ジッターあり（±0.02）')
for ax in axes:
    ax.set_xlabel('合計特殊出生率（小数第 1 位に丸め）')
axes[0].set_ylabel('高齢化率（%）')
plt.tight_layout(); plt.savefig('scatter_jitter.png', dpi=120)

# 揺らしても値の意味が変わっていないことを確認する
print(f'ジッター幅 ±0.02 は刻み幅 0.1 の {0.02/0.1*100:.0f}% → 隣の目盛りとは混ざらない')
print(f'ジッター後の平均 {(d["tfr_round"]+jitter).mean():.4f} / 元の平均 {d["tfr_round"].mean():.4f}')

import pandas as pd
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from scipy import stats

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
d = df[df['SSDSE-B-2026'] == 2023]
x = d['A1101'].values / 1e4          # 総人口（万人）
y = d['A4101'].values / 1e3          # 出生数（千人）

res = stats.linregress(x, y)
print(f'回帰式  y = {res.slope:.4f} x + {res.intercept:.3f}')
print(f'相関係数 r = {res.rvalue:.4f}   決定係数 R² = {res.rvalue**2:.4f}')
print(f'傾きの標準誤差 = {res.stderr:.5f}   p 値 = {res.pvalue:.2e}')

# 回帰線と 95% 信頼帯
xs = np.linspace(x.min(), x.max(), 100)
ys = res.slope * xs + res.intercept
n = len(x)
se = np.sqrt(((y - (res.slope*x + res.intercept))**2).sum() / (n-2))
band = stats.t.ppf(0.975, n-2) * se * np.sqrt(1/n + (xs-x.mean())**2 / ((x-x.mean())**2).sum())

fig, ax = plt.subplots(figsize=(7, 5))
ax.scatter(x, y, alpha=0.7, edgecolor='white')
ax.plot(xs, ys, color='crimson', lw=2, label=f'回帰直線  R²={res.rvalue**2:.3f}')
ax.fill_between(xs, ys-band, ys+band, color='crimson', alpha=0.15, label='回帰直線の 95% 信頼帯')
ax.set_xlabel('総人口（万人）'); ax.set_ylabel('出生数（千人）')
ax.legend(); plt.tight_layout(); plt.savefig('scatter_regline.png', dpi=120)

print(f'人口 500 万人（x=500）の県の予測出生数 = {res.slope*500 + res.intercept:.2f} 千人')
print(f'信頼帯の幅（データ中心 x={x.mean():.1f} 付近）= ±{band[np.argmin(abs(xs-x.mean()))]:.3f} 千人')
print(f'信頼帯の幅（右端 x={x.max():.1f} 付近）    = ±{band[-1]:.3f} 千人')

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

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
# 年度で絞らず 47 県 × 12 年 = 564 点を重ねて描く（点が重なる状況を作る）
x = df['A1101'].values / 1e4
y = df['A4101'].values / 1e3
print(f'描画する点の数: {len(x)} 点（47 県 × 12 年）')

fig, axes = plt.subplots(1, 3, figsize=(14, 4.2), sharex=True, sharey=True)
for ax, a in zip(axes, [1.0, 0.3, 0.08]):
    ax.scatter(x, y, alpha=a, s=28, color='#1E88E5', edgecolor='none')
    ax.set_title(f'透明度 α = {a}')
    ax.set_xlabel('総人口（万人）')
axes[0].set_ylabel('出生数（千人）')
plt.tight_layout(); plt.savefig('scatter_alpha.png', dpi=120)

# 「画面上で何点が重なるか」を数える。 データ値の一致ではなく、
# 描画領域を 60x60 の升目に切って「同じ升に落ちる点」を重なりとみなす。
import numpy as np
gx = np.floor((x - x.min()) / (x.max() - x.min()) * 59).astype(int)
gy = np.floor((y - y.min()) / (y.max() - y.min()) * 59).astype(int)
cells = list(zip(gx, gy))
uniq = len(set(cells))
from collections import Counter
c = Counter(cells)
print(f'60x60 の升目に落とすと異なる位置は {uniq} 箇所 → {len(x) - uniq} 点が他の点に隠れる')
print(f'重なり率 = {(len(x)-uniq)/len(x)*100:.1f}%')
print(f'最も混み合う升には {max(c.values())} 点が重なっている')
# alpha を重ねたときの見かけの濃さ
for a, k in [(1.0, 3), (0.3, 3), (0.08, 3)]:
    print(f'alpha={a}: 同じ位置に {k} 点重なると不透明度 {1-(1-a)**k:.3f}')

import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from scipy import stats

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
d = df[df['SSDSE-B-2026'] == 2023]
x = d['A1101'].values / 1e4      # 総人口（万人）
y = d['A4101'].values / 1e3      # 出生数（千人）

# 散布図の上と右に、 各軸の分布（周辺分布）を並べる
fig = plt.figure(figsize=(7.5, 6))
gs = fig.add_gridspec(2, 2, width_ratios=(4, 1), height_ratios=(1, 4),
                      wspace=0.05, hspace=0.05)
ax = fig.add_subplot(gs[1, 0])
ax_top = fig.add_subplot(gs[0, 0], sharex=ax)
ax_right = fig.add_subplot(gs[1, 1], sharey=ax)
ax.scatter(x, y, alpha=0.75, edgecolor='white')
ax_top.hist(x, bins=15, color='#90CAF9')
ax_right.hist(y, bins=15, orientation='horizontal', color='#90CAF9')
ax_top.tick_params(labelbottom=False); ax_right.tick_params(labelleft=False)
ax.set_xlabel('総人口（万人）'); ax.set_ylabel('出生数（千人）')
plt.savefig('scatter_marginal.png', dpi=120, bbox_inches='tight')

print('=== 散布図だけでは見えない、 各軸それ自体の形 ===')
for name, v in [('総人口', x), ('出生数', y)]:
    print(f'{name}: 平均 {v.mean():7.2f}  中央値 {pd.Series(v).median():7.2f}  '
          f'歪度 {stats.skew(v):5.2f}')
print()
print(f'総人口: 平均 > 中央値 の差は {x.mean() - pd.Series(x).median():.2f} 万人')
print('→ どちらの周辺分布も強い右裾を持つ（歪度 2 以上）。')
print('   散布図の点は左下に密集し、 右上に少数が伸びる形になる。')
print(f'相関係数 r = {stats.pearsonr(x, y)[0]:.4f}（散布図が示す 2 変数の関係）')

import pandas as pd
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from scipy import stats

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
d = df[df['SSDSE-B-2026'] == 2023]
x = d['A1101'].values.astype(float)      # 総人口（人）
y = d['I5102'].values.astype(float)      # 一般診療所数

fig, axes = plt.subplots(1, 2, figsize=(11, 4.5))
axes[0].scatter(x / 1e4, y, alpha=0.75, edgecolor='white')
axes[0].set_title('通常の軸'); axes[0].set_xlabel('総人口（万人）')
axes[1].scatter(x / 1e4, y, alpha=0.75, edgecolor='white')
axes[1].set_xscale('log'); axes[1].set_yscale('log')
axes[1].set_title('両対数の軸'); axes[1].set_xlabel('総人口（万人・対数）')
axes[0].set_ylabel('一般診療所数'); axes[1].set_ylabel('一般診療所数（対数）')
plt.tight_layout(); plt.savefig('scatter_logscale.png', dpi=120)

r_lin = stats.pearsonr(x, y)[0]
r_log = stats.pearsonr(np.log10(x), np.log10(y))[0]
print(f'元のスケールでの相関 r        = {r_lin:.4f}')
print(f'両対数にしたあとの相関 r      = {r_log:.4f}')
print(f'人口の歪度  元 {stats.skew(x):.2f} → log10 {stats.skew(np.log10(x)):.2f}')
print(f'診療所数の歪度 元 {stats.skew(y):.2f} → log10 {stats.skew(np.log10(y)):.2f}')

# 両対数の傾き = べき乗則の指数
sl = stats.linregress(np.log10(x), np.log10(y))
print()
print(f'両対数での回帰: log10(診療所数) = {sl.slope:.3f} * log10(人口) + {sl.intercept:.3f}')
print(f'→ 診療所数 ∝ 人口^{sl.slope:.3f}   (指数 1 なら人口に完全比例)')
print(f'人口が 2 倍になると診療所数は {2**sl.slope:.3f} 倍')