import pandas as pd, numpy as np
from sklearn.linear_model import LinearRegression, Ridge, Lasso, ElasticNet
from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import KFold, cross_val_score, RepeatedKFold

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='utf-8', skiprows=1)
X = df[['一人当たり県民所得','世帯人員','高齢化率','人口密度','就業率']]
y = df['持ち家比率']

models = {
    'OLS': LinearRegression(),
    'Ridge(α=1)': Ridge(alpha=1.0),
    'Ridge(α=10)': Ridge(alpha=10.0),
    'Lasso(α=0.1)': Lasso(alpha=0.1, max_iter=10000),
    'ElasticNet': ElasticNet(alpha=0.1, l1_ratio=0.5, max_iter=10000),
    'RF': RandomForestRegressor(n_estimators=300, random_state=42),
    'GBM': GradientBoostingRegressor(n_estimators=200, random_state=42),
}

cv = RepeatedKFold(n_splits=5, n_repeats=10, random_state=42)
print(f'{"model":15} {"RMSE_mean":>10} {"RMSE_std":>10}')
for name, m in models.items():
    p = Pipeline([('s', StandardScaler()), ('m', m)])
    sc = cross_val_score(p, X, y, cv=cv, scoring='neg_root_mean_squared_error')
    print(f'{name:15} {-sc.mean():10.3f} {sc.std():10.3f}')

from scipy import stats
from sklearn.model_selection import cross_val_score

cv = RepeatedKFold(n_splits=5, n_repeats=20, random_state=42)
s_ridge = -cross_val_score(Pipeline([('s', StandardScaler()), ('m', Ridge(alpha=10))]),
                            X, y, cv=cv, scoring='neg_root_mean_squared_error')
s_ols   = -cross_val_score(Pipeline([('s', StandardScaler()), ('m', LinearRegression())]),
                            X, y, cv=cv, scoring='neg_root_mean_squared_error')

# 同じ fold での対応比較
t, p = stats.ttest_rel(s_ols, s_ridge)
print(f't={t:.3f}, p={p:.4f}')
print(f'Ridge は OLS より RMSE が小さい？: {(s_ridge < s_ols).mean():.0%} の fold で勝利')

from sklearn.model_selection import RepeatedKFold, cross_validate

cv = RepeatedKFold(n_splits=5, n_repeats=20, random_state=42)
results = cross_validate(
    Pipeline([('s', StandardScaler()), ('m', Ridge(alpha=10))]),
    X, y, cv=cv,
    scoring=['neg_root_mean_squared_error', 'r2', 'neg_mean_absolute_error'],
    return_train_score=True, n_jobs=-1
)
import numpy as np
print('RMSE:', -np.mean(results['test_neg_root_mean_squared_error']))
print('R²  :',  np.mean(results['test_r2']))
print('MAE :', -np.mean(results['test_neg_mean_absolute_error']))
print('train-test gap (R²):', np.mean(results['train_r2']) - np.mean(results['test_r2']))

import matplotlib.pyplot as plt
from sklearn.model_selection import cross_val_score

cv = RepeatedKFold(n_splits=5, n_repeats=30, random_state=42)
scores = {}
for name, m in [('OLS', LinearRegression()), ('Ridge', Ridge(alpha=10)),
                ('Lasso', Lasso(alpha=0.1, max_iter=10000)), ('RF', RandomForestRegressor(300, random_state=42))]:
    p = Pipeline([('s', StandardScaler()), ('m', m)])
    scores[name] = -cross_val_score(p, X, y, cv=cv, scoring='neg_root_mean_squared_error')

fig, ax = plt.subplots(figsize=(8, 4.5))
ax.boxplot(scores.values(), labels=scores.keys())
ax.set_ylabel('CV RMSE'); ax.set_title('Repeated 5-fold × 30 — 持ち家比率予測')
plt.tight_layout(); plt.savefig('cv_compare.png', dpi=110)

import statsmodels.api as sm

candidates = {
    'M1: 県民所得のみ':       ['一人当たり県民所得'],
    'M2: + 高齢化率':         ['一人当たり県民所得', '高齢化率'],
    'M3: + 人口密度':         ['一人当たり県民所得', '高齢化率', '人口密度'],
    'M4: フル':               ['一人当たり県民所得', '世帯人員', '高齢化率', '人口密度', '就業率'],
}
print(f'{"model":25} {"AIC":>8} {"BIC":>8} {"adj R²":>8}')
for name, cols in candidates.items():
    Xc = sm.add_constant(df[cols])
    res = sm.OLS(df['持ち家比率'], Xc).fit()
    print(f'{name:25} {res.aic:8.1f} {res.bic:8.1f} {res.rsquared_adj:8.3f}')

import optuna
from sklearn.ensemble import RandomForestRegressor

def objective(trial):
    params = {
        'n_estimators': trial.suggest_int('n_estimators', 100, 800),
        'max_depth':    trial.suggest_int('max_depth', 2, 20),
        'min_samples_leaf': trial.suggest_int('min_samples_leaf', 1, 10),
        'max_features': trial.suggest_float('max_features', 0.3, 1.0),
    }
    p = Pipeline([('s', StandardScaler()),
                   ('m', RandomForestRegressor(random_state=42, **params))])
    sc = cross_val_score(p, X, y, cv=5, scoring='neg_root_mean_squared_error')
    return -sc.mean()

study = optuna.create_study(direction='minimize', sampler=optuna.samplers.TPESampler(seed=42))
study.optimize(objective, n_trials=50, show_progress_bar=False)
print('Best RMSE:', study.best_value, ' params:', study.best_params)

from sklearn.experimental import enable_halving_search_cv  # noqa
from sklearn.model_selection import HalvingGridSearchCV
import numpy as np

p = Pipeline([('s', StandardScaler()), ('m', Ridge())])
search = HalvingGridSearchCV(
    p, param_grid={'m__alpha': np.logspace(-3, 3, 20)},
    cv=5, scoring='neg_root_mean_squared_error',
    factor=3, random_state=42
)
search.fit(X, y)
print('Best α:', search.best_params_, 'CV RMSE:', -search.best_score_)

from sklearn.model_selection import GridSearchCV
import numpy as np

alphas = np.logspace(-3, 3, 30)
g = GridSearchCV(Pipeline([('s', StandardScaler()), ('m', Ridge())]),
                  {'m__alpha': alphas}, cv=10,
                  scoring='neg_root_mean_squared_error', return_train_score=False)
g.fit(X, y)

# 最良 α より RMSE が「最良 + 1SE」以内で、 α が最大（最も正則化が強い＝最も単純）を選ぶ
means = -g.cv_results_['mean_test_score']
sds   =  g.cv_results_['std_test_score']
best_idx = means.argmin()
threshold = means[best_idx] + sds[best_idx]
eligible = np.where(means <= threshold)[0]
one_se_idx = eligible[np.argmax(alphas[eligible])]
print(f'best α = {alphas[best_idx]:.3f}, RMSE = {means[best_idx]:.3f}')
print(f'1-SE α = {alphas[one_se_idx]:.3f}, RMSE = {means[one_se_idx]:.3f}（より頑健）')

from sklearn.ensemble import HistGradientBoostingRegressor
from sklearn.model_selection import train_test_split

X_tr, X_va, y_tr, y_va = train_test_split(X, y, test_size=0.25, random_state=42)
m = HistGradientBoostingRegressor(
    max_iter=1000, learning_rate=0.05,
    early_stopping=True, validation_fraction=0.2,
    n_iter_no_change=20, random_state=42
)
m.fit(X_tr, y_tr)
print('停止 iter:', m.n_iter_)
print('val RMSE :', ((m.predict(X_va)-y_va)**2).mean()**0.5)