import numpy as np
import pandas as pd

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
df['age'] = (df['A1303'] / df['A1101'] * 100).round(1)   # 高齢化率（%、小数 1 桁）
df['dr'] = (df['A4200'] / df['A1101'] * 1000).round(2)   # 死亡率（‰、小数 2 桁）
pref = ['東京都', '沖縄県', '大阪府', '島根県', '秋田県']
s = df[df['SSDSE-B-2026'] == 2019].set_index('Prefecture').loc[pref]
t = df[df['SSDSE-B-2026'] == 2023].set_index('Prefecture').loc[pref]

# Step 1: source（2019 年度の 5 県）で直線 y = a + b x
x, y = s['age'].to_numpy(), s['dr'].to_numpy()
b = ((x - x.mean()) * (y - y.mean())).sum() / ((x - x.mean()) ** 2).sum()
a = y.mean() - b * x.mean()
print(f'x̄ = {x.mean():.2f}, ȳ = {y.mean():.3f}, b = {b:.4f}, a = {a:.3f}')

# Step 2: そのまま 2023 年度に当てる
pred = a + b * t['age'].to_numpy()
err = pred - t['dr'].to_numpy()
print(pd.DataFrame({'高齢化率': t['age'], '死亡率': t['dr'], '予測': pred.round(2), '誤差': err.round(2)}))
print(f'平均誤差 {err.mean():+.3f}  RMSE {np.sqrt((err**2).mean()):.3f}')

# Step 3: 2023 年度のラベルを 2 県（東京都・秋田県）だけ使って切片を補正
shift = -err[[0, 4]].mean()
err2 = pred + shift - t['dr'].to_numpy()
print(f'切片の補正量 {shift:+.3f}  残り 3 県の RMSE: 補正前 {np.sqrt((err[1:4]**2).mean()):.3f} → 補正後 {np.sqrt((err2[1:4]**2).mean()):.3f}')

from transformers import AutoModelForSequenceClassification, AutoTokenizer
import torch
m = AutoModelForSequenceClassification.from_pretrained('bert-base-multilingual-cased', num_labels=3)
tok = AutoTokenizer.from_pretrained('bert-base-multilingual-cased')
# ヘッド以外を凍結する例
for p in m.bert.parameters():
    p.requires_grad = False
print('trainable:', sum(p.numel() for p in m.parameters() if p.requires_grad))

import pandas as pd
from sklearn.linear_model import Ridge
from sklearn.metrics import r2_score

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='shift_jis', skiprows=[1])
d22 = df[df['SSDSE-B-2026'] == 2022]
d23 = df[df['SSDSE-B-2026'] == 2023]

m_source = Ridge(alpha=1.0).fit(d22[['A1101']], d22['H1800'])
print('source 2022 R²:', round(m_source.score(d22[['A1101']], d22['H1800']), 3))
print('そのまま 2023 に適用 R²:', round(r2_score(d23['H1800'], m_source.predict(d23[['A1101']])), 3))

import pandas as pd
from sklearn.linear_model import SGDRegressor
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import make_pipeline
from sklearn.metrics import r2_score

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='shift_jis', skiprows=[1])
d22 = df[df['SSDSE-B-2026'] == 2022]
d23 = df[df['SSDSE-B-2026'] == 2023]
scaler = StandardScaler().fit(d22[['A1101']])
m = SGDRegressor(warm_start=True, max_iter=200, random_state=0)
m.fit(scaler.transform(d22[['A1101']]), d22['H1800'])
base = r2_score(d23['H1800'], m.predict(scaler.transform(d23[['A1101']])))
m.partial_fit(scaler.transform(d23[['A1101']].iloc[:10]), d23['H1800'].iloc[:10])
after = r2_score(d23['H1800'], m.predict(scaler.transform(d23[['A1101']])))
print(f'fine-tune 前 R²: {base:.3f}')
print(f'fine-tune 後 R²: {after:.3f}')

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

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='shift_jis', skiprows=[1])
d = df[df['SSDSE-B-2026'] == 2023]
big = d[d['Prefecture'].isin(['東京都','大阪府','愛知県','神奈川県'])]
small = d[~d['Prefecture'].isin(['東京都','大阪府','愛知県','神奈川県'])]

m = LinearRegression().fit(big[['A1101']], big['H1800'])
print('source 大都市圏 R²:', round(m.score(big[['A1101']], big['H1800']), 3))
print('target 地方圏に転移 R²:', round(r2_score(small['H1800'], m.predict(small[['A1101']])), 3))

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

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='shift_jis', skiprows=[1])
d = df[df['SSDSE-B-2026'] == 2023]
src = LinearRegression().fit(d[['A1101']], d['H1800'])  # 人口 → 着工新設住宅戸数
print('Source R²:', round(src.score(d[['A1101']], d['H1800']), 3))
# 同モデル（係数・切片）の predict を気温予測に「無理やり流用」
pred_temp = src.predict(d[['A1101']])
print('Target R² (人口モデルで気温予測): ', round(r2_score(d['B4102'], pred_temp), 3))

import numpy as np
import pandas as pd

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
df['age'] = df['A1303'] / df['A1101'] * 100        # 高齢化率（%）
df['dr'] = df['A4200'] / df['A1101'] * 1000        # 死亡率（人口千人当たり）
df['lp'] = np.log10(df['A1101'])                   # log10 総人口
F = ['age', 'lp', 'A4103', 'B4101']                # 高齢化率・log 人口・合計特殊出生率・年平均気温

src = df[df['SSDSE-B-2026'] <= 2019]               # source: 2012〜2019 年度（47 県 × 8 年度）
tgt = df[df['SSDSE-B-2026'] == 2023]               # target: 2023 年度の 47 県
mu, sd = src[F].mean(), src[F].std()               # 標準化は source の平均・標準偏差で固定
Xs = np.c_[((src[F] - mu) / sd).to_numpy(), np.ones(len(src))]
Xt = np.c_[((tgt[F] - mu) / sd).to_numpy(), np.ones(len(tgt))]
ys, yt = src['dr'].to_numpy(), tgt['dr'].to_numpy()

w_src = np.linalg.lstsq(Xs, ys, rcond=None)[0]      # source で最小二乗
print('source の行数:', len(src), ' target の県数:', len(tgt))
print('重み（高齢化率, log人口, 出生率, 気温, 切片）:', w_src.round(3))
res_s = Xs @ w_src - ys
res_t = Xt @ w_src - yt
print(f'source 内の RMSE: {np.sqrt((res_s**2).mean()):.3f}')
print(f'2023 年度にそのまま当てた RMSE: {np.sqrt((res_t**2).mean()):.3f}  平均の予測誤差: {res_t.mean():+.3f}')
print(f'2023 年度の死亡率 平均 {yt.mean():.2f}  予測の平均 {(Xt @ w_src).mean():.2f}')

# tl1.py の続き（Xs, Xt, ys, yt, w_src を使う）
from sklearn.linear_model import Ridge

def rmse(w, i):
    return np.sqrt(((Xt[i] @ w - yt[i]) ** 2).mean())

rng = np.random.default_rng(0)
print(' k | source のまま | 切片だけ補正 | source へ縮小 | target だけ')
for k in [3, 5, 8, 12, 20, 30]:
    r = {'src': [], 'bias': [], 'shrink': [], 'tgt': []}
    for _ in range(300):
        lab = rng.choice(47, k, replace=False)
        val = np.setdiff1d(np.arange(47), lab)
        r['src'].append(rmse(w_src, val))
        w_b = w_src.copy(); w_b[-1] += (yt[lab] - Xt[lab] @ w_src).mean()   # 切片だけ動かす
        r['bias'].append(rmse(w_b, val))
        d = Ridge(alpha=5.0).fit(Xt[lab][:, :4], yt[lab] - Xt[lab] @ w_src)   # 残差を小さな修正で学習
        w_s = w_src + np.r_[d.coef_, d.intercept_]
        r['shrink'].append(rmse(w_s, val))
        w_t = np.linalg.lstsq(Xt[lab], yt[lab], rcond=None)[0]              # target の k 県だけで最小二乗
        r['tgt'].append(rmse(w_t, val))
    print(f"{k:2d} | {np.median(r['src']):12.3f} | {np.median(r['bias']):12.3f} | {np.median(r['shrink']):12.3f} | {np.median(r['tgt']):10.3f}")

# tl1.py の続き（df, F, mu, sd, Xt, yt を使う）
print('source 年度 | 2023 への RMSE | 平均の予測誤差 | その年度の死亡率平均')
for yr in range(2012, 2023):
    s = df[df['SSDSE-B-2026'] == yr]
    X1 = np.c_[((s[F] - mu) / sd).to_numpy(), np.ones(len(s))]
    w = np.linalg.lstsq(X1, s['dr'].to_numpy(), rcond=None)[0]
    e = Xt @ w - yt
    print(f'{yr} | {np.sqrt((e**2).mean()):6.3f} | {e.mean():+6.3f} | {s["dr"].mean():6.2f}')

# tl1.py の続き（Xt, yt, w_src を使う）
rng = np.random.default_rng(0)
lab = np.sort(rng.choice(47, 12, replace=False))   # ラベルを使える 12 県
val = np.setdiff1d(np.arange(47), lab)             # 検証 35 県
print('ラベルあり 12 県:', ' '.join(tgt['Prefecture'].to_numpy()[lab]))

def mse(w, i):
    return ((Xt[i] @ w - yt[i]) ** 2).mean()

def fine_tune(w0, epochs, lr=0.05):
    w, hist = w0.copy(), []
    for _ in range(epochs + 1):
        hist.append((mse(w, lab), mse(w, val)))
        w = w - lr * 2 * Xt[lab].T @ (Xt[lab] @ w - yt[lab]) / len(lab)   # 12 県で 1 回更新 = 1 エポック
    return np.array(hist)

ft = fine_tune(w_src, 300)                          # source の重みから始める
sc = fine_tune(np.r_[np.zeros(4), yt[lab].mean()], 300)   # 重み 0 から始める（スクラッチ）
for e in [0, 5, 20, 50, 100, 300]:
    print(f'エポック {e:3d}: fine-tune 学習 {ft[e,0]:.3f} 検証 {ft[e,1]:.3f} | スクラッチ 検証 {sc[e,1]:.3f}')
b = ft[:, 1].argmin()
print(f'fine-tune の検証 MSE 最小: エポック {b} で {ft[b,1]:.3f}')

# tl1.py・tl4.py の続き（lab, val を使う）
big = src.groupby('Prefecture')['A1101'].mean().nlargest(10).index      # 2012〜2019 年度の人口上位 10 都道府県
for name, s in [('全 47 県', src), ('人口上位 10 県だけ', src[src['Prefecture'].isin(big)])]:
    X1 = np.c_[((s[F] - mu) / sd).to_numpy(), np.ones(len(s))]
    w = np.linalg.lstsq(X1, s['dr'].to_numpy(), rcond=None)[0]
    w_b = w.copy(); w_b[-1] += (yt[lab] - Xt[lab] @ w).mean()          # 12 県で切片だけ補正
    print(f'{name:10s} 行数 {len(s):3d}  高齢化率の重み {w[0]:+.3f}  検証 RMSE: そのまま {np.sqrt(mse(w, val)):.3f} / 切片補正 {np.sqrt(mse(w_b, val)):.3f}')
w_t = np.linalg.lstsq(Xt[lab], yt[lab], rcond=None)[0]
print(f'target 12 県だけの最小二乗: 検証 RMSE {np.sqrt(mse(w_t, val)):.3f}')
print('source の高齢化率の範囲:', round(src[src['Prefecture'].isin(big)]['age'].min(), 1), '〜', round(src[src['Prefecture'].isin(big)]['age'].max(), 1), '%  2023 年度 47 県:', round(tgt['age'].min(), 1), '〜', round(tgt['age'].max(), 1), '%')

# tl1.py の続き（src, tgt, Xs, Xt, ys, yt を使う）
from scipy import stats
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import cross_val_score

print('特徴量      source 平均 | 2023 平均 | 差/source の標準偏差 | KS 統計量 (p 値)')
for c, name in zip(F, ['高齢化率', 'log10人口', '出生率', '年平均気温']):
    ks = stats.ks_2samp(src[c], tgt[c])
    print(f'{name:8s} {src[c].mean():9.3f} | {tgt[c].mean():8.3f} | {(tgt[c].mean() - src[c].mean()) / src[c].std():+6.2f} | {ks.statistic:.3f} ({ks.pvalue:.4f})')

# source の行か target の行かを当てる分類器（adversarial validation）
Xa = np.r_[Xs[:, :4], Xt[:, :4]]
ya = np.r_[np.zeros(len(Xs)), np.ones(len(Xt))]
auc = cross_val_score(LogisticRegression(max_iter=1000), Xa, ya, cv=5, scoring='roc_auc')
print(f'source と target を見分ける AUC（5 分割）: {auc.mean():.3f}')

# tl1.py・tl4.py の続き（fine_tune, w_src, lab, val, mse を使う）
for lr in [0.005, 0.05, 0.1, 0.2]:
    h = fine_tune(w_src, 300, lr=lr)
    b = h[:, 1].argmin()
    print(f'学習率 {lr:5.3f}: 検証 MSE 最小 {h[b,1]:.3f}（エポック {b:3d}）  300 エポック後 {h[300,1]:.3g}')
lam = np.linalg.eigvalsh(2 * Xt[lab].T @ Xt[lab] / len(lab)).max()   # 損失の曲がり具合の最大値
print(f'発散しない学習率の上限 2 / λmax = 2 / {lam:.2f} = {2 / lam:.3f}')
# source の重みがどれだけ動いたか（学習率 0.05、エポック 108 と 300）
w = w_src.copy()
for e in range(300):
    w = w - 0.05 * 2 * Xt[lab].T @ (Xt[lab] @ w - yt[lab]) / len(lab)
    if e + 1 in (108, 300):
        print(f'エポック {e+1}: 重み {w.round(3)}  source からの変化 {np.abs(w - w_src).round(3)}')

# tl1.py の続き（Xt, yt, w_src を使う）
names = ['切片だけ更新（他は凍結）', '切片＋高齢化率を更新', '高齢化率・log人口・切片を更新', '5 つ全部を更新（最小二乗）']
free = [[4], [0, 4], [0, 1, 4], [0, 1, 2, 3, 4]]     # 更新する重みの位置（0=高齢化率 … 4=切片）
rng = np.random.default_rng(1)
res = {n: [] for n in names}
for _ in range(300):
    lab = rng.choice(47, 12, replace=False)
    val = np.setdiff1d(np.arange(47), lab)
    for n, f in zip(names, free):
        w = w_src.copy()
        fixed = [j for j in range(5) if j not in f]
        r = yt[lab] - Xt[lab][:, fixed] @ w[fixed]          # 凍結した重みの寄与を引いた残り
        w[f] = np.linalg.lstsq(Xt[lab][:, f], r, rcond=None)[0]
        res[n].append(np.sqrt(((Xt[val] @ w - yt[val]) ** 2).mean()))
for n in names:
    v = np.array(res[n])
    print(f'{n:20s} 検証 RMSE 中央値 {np.median(v):.3f}  90% 点 {np.quantile(v, 0.9):.3f}')

# tl1.py の続き（df, src, tgt, yt を使う）
from sklearn.decomposition import PCA
from sklearn.linear_model import LinearRegression

P_ = df['A1101']
df2 = pd.DataFrame({'年少人口比': df['A1301'] / P_, '高齢化率': df['age'], '出生率': df['A4103'],
    '転入率': df['A5101'] / P_, '転出率': df['A5102'] / P_, '婚姻率': df['A9101'] / P_, '離婚率': df['A9201'] / P_,
    '気温': df['B4101'], '降水量': df['B4109'], '病院/人口': df['I510120'] / P_, '診療所/人口': df['I5102'] / P_,
    '保育所/人口': df['J2503'] / P_, 'log人口': df['lp']})
Z_src = df2.loc[src.index]; Z_tgt = df2.loc[tgt.index]
m2, s2 = Z_src.mean(), Z_src.std()
Zs = ((Z_src - m2) / s2).to_numpy(); Zt = ((Z_tgt - m2) / s2).to_numpy()
pca = PCA(3).fit(Zs)                                   # 特徴の「まとめ方」を source 376 行で学習して凍結
print('source で学習した PCA の累積寄与率:', pca.explained_variance_ratio_.cumsum().round(3))
rng = np.random.default_rng(2)
print(' k | 13 列をそのまま最小二乗 | target k 県で PCA | source の PCA を凍結')
for k in [5, 8, 12, 20]:
    r = [[], [], []]
    for _ in range(300):
        lab = rng.choice(47, k, replace=False); val = np.setdiff1d(np.arange(47), lab)
        w = np.linalg.lstsq(np.c_[Zt[lab], np.ones(k)], yt[lab], rcond=None)[0]
        r[0].append(np.sqrt(((np.c_[Zt[val], np.ones(len(val))] @ w - yt[val]) ** 2).mean()))
        p_t = PCA(min(3, k - 1)).fit(Zt[lab])
        h = LinearRegression().fit(p_t.transform(Zt[lab]), yt[lab])
        r[1].append(np.sqrt(((h.predict(p_t.transform(Zt[val])) - yt[val]) ** 2).mean()))
        h = LinearRegression().fit(pca.transform(Zt[lab]), yt[lab])
        r[2].append(np.sqrt(((h.predict(pca.transform(Zt[val])) - yt[val]) ** 2).mean()))
    print(f'{k:2d} | {np.median(r[0]):22.3f} | {np.median(r[1]):16.3f} | {np.median(r[2]):18.3f}')

# tl1.py の続き（Xt, yt, w_src を使う）
from sklearn.metrics import roc_auc_score
rng = np.random.default_rng(0)
lab = np.sort(rng.choice(47, 12, replace=False))     # tl4.py と同じ 12 県
val = np.setdiff1d(np.arange(47), lab)
thr = np.median(yt)                                   # 2023 年度 47 県の死亡率の中央値
yb = (yt[val] >= thr).astype(int)
w = w_src.copy()
for _ in range(100):                                  # 学習率 0.05 で 100 エポック fine-tune
    w = w - 0.05 * 2 * Xt[lab].T @ (Xt[lab] @ w - yt[lab]) / len(lab)
print(f'しきい値（中央値）: {thr:.2f}‰  検証 35 県のうち中央値以上: {yb.sum()} 県')
for name, ww in [('source のまま', w_src), ('fine-tune 100 エポック', w)]:
    p = Xt[val] @ ww
    print(f'{name:14s} AUC {roc_auc_score(yb, p):.3f}  正解率 {((p >= thr) == yb).mean():.3f}  「以上」と予測 {(p >= thr).sum()} 県  予測の平均 {p.mean():.2f}')

# tl1.py・tl2.py の続き（tgt, Xt, yt, w_src, rmse を使う）
east = np.where(tgt['Code'].str[1:3].astype(int).to_numpy() <= 15)[0]   # 北海道・東北・関東・新潟（地域コード 01〜15）
rng = np.random.default_rng(3)
for name, pool in [('47 県から無作為', np.arange(47)), ('東日本 15 都道県からだけ', east)]:
    r_b, r_t = [], []
    for _ in range(300):
        lab = rng.choice(pool, 12, replace=False); val = np.setdiff1d(np.arange(47), lab)
        w_b = w_src.copy(); w_b[-1] += (yt[lab] - Xt[lab] @ w_src).mean()
        w_t = np.linalg.lstsq(Xt[lab], yt[lab], rcond=None)[0]
        r_b.append(rmse(w_b, val)); r_t.append(rmse(w_t, val))
    print(f'{name:14s}: 切片だけ補正 {np.median(r_b):.3f}  target 12 県だけ {np.median(r_t):.3f}')
print('東日本 15 都道県の高齢化率:', round(tgt['age'].to_numpy()[east].min(), 1), '〜', round(tgt['age'].to_numpy()[east].max(), 1), '%')
west = np.setdiff1d(np.arange(47), east)
res = yt - Xt @ w_src
print(f'source モデルの過小予測（実際 − 予測）の平均: 東日本 {res[east].mean():+.3f}  それ以外 32 府県 {res[west].mean():+.3f}')

# tl1.py・tl2.py の続き（Xt, yt, w_src, rmse, Ridge を使う）
rng = np.random.default_rng(4)
splits = [rng.choice(47, 12, replace=False) for _ in range(300)]
print('  α     | 検証 RMSE 中央値 | 高齢化率の重みの平均')
for alpha in [0.01, 0.1, 1, 5, 20, 100, 1000]:
    r, wa = [], []
    for lab in splits:
        val = np.setdiff1d(np.arange(47), lab)
        d = Ridge(alpha=alpha).fit(Xt[lab][:, :4], yt[lab] - Xt[lab] @ w_src)   # source からの修正量を学習
        w = w_src + np.r_[d.coef_, d.intercept_]
        r.append(rmse(w, val)); wa.append(w[0])
    print(f'{alpha:7g} | {np.median(r):16.3f} | {np.mean(wa):.3f}')

# tl1.py の続き（src, tgt, F, mu, sd, Xt, yt を使う）
rng = np.random.default_rng(5)
codes = tgt['Code'].to_numpy()
res = {'source に評価の 12 県を含む': [], 'source から評価の 12 県を抜く': []}
for _ in range(200):
    val = rng.choice(47, 12, replace=False)                        # 評価する 12 県
    lab = np.setdiff1d(np.arange(47), val)                          # 2023 年度のラベルは残り 35 県
    for name, keep in [('source に評価の 12 県を含む', np.ones(len(src), bool)),
                       ('source から評価の 12 県を抜く', ~src['Code'].isin(codes[val]).to_numpy())]:
        s = src[keep]
        X1 = np.c_[((s[F] - mu) / sd).to_numpy(), np.ones(len(s))]
        w = np.linalg.lstsq(X1, s['dr'].to_numpy(), rcond=None)[0]
        w[-1] += (yt[lab] - Xt[lab] @ w).mean()                     # 35 県で切片を補正
        res[name].append(np.sqrt(((Xt[val] @ w - yt[val]) ** 2).mean()))
for name, v in res.items():
    print(f'{name}: 検証 RMSE 中央値 {np.median(v):.3f}')
print('source の行数: 含む', len(src), '行 / 抜く', len(src) - 12 * 8, '行')

# tl1.py の続き（src, Xs, ys, Xt, yt を使う）
rng = np.random.default_rng(6)
pref_s = src['Code'].to_numpy(); up = np.unique(pref_s)
B_s, B_t = [], []
for _ in range(1000):
    pick = rng.choice(up, 47, replace=True)                       # source は県ごと（8 年度まとめて）復元抽出
    idx = np.concatenate([np.where(pref_s == c)[0] for c in pick])
    B_s.append(np.linalg.lstsq(Xs[idx], ys[idx], rcond=None)[0])
    j = rng.choice(47, 47, replace=True)                           # target は 47 県を復元抽出
    B_t.append(np.linalg.lstsq(Xt[j], yt[j], rcond=None)[0])
B_s, B_t = np.array(B_s), np.array(B_t)
for i, n in enumerate(['高齢化率', 'log人口', '出生率', '気温', '切片']):
    ls, hs = np.quantile(B_s[:, i], [0.025, 0.975]); lt, ht = np.quantile(B_t[:, i], [0.025, 0.975])
    print(f'{n:6s} source [{ls:+.3f}, {hs:+.3f}]  2023 年度 [{lt:+.3f}, {ht:+.3f}]  重なり {"あり" if max(ls, lt) <= min(hs, ht) else "なし"}')