y = 10
pred = 5
lr = 0.5
for i in range(4):
    print(f"段 {i}: pred={pred}, 残差={y-pred}")
    pred = pred + lr * (y - pred)

import pandas as pd
import xgboost as xgb
from sklearn.model_selection import train_test_split

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

feats = ['A1301', 'A4101', 'A4200', 'A5101', 'A9101', 'E1101', 'C3301', 'F3101']
X = last[feats]
y = pd.qcut(last['A1101'], q=3, labels=[0, 1, 2]).astype(int)

Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.2, random_state=0, stratify=y)

model = xgb.XGBClassifier(
    n_estimators=200, max_depth=4, learning_rate=0.05,
    objective='multi:softprob', num_class=3, eval_metric='mlogloss',
)
model.fit(Xtr, ytr)
print(f'Train acc: {model.score(Xtr,ytr):.4f}')
print(f'Test  acc: {model.score(Xte,yte):.4f}')

import pandas as pd
import xgboost as xgb
from sklearn.model_selection import train_test_split

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
last = df.query('`SSDSE-B-2026` == `SSDSE-B-2026`.max()').copy()
feats = ['A1301', 'A4101', 'A4200', 'A5101', 'A9101', 'E1101', 'C3301', 'F3101']
X, y = last[feats], pd.qcut(last['A1101'], 3, labels=[0,1,2]).astype(int)
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.3, random_state=0, stratify=y)

model = xgb.XGBClassifier(
    n_estimators=500, max_depth=3, learning_rate=0.05,
    objective='multi:softprob', num_class=3, eval_metric='mlogloss',
    early_stopping_rounds=20,
)
model.fit(Xtr, ytr, eval_set=[(Xtr,ytr), (Xte,yte)], verbose=False)
print(f'Best iter: {model.best_iteration}')
print(f'Train acc: {model.score(Xtr,ytr):.4f}  Test acc: {model.score(Xte,yte):.4f}')

import pandas as pd
import numpy as np
import xgboost as xgb
from sklearn.model_selection import cross_val_score, StratifiedKFold

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
last = df.query('`SSDSE-B-2026` == `SSDSE-B-2026`.max()').copy()
feats = ['A1301', 'A4101', 'A4200', 'A5101', 'A9101', 'E1101', 'C3301', 'F3101']
X, y = last[feats], pd.qcut(last['A1101'], 3, labels=[0,1,2]).astype(int)

model = xgb.XGBClassifier(n_estimators=100, max_depth=3, learning_rate=0.05,
                          objective='multi:softprob', num_class=3, eval_metric='mlogloss')
cv     = StratifiedKFold(n_splits=5, shuffle=True, random_state=0)
scores = cross_val_score(model, X, y, cv=cv, scoring='accuracy')
print(f'CV accuracy: {scores.mean():.4f} ± {scores.std():.4f}')
print(f'各 fold: {scores.round(3)}')

import pandas as pd
import xgboost as xgb
from sklearn.model_selection import GridSearchCV

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
last = df.query('`SSDSE-B-2026` == `SSDSE-B-2026`.max()').copy()
feats = ['A1301', 'A4101', 'A4200', 'A5101', 'A9101', 'E1101', 'C3301', 'F3101']
X, y = last[feats], pd.qcut(last['A1101'], 3, labels=[0,1,2]).astype(int)

grid = {
    'n_estimators':  [50, 100, 200],
    'max_depth':     [2, 3, 4, 6],
    'learning_rate': [0.01, 0.05, 0.1],
}
model = xgb.XGBClassifier(objective='multi:softprob', num_class=3, eval_metric='mlogloss', random_state=0)
gs = GridSearchCV(model, grid, cv=3, scoring='accuracy', n_jobs=-1)
gs.fit(X, y)
print(f'Best params: {gs.best_params_}')
print(f'Best score : {gs.best_score_:.4f}')

import pandas as pd
import xgboost as xgb
from sklearn.inspection import permutation_importance

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
last = df.query('`SSDSE-B-2026` == `SSDSE-B-2026`.max()').copy()
feats = ['A1301', 'A4101', 'A4200', 'A5101', 'A9101', 'E1101', 'C3301', 'F3101']
X, y = last[feats], pd.qcut(last['A1101'], 3, labels=[0,1,2]).astype(int)

model = xgb.XGBClassifier(n_estimators=100, max_depth=3, learning_rate=0.1,
                          objective='multi:softprob', num_class=3, eval_metric='mlogloss').fit(X,y)

gain_imp = pd.Series(model.feature_importances_, index=feats).sort_values(ascending=False)
perm     = permutation_importance(model, X, y, n_repeats=10, random_state=0)
perm_imp = pd.Series(perm.importances_mean, index=feats).sort_values(ascending=False)

print('feature_importances_ (gain):')
print(gain_imp.to_string())
print()
print('permutation_importance:')
print(perm_imp.to_string())

import pandas as pd
import numpy as np
import xgboost as xgb
import shap

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
last = df.query('`SSDSE-B-2026` == `SSDSE-B-2026`.max()').reset_index(drop=True)
feats = ['A1301', 'A4101', 'A4200', 'A5101', 'A9101', 'E1101', 'C3301', 'F3101']
X, y = last[feats], pd.qcut(last['A1101'], 3, labels=[0,1,2]).astype(int)

model     = xgb.XGBClassifier(n_estimators=100, max_depth=3, learning_rate=0.1,
                              objective='multi:softprob', num_class=3, eval_metric='mlogloss', random_state=0).fit(X,y)
explainer = shap.TreeExplainer(model)
# shap 0.4x 以降、多クラスの shap_values は list ではなく (サンプル, 特徴量, クラス) の 3 次元配列
sv        = np.asarray(explainer.shap_values(X))
print('sv.shape =', sv.shape)

# skiprows=[1] で読むと列名は英語コードなので、県名の列は 'Prefecture'
tokyo_pos = last.index[last['Prefecture'] == '東京都'][0]
print('東京の クラス 2 (大県) SHAP 寄与:')
for i, f in enumerate(feats):
    print(f'  {f:8s} SHAP={sv[tokyo_pos, i, 2]:+.3f}  値={X.iloc[tokyo_pos, i]:.1f}')

import pandas as pd
import xgboost as xgb
from sklearn.model_selection import train_test_split

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
last = df.query('`SSDSE-B-2026` == `SSDSE-B-2026`.max()').copy()
feats = ['A1301', 'A4101', 'A4200', 'A5101', 'A9101', 'E1101', 'C3301', 'F3101']
X, y = last[feats], pd.qcut(last['A1101'], 3, labels=[0,1,2]).astype(int)
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.3, stratify=y, random_state=0)

model = xgb.XGBClassifier(n_estimators=150, max_depth=3, learning_rate=0.1,
                          objective='multi:softprob', num_class=3, eval_metric='mlogloss')
model.fit(Xtr, ytr, eval_set=[(Xtr,ytr), (Xte,yte)], verbose=False)

hist = model.evals_result_
print(f'最初 5 round:')
for i in range(5):
    print(f'  iter {i:3d}  train={hist["validation_0"]["mlogloss"][i]:.4f}  test={hist["validation_1"]["mlogloss"][i]:.4f}')
print(f'最終 round 150:  train={hist["validation_0"]["mlogloss"][-1]:.4f}  test={hist["validation_1"]["mlogloss"][-1]:.4f}')

import pandas as pd, numpy as np, xgboost as xgb
from sklearn.model_selection import KFold, cross_val_score
df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
df = df[df['SSDSE-B-2026'] == 2023].reset_index(drop=True)   # 最新年度(2023)の47都道府県
cols = ['A1101', 'A1303', 'I510120', 'J2503']   # 総人口/65歳以上人口/一般病院数/保育所等数
X = df[cols].astype(float)
y = df['A4103'].astype(float)                        # 目的変数: 合計特殊出生率
best = None
for d in [3, 5, 7, 9]:
    for e in [0.05, 0.1, 0.2]:
        m = xgb.XGBRegressor(max_depth=d, learning_rate=e, n_estimators=400,
                             subsample=0.8, colsample_bytree=0.8, random_state=0)
        cv = cross_val_score(m, X, y, cv=KFold(5, shuffle=True, random_state=0), scoring='r2').mean()
        if best is None or cv > best[0]:
            best = (cv, d, e)
print(f'best R2 = {best[0]:.3f} @ max_depth={best[1]}, eta={best[2]}')

import xgboost as xgb, pandas as pd, shap
from sklearn.model_selection import train_test_split
df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
df = df[df['SSDSE-B-2026'] == 2023].reset_index(drop=True)   # 最新年度(2023)の47都道府県
cols = ['A1101', 'A1303', 'I510120', 'J2503']   # 総人口/65歳以上人口/一般病院数/保育所等数
X = df[cols].astype(float)
y = df['A4103'].astype(float)                        # 合計特殊出生率
Xtr, Xva, ytr, yva = train_test_split(X, y, test_size=0.2, random_state=42)
model = xgb.XGBRegressor(n_estimators=2000, max_depth=5, learning_rate=0.05,
                         subsample=0.8, colsample_bytree=0.8,
                         early_stopping_rounds=30, eval_metric='rmse', random_state=0)
model.fit(Xtr, ytr, eval_set=[(Xva, yva)], verbose=False)
print(f'最適な木の本数: {model.best_iteration}')
print(f'検証 RMSE       : {model.best_score:.3f}')
explainer = shap.TreeExplainer(model)
sv = explainer.shap_values(Xva)
print('SHAP 平均 |寄与|:')
for c, v in zip(X.columns, abs(sv).mean(axis=0)):
    print(f'  {c}: {v:.3f}')

import os
os.makedirs('figures', exist_ok=True)  # 保存先のフォルダを作っておく

import pandas as pd
import numpy as np
import xgboost as xgb
import matplotlib.pyplot as plt

# 先頭列名は "SSDSE-B-2026"（=年度）、2行目の日本語見出しは skiprows=[1] で除去
df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
df = df[df['SSDSE-B-2026'] == 2023].reset_index(drop=True)   # 最新年度(2023)の47都道府県だけ残す
X = df[['A1101']].astype(float).values     # 説明変数: A1101 総人口
y = df['A4101'].astype(float).values       # 目的変数: A4101 出生数

model = xgb.XGBRegressor(n_estimators=200, max_depth=3, learning_rate=0.1,
                         objective='reg:squarederror', random_state=0)
model.fit(X, y)
y_pred = model.predict(X)

order = np.argsort(X.ravel())
plt.figure(figsize=(7, 5))
plt.scatter(X, y, color='#1976D2', label='actual')
plt.plot(X[order], y_pred[order], color='#E64A19', linewidth=2, label='XGBoost predicted')
plt.xlabel('total population'); plt.ylabel('births')
plt.legend(); plt.tight_layout(); plt.savefig('figures/xgb_scatter_fit.png', dpi=120)

# 当てはまりと残差の大きい県を数値で確かめる（学習データそのものへの当てはまり）
resid = y - y_pred
print(f'RMSE (train) = {np.sqrt(np.mean(resid**2)):.1f} 人')
print(f'R^2 (train)  = {1 - np.sum(resid**2) / np.sum((y - y.mean())**2):.4f}')
print(f'最大残差県 = {df["Prefecture"][resid.argmax()]} ({resid.max():+.0f} 人)')
print(f'最小残差県 = {df["Prefecture"][resid.argmin()]} ({resid.min():+.0f} 人)')

import os
os.makedirs('figures', exist_ok=True)  # 保存先のフォルダを作っておく

import pandas as pd, numpy as np, xgboost as xgb, matplotlib.pyplot as plt
from sklearn.model_selection import cross_val_predict, KFold
df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
df = df[df['SSDSE-B-2026'] == 2023].reset_index(drop=True)   # 最新年度(2023)の47都道府県
cols = ['A1101', 'A1303', 'I510120', 'J2503']   # 総人口/65歳以上人口/一般病院数/保育所等数
target = 'A4103'                                     # 合計特殊出生率
X = df[cols].values
y = df[target].astype(float).values
model = xgb.XGBRegressor(n_estimators=300, max_depth=4, learning_rate=0.05,
                         subsample=0.8, colsample_bytree=0.8, random_state=0)
# 分布縮小はホールドアウト予測で確認する(学習データ内予測は過学習で縮まない)
y_pred = cross_val_predict(model, X, y, cv=KFold(5, shuffle=True, random_state=0))
plt.figure(figsize=(7, 5))
plt.hist(y, bins=12, alpha=0.6, color='#1976D2', label='actual')
plt.hist(y_pred, bins=12, alpha=0.6, color='#E64A19', label='predicted')
plt.xlabel('total fertility rate'); plt.ylabel('count')
plt.legend(); plt.tight_layout(); plt.savefig('figures/xgb_hist.png', dpi=120)

print(f'真値分布 : 平均 {y.mean():.3f}、 標準偏差 {y.std():.3f}、 最大 {y.max():.3f}、 最小 {y.min():.3f}')
print(f'予測分布 : 平均 {y_pred.mean():.3f}、 標準偏差 {y_pred.std():.3f}、 最大 {y_pred.max():.3f}、 最小 {y_pred.min():.3f}')
print(f'分布縮小率 : 標準偏差ベースで 約 {(1 - y_pred.std() / y.std()) * 100:.0f}% 圧縮')

import os
os.makedirs('figures', exist_ok=True)  # 保存先のフォルダを作っておく

import pandas as pd, numpy as np, xgboost as xgb, matplotlib.pyplot as plt
from sklearn.model_selection import cross_val_predict, KFold
df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
df = df[df['SSDSE-B-2026'] == 2023].reset_index(drop=True)   # 最新年度(2023)の47都道府県
block_map = {'北海道':'北海道','青森県':'東北','岩手県':'東北','宮城県':'東北','秋田県':'東北','山形県':'東北','福島県':'東北',
             '茨城県':'関東','栃木県':'関東','群馬県':'関東','埼玉県':'関東','千葉県':'関東','東京都':'関東','神奈川県':'関東',
             '新潟県':'中部','富山県':'中部','石川県':'中部','福井県':'中部','山梨県':'中部','長野県':'中部','岐阜県':'中部','静岡県':'中部','愛知県':'中部',
             '三重県':'近畿','滋賀県':'近畿','京都府':'近畿','大阪府':'近畿','兵庫県':'近畿','奈良県':'近畿','和歌山県':'近畿',
             '鳥取県':'中国','島根県':'中国','岡山県':'中国','広島県':'中国','山口県':'中国',
             '徳島県':'四国','香川県':'四国','愛媛県':'四国','高知県':'四国',
             '福岡県':'九州沖縄','佐賀県':'九州沖縄','長崎県':'九州沖縄','熊本県':'九州沖縄','大分県':'九州沖縄','宮崎県':'九州沖縄','鹿児島県':'九州沖縄','沖縄県':'九州沖縄'}
df['block'] = df['Prefecture'].map(block_map)
cols = ['A1101', 'A1303', 'I510120', 'J2503']   # 総人口/65歳以上人口/一般病院数/保育所等数
target = 'A4103'
model = xgb.XGBRegressor(n_estimators=300, max_depth=4, learning_rate=0.05, random_state=0)
# ブロック別の系統誤差はホールドアウト残差で評価する
df['resid'] = df[target].values - cross_val_predict(model, df[cols].values,
              df[target].astype(float).values, cv=KFold(5, shuffle=True, random_state=0))
order = ['北海道','東北','関東','中部','近畿','中国','四国','九州沖縄']
groups = [df.loc[df['block']==b, 'resid'].values for b in order]
plt.figure(figsize=(8, 5))
plt.boxplot(groups, tick_labels=order)
plt.ylabel('residual (true - pred)')
plt.tight_layout(); plt.savefig('figures/xgb_box.png', dpi=120)

# ブロックごとの残差中央値と IQR（真値 − 予測。正なら過小予測）
q = df.groupby('block')['resid'].describe()[['count', '25%', '50%', '75%']]
q['IQR'] = q['75%'] - q['25%']
print(q.loc[order, ['count', '50%', 'IQR']].round(3))

import pandas as pd, numpy as np, xgboost as xgb
from sklearn.model_selection import KFold
from sklearn.metrics import r2_score
df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
df = df[df['SSDSE-B-2026'] == 2023].reset_index(drop=True)   # 最新年度(2023)の47都道府県
cols = ['A1101', 'A1303', 'I510120', 'J2503']   # 総人口/65歳以上人口/一般病院数/保育所等数
target = 'A4103'
X, y = df[cols].values, df[target].astype(float).values
kf = KFold(n_splits=5, shuffle=True, random_state=0)
train_scores, valid_scores = [], []
for tr, va in kf.split(X):
    m = xgb.XGBRegressor(n_estimators=300, max_depth=4, learning_rate=0.05, random_state=0)
    m.fit(X[tr], y[tr])
    train_scores.append(r2_score(y[tr], m.predict(X[tr])))
    valid_scores.append(r2_score(y[va], m.predict(X[va])))
print('train R^2 mean =', np.mean(train_scores))
print('valid R^2 mean =', np.mean(valid_scores))

import os
os.makedirs('figures', exist_ok=True)  # 保存先のフォルダを作っておく

import pandas as pd, numpy as np, xgboost as xgb, matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.metrics import r2_score
df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
df = df[df['SSDSE-B-2026'] == 2023].reset_index(drop=True)   # 最新年度(2023)の47都道府県
cols = ['A1101', 'A1303', 'I510120', 'J2503']   # 総人口/65歳以上人口/一般病院数/保育所等数
target = 'A4103'
X, y = df[cols].values, df[target].astype(float).values
X_tr_full, X_te, y_tr_full, y_te = train_test_split(X, y, test_size=0.3, random_state=0)
sizes = [10, 15, 20, 25, 30, len(X_tr_full)]
tr_scores, te_scores = [], []
for s in sizes:
    m = xgb.XGBRegressor(n_estimators=200, max_depth=3, learning_rate=0.05, random_state=0)
    m.fit(X_tr_full[:s], y_tr_full[:s])
    tr_scores.append(r2_score(y_tr_full[:s], m.predict(X_tr_full[:s])))
    te_scores.append(r2_score(y_te, m.predict(X_te)))
plt.plot(sizes, tr_scores, 'o-', label='train R^2', color='#1976D2')
plt.plot(sizes, te_scores, 's-', label='test R^2',  color='#E64A19')
plt.xlabel('train size'); plt.ylabel('R^2'); plt.legend(); plt.tight_layout()
plt.savefig('figures/xgb_learning_curve.png', dpi=120)

print('size  train R^2  test R^2')
for s, a, b in zip(sizes, tr_scores, te_scores):
    print(f'{s:4d}  {a:9.3f}  {b:8.2f}')

import pandas as pd, numpy as np, xgboost as xgb, shap
df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
df = df[df['SSDSE-B-2026'] == 2023].reset_index(drop=True)   # 最新年度(2023)の47都道府県
cols = ['A1101', 'A1303', 'I510120', 'J2503']   # 総人口/65歳以上人口/一般病院数/保育所等数
target = 'A4103'
X, y = df[cols].values, df[target].astype(float).values
model = xgb.XGBRegressor(n_estimators=300, max_depth=4, learning_rate=0.05, random_state=0)
model.fit(X, y)
explainer = shap.TreeExplainer(model)
shap_vals = explainer.shap_values(X)

# 沖縄県の SHAP 内訳
i = df.index[df['Prefecture'] == '沖縄県'][0]
print('high-contribution variables (Okinawa):')
for col, v in sorted(zip(cols, shap_vals[i]), key=lambda x: -abs(x[1])):
    sign = '+' if v >= 0 else '-'
    print(f'  {col:8s} {sign}{abs(v):.3f}')
print('base value (mean of train target) =', explainer.expected_value)
print('predicted =', model.predict(X[i:i+1])[0])
print('actual    =', y[i])

import pandas as pd, xgboost as xgb
from sklearn.model_selection import train_test_split
df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
df = df[df['SSDSE-B-2026'] == 2023].reset_index(drop=True)   # 最新年度(2023)の47都道府県
cols = ['A1101', 'A1303', 'I510120', 'J2503']   # 総人口/65歳以上人口/一般病院数/保育所等数
target = 'A4103'
X_tr, X_va, y_tr, y_va = train_test_split(df[cols].values, df[target].astype(float).values,
                                          test_size=0.3, random_state=0)
model = xgb.XGBRegressor(n_estimators=2000, max_depth=3, learning_rate=0.02,
                         early_stopping_rounds=50, random_state=0)
model.fit(X_tr, y_tr, eval_set=[(X_va, y_va)], verbose=False)
print('best iteration =', model.best_iteration)
print('best score     =', model.best_score)

import os, hashlib
import numpy as np, pandas as pd, xgboost as xgb

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
d = df[df['SSDSE-B-2026'] == 2023]                      # 2023 年度の 47 都道府県
X = d[['A1101', 'A1303', 'I510120', 'J2503']].values        # 総人口・65歳以上人口・一般病院数・保育所等数
y = d['A4103'].values                                     # 合計特殊出生率

model = xgb.XGBRegressor(n_estimators=300, max_depth=4, learning_rate=0.05, random_state=0)
model.fit(X, y)
os.makedirs('models', exist_ok=True)
model.save_model('models/xgb_ssdse_b_2026.json')

# 再読込して、保存前と同じ予測が出るか確かめる
loaded = xgb.XGBRegressor()
loaded.load_model('models/xgb_ssdse_b_2026.json')
print('予測の最大差 =', np.abs(model.predict(X) - loaded.predict(X)).max())

# 配布時にハッシュを明記
h = hashlib.sha256(open('models/xgb_ssdse_b_2026.json', 'rb').read()).hexdigest()
print('sha256 =', h[:16], '...')
print(f"file size = {os.path.getsize('models/xgb_ssdse_b_2026.json') / 1024:.1f} KB")

import pandas as pd, numpy as np, xgboost as xgb
from sklearn.model_selection import KFold
from sklearn.metrics import r2_score

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', skiprows=1, encoding='cp932')
df = df[df['年度'] == 2023].reset_index(drop=True)   # 2023 年度の 47 都道府県

# 県民幸福度スコア (簡易版): SSDSE-B にある率で、正の要素 - 負の要素を z-score 化して平均
df['高校進学率'] = df['高等学校卒業者のうち進学者数'] / df['高等学校卒業者数']
df['診療所密度'] = df['一般診療所数'] / df['総人口'] * 1e5      # 人口 10 万人あたり
df['離婚率'] = df['離婚件数'] / df['総人口'] * 1000              # 人口千人あたり
df['待機児童率'] = df['保育所等利用待機児童数'] / df['15歳未満人口'] * 1000
positives = ['合計特殊出生率', '高校進学率', '診療所密度']
negatives = ['離婚率', '待機児童率', '1人1日当たりの排出量']
def zscore(s): return (s - s.mean()) / s.std()
df['happiness'] = (
    sum(zscore(df[c]) for c in positives) -
    sum(zscore(df[c]) for c in negatives)
) / (len(positives) + len(negatives))

# 説明変数はスコアの材料と重ならない列にする
feature_cols = ['総人口', '一般病院数', '保育所等数', '小学校教員数',
                '延べ宿泊者数', '年平均気温', '消費支出（二人以上の世帯）']
X, y = df[feature_cols].values, df['happiness'].values

kf = KFold(n_splits=5, shuffle=True, random_state=42)
r2s = []
for tr, va in kf.split(X):
    m = xgb.XGBRegressor(n_estimators=400, max_depth=3, learning_rate=0.03,
                         reg_lambda=1.0, random_state=42)
    m.fit(X[tr], y[tr])
    r2s.append(r2_score(y[va], m.predict(X[va])))
print('mean valid R^2 =', round(np.mean(r2s), 3))
print('std            =', round(np.std(r2s), 3))

# 全データで学習 → トップ 5 県とワースト 5 県
m_all = xgb.XGBRegressor(n_estimators=400, max_depth=3, learning_rate=0.03,
                         reg_lambda=1.0, random_state=42).fit(X, y)
df['pred_happiness'] = m_all.predict(X)
print('--- 予測幸福度トップ 5 ---')
print(df.nlargest(5, 'pred_happiness')[['都道府県','happiness','pred_happiness']].round(2).to_string(index=False))
print('--- 予測幸福度ワースト 5 ---')
print(df.nsmallest(5, 'pred_happiness')[['都道府県','happiness','pred_happiness']].round(2).to_string(index=False))