import pandas as pd
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error

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

# サンプリングバイアス: 大都市 3 県だけで学習
train_pref = ['東京都','大阪府','愛知県']
test_pref  = ['秋田県','高知県','島根県']

train = df_2023[df_2023['Prefecture'].isin(train_pref)]
test_local = df_2023[df_2023['Prefecture'].isin(test_pref)]
test_all = df_2023[~df_2023['Prefecture'].isin(train_pref)]

feats = ['A1101','A1303','A9101']
model = LinearRegression().fit(train[feats], train['A4103'])

mse_metro = mean_squared_error(train['A4103'], model.predict(train[feats]))
mse_local = mean_squared_error(test_local['A4103'], model.predict(test_local[feats]))
mse_all   = mean_squared_error(test_all['A4103'], model.predict(test_all[feats]))
print(f"訓練(大都市3県) MSE = {mse_metro:.5f}")
print(f"地方3県         MSE = {mse_local:.5f}")
print(f"全44県          MSE = {mse_all:.5f}")
print(f"\n大都市->地方 性能劣化倍率: {mse_local/max(mse_metro,1e-9):.1f}倍")

# 部分集団バイアス: 「人口規模」グループごとに性能差を測る
import pandas as pd
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import KFold
import numpy as np

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
df_2023 = df[df['SSDSE-B-2026']==2023].copy()
df_2023['size_g'] = pd.qcut(df_2023['A1101'], 3, labels=['小','中','大'])

feats = ['A1303','A9101','B4101']
y = df_2023['A4103']
oof = np.zeros(len(df_2023))
kf = KFold(n_splits=5, shuffle=True, random_state=42)
for tr,te in kf.split(df_2023):
    m = LinearRegression().fit(df_2023[feats].iloc[tr], y.iloc[tr])
    oof[te] = m.predict(df_2023[feats].iloc[te])

df_2023['pred'] = oof
for g in ['小','中','大']:
    sub = df_2023[df_2023['size_g']==g]
    mse = mean_squared_error(sub['A4103'], sub['pred'])
    print(f"人口 {g} 県 ({len(sub)}件): MSE = {mse:.5f}")

import pandas as pd
import numpy as np

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

# 47 県を 1 票ずつ数える平均（県の平均）と、人口で重み付けした平均（全国の値）
print(f"高齢化率   県の単純平均 {aging.mean():.1f}%   人口で重み付け {d['A1303'].sum() / d['A1101'].sum() * 100:.1f}%")
print(f"合計特殊出生率 県の単純平均 {d['A4103'].mean():.3f}   総人口で重み付け {np.average(d['A4103'], weights=d['A1101']):.3f}")

# 人口の少ない半分の県が、県の平均と人口の中でどれだけの重みを持つか
small = d['A1101'] <= d['A1101'].median()
print(f"人口が中央値以下の {small.sum()} 県: 県の数では {small.mean() * 100:.0f}%、人口では {d.loc[small, 'A1101'].sum() / d['A1101'].sum() * 100:.0f}%")

import pandas as pd
import numpy as np

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

def r(sub, x):
    return np.corrcoef(sub[x], sub['A4103'])[0, 1]

big = d.nlargest(20, 'A1101')      # 人口の多い 20 都道府県だけ観測できたとする
small = d.nsmallest(20, 'A1101')   # 人口の少ない 20 県だけ観測できたとする
for x, name in [('高齢化率', '高齢化率'), ('B4101', '年平均気温')]:
    print(f"{name} × 合計特殊出生率:  47 県 r = {r(d, x):+.2f}   人口上位 20 r = {r(big, x):+.2f}   人口下位 20 r = {r(small, x):+.2f}")
print('人口下位 20 の高齢化率の範囲:', round(small['高齢化率'].min(), 1), '〜', round(small['高齢化率'].max(), 1), '%')
print('人口上位 20 の高齢化率の範囲:', round(big['高齢化率'].min(), 1), '〜', round(big['高齢化率'].max(), 1), '%')

import pandas as pd
import numpy as np

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
d = df[df['SSDSE-B-2026'] == 2023].copy()
d['高齢化率'] = d['A1303'] / d['A1101'] * 100
y = d['A4103']                          # 合計特殊出生率（47 県とも実測がある）
true_mean = y.mean()
print(f'47 県の真の平均: {true_mean:.3f}')

# (1) 無作為に 12 県が欠ける（MCAR）: 1000 回くり返してずれの分布を見る
rng = np.random.default_rng(0)
shift = np.array([y.drop(rng.choice(d.index, 12, replace=False)).mean() - true_mean for _ in range(1000)])
print(f'無作為に欠ける（MCAR）   ずれの平均 {shift.mean():+.3f}  95% の範囲 {np.percentile(shift, 2.5):+.3f} 〜 {np.percentile(shift, 97.5):+.3f}')

# (2) 高齢化率の高い 12 県が欠ける（MAR）、(3) 出生率の低い 12 県が欠ける（MNAR）
mar = d.nlargest(12, '高齢化率').index
mnar = d.nsmallest(12, 'A4103').index
print(f'高齢化率の高い県が欠ける（MAR） ずれ {y.drop(mar).mean() - true_mean:+.3f}')
print(f'出生率の低い県が欠ける（MNAR）  ずれ {y.drop(mnar).mean() - true_mean:+.3f}')

# MAR の欠けを、残り 35 県で当てはめた「高齢化率 → 出生率」の回帰で埋める
obs = d.drop(mar)
b1, b0 = np.polyfit(obs['高齢化率'], obs['A4103'], 1)
filled = y.copy()
filled[mar] = b0 + b1 * d.loc[mar, '高齢化率']
print(f'MAR を回帰で補完（傾き {b1:+.4f}）した平均のずれ {filled.mean() - true_mean:+.3f}')
print(f'欠けた 12 県の実測平均 {y[mar].mean():.3f}  補完値の平均 {filled[mar].mean():.3f}')

# 時系列バイアス: 古い年度で学習し新しい年度で評価
import pandas as pd
from sklearn.linear_model import LinearRegression
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import mean_squared_error

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

feats = ['A1101','A1303','A9101','B4101']
y_col = 'A4103'

# 古い年度 (2018-2020) で学習、 新しい年度 (2021-2023) で評価
train = df[df['year'].isin([2018,2019,2020])].dropna(subset=feats+[y_col])
test  = df[df['year'].isin([2021,2022,2023])].dropna(subset=feats+[y_col])

# 人口（〜10^6）と気温（〜10）は桁が違うので、標準化してから回帰する
m = make_pipeline(StandardScaler(), LinearRegression()).fit(train[feats], train[y_col])
pred_tr = m.predict(train[feats])
pred_te = m.predict(test[feats])
print(f"訓練(2018-2020) MSE = {mean_squared_error(train[y_col], pred_tr):.5f}")
print(f"テスト(2021-2023) MSE = {mean_squared_error(test[y_col], pred_te):.5f}")
print(f"\n劣化率: {mean_squared_error(test[y_col], pred_te) / mean_squared_error(train[y_col], pred_tr):.2f}倍")

# 部分集団: 各年度ごとに性能を見る
for yr in [2018, 2019, 2020, 2021, 2022, 2023]:
    sub = df[df['year']==yr].dropna(subset=feats+[y_col])
    if len(sub) > 0:
        mse = mean_squared_error(sub[y_col], m.predict(sub[feats]))
        print(f" 年度 {yr}: MSE = {mse:.5f}")

# 必要なパッケージをインストール
# pip install pandas numpy scikit-learn matplotlib seaborn

# データを読み込む
import pandas as pd
df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
print(f"形状: {df.shape}")
print(f"年度: {df['SSDSE-B-2026'].unique()}")

import pandas as pd
from sklearn.linear_model import Ridge
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.model_selection import cross_val_score, KFold

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
df_2023 = df[df['SSDSE-B-2026'] == 2023].sort_values('Code')   # 北海道→沖縄県の順
X = df_2023[['A1101', 'A1303', 'A9101', 'B4101', 'B4106']]
y = df_2023['A4103']   # 合計特殊出生率

pipe = Pipeline([('sc', StandardScaler()), ('reg', Ridge(alpha=1.0))])

# 行順のまま 5 分割 = 各 fold が「地方のまとまり」になる
s_seq = cross_val_score(pipe, X, y, cv=5, scoring='r2')
# シャッフルして 5 分割 = 各 fold に全国の県が混ざる
s_shf = cross_val_score(pipe, X, y, cv=KFold(5, shuffle=True, random_state=0), scoring='r2')

print('行順の 5 分割   R^2 =', s_seq.round(3), f' 平均 {s_seq.mean():.3f}')
print('シャッフル 5 分割 R^2 =', s_shf.round(3), f' 平均 {s_shf.mean():.3f}')
print('行順の各 fold の県:', [df_2023['Prefecture'].iloc[te].iloc[[0, -1]].tolist() for _, te in KFold(5).split(X)])

# 部分集団ごとの誤差: 都市・地方で予測の当たり方に差があるかを測る
import pandas as pd
import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import LeaveOneOut

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

# 都市 (人口 200万以上) vs 地方
df_2023['group'] = (df_2023['A1101']>=2_000_000).map({True:'都市', False:'地方'})

feats = ['A1303','A9101','B4101']
y = df_2023['A4103'].values
preds = np.zeros(len(df_2023))
for tr, te in LeaveOneOut().split(df_2023):
    # 人口（〜10^6）と気温（〜10）は桁が違うので標準化してから回帰する
    m = make_pipeline(StandardScaler(), LinearRegression()).fit(df_2023[feats].iloc[tr], y[tr])
    preds[te] = m.predict(df_2023[feats].iloc[te])

mse = {}
for g in ['都市','地方']:
    sub = df_2023[df_2023['group']==g]
    mse[g] = mean_squared_error(sub['A4103'], preds[sub.index])
    print(f"{g} ({len(sub)}件): MSE = {mse[g]:.5f}")

# 誤差の比 = 都市 MSE / 地方 MSE（1 なら両群で同じくらい当たっている）
print(f"\n誤差の比 (都市/地方): {mse['都市']/mse['地方']:.2f}")

# ラベルノイズバイアス: 出生率に意図的ノイズを入れて再学習
import pandas as pd
import numpy as np
from sklearn.linear_model import Ridge
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import LeaveOneOut
np.random.seed(0)   # 実行のたびに同じ結果が出るようにする

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

feats = ['A1101','A1303','A9101','B4101']
X = df_2023[feats].values
y = df_2023['A4103'].values

# ラベルノイズ強度を変えて性能を測る (合成データ生成ではなく、 実値 y に観測ノイズを 加える)
for noise in [0.0, 0.05, 0.10, 0.20, 0.40]:
    rs = np.random.RandomState(0)
    y_noisy = y + rs.normal(0, noise, size=y.shape)
    loo = LeaveOneOut()
    preds = np.zeros(len(y))
    for tr, te in loo.split(X):
        m = Ridge(alpha=1.0).fit(X[tr], y_noisy[tr])
        preds[te] = m.predict(X[te])
    mse = mean_squared_error(y, preds)
    print(f"ラベルノイズ σ={noise:.2f}: 元 y との MSE = {mse:.4f}")