import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import os
import random
import copy
from itertools import repeat
import shutil
import time
import math
import array

df = pd.read_csv('train.csv')

# drop columns with NaN
df.dropna(axis=1, inplace=True)
# drop irrelevant columns
df.drop(columns=['Id'], inplace=True)
# drop major outliers
df = df[df['BsmtFinSF1'] <= 5000]
df = df[df['MiscVal'] <= 5000]
df = df[df['LotArea'] <= 100000]

columns_non_categorical = df.select_dtypes(exclude='object').columns.to_list()
columns_non_categorical.sort()
columns_non_categorical.remove('SalePrice')
columns_non_categorical = ['SalePrice'] + columns_non_categorical


# alphabetically sort columns, keep target first
df_temp = df[['SalePrice']]
df.drop(columns=['SalePrice'], inplace=True)
df.sort_index(axis=1, inplace=True)
df = pd.concat([df_temp, df], axis=1)
del df_temp

X = df.drop(columns=['SalePrice']).copy()
y = df['SalePrice'].copy()
X = add_constant(X)
linear_model = sm.OLS(y, X).fit()
print(f'baseline BIC: {linear_model.bic}')

baseline BIC: 34570.166173470934

import statsmodels.api as sm
from mlxtend.feature_selection import SequentialFeatureSelector as SFS
from sklearn.base import BaseEstimator

class DummyEstimator(BaseEstimator):
    # mlxtend 希望使用 sklearn 估計器，但這里不需要（使用 statsmodels OLS 代替）。
    def fit(self, X, y=None, **kwargs):
        return self

def neg_bic(m, X, y):
    # objective function
    lin_mod_res = sm.OLS(y, X, hascnotallow=True).fit()
    return -lin_mod_res.bic

seq_selector = SFS(
    DummyEstimator(),
    k_features=(1, X.shape[1]),
    forward=True,
    floating=False,
    scoring=neg_bic,
    cv=None,
    n_jobs=-1,
    verbose=0,
    # make sure the intercept is not dropped
    fixed_features=['const'],
)

n_features = X.shape[1] - 1
objective_runs_sfs = round(n_features * (n_features + 1) / 2)
t_start_seq = time.time()
# mlxtend will mess with your dataframes if you don't .copy()
seq_res = seq_selector.fit(X.copy(), y.copy())
t_end_seq = time.time()
run_time_seq = t_end_seq - t_start_seq
seq_metrics = seq_res.get_metric_dict()

best k:         36
best objective: 33708.98602877906
R2 @ best k:    0.9075677543649224
objective runs: 22791
total run time: 42.448 sec

best_objective_seq = -np.inf
r2_of_best_k = 0
r2_list = []
best_k = 1
best_features_seq = []

# also extract R2 from the feature selection search
for k in tqdm(seq_metrics.keys()):
    r2_eval_mod = sm.OLS(
        y, X[list(seq_metrics[k]['feature_names'])], hascnotallow=True
    )
    r2_eval_mod_res = r2_eval_mod.fit()
    r2 = r2_eval_mod_res.rsquared
    r2_list.append(r2)
    score_k = seq_metrics[k]['avg_score']
    if score_k > best_objective_seq:
        best_objective_seq = score_k
        best_k = k
        best_features_seq = list(seq_metrics[k]['feature_names'])
        r2_of_best_k = r2

print(f'best k:         {best_k}')
print(f'best objective: {-best_objective_seq}')
print(f'R2 @ best k:    {r2_of_best_k}')
print(f'objective runs: {objective_runs_sfs}')
print(f'total run time: {run_time_seq:.3f} sec')
print(f'best features:  {best_features_seq}')

sfs_avg = [-seq_metrics[k]['avg_score'] for k in sorted(seq_metrics.keys())]
fig, ax = plt.subplots(1, 2, figsize=(10, 4), layout='constrained')
ax[0].scatter(sorted(seq_metrics.keys()), sfs_avg, s=1)
ax[0].set_xticks(sorted(seq_metrics.keys()), minor=True)
ax[0].set_title('BIC')
ax[0].set_xlabel('number of features')
ax[0].set_ylabel('BIC')
ax[1].scatter(sorted(seq_metrics.keys()), r2_list, s=1)
ax[1].set_title('R2')
ax[1].set_xlabel('number of features')
ax[1].set_ylabel('R2')
fig.suptitle('sequential forward search')
fig.savefig('sfs-performance.png')
fig.show()

best k:         36
best objective: 33708.98602877906
R2 @ best k:    0.9075677543649224
objective runs: 22791
total run time: 42.448 sec
best features:  ['const', 'BedroomAbvGr', 'BldgType_1Fam', 'BsmtFinSF1', 'Condition1_Norm', 
'Condition2_PosN', 'ExterQual_Ex', 'Exterior1st_BrkFace', 'Functional_Sev', 
'Functional_Typ', 'GarageArea', 'GrLivArea', 'HeatingQC_Ex', 'KitchenQual_Ex', 
'LandContour_HLS', 'LotArea', 'LotConfig_CulDSac', 'LowQualFinSF', 'MSSubClass', 
'Neighborhood_BrkSide', 'Neighborhood_Crawfor', 'Neighborhood_NWAmes', 
'Neighborhood_NoRidge', 'Neighborhood_NridgHt', 'Neighborhood_Somerst', 
'Neighborhood_StoneBr', 'OverallCond', 'OverallQual', 'PoolArea', 'RoofMatl_WdShngl', 
'SaleCondition_Normal', 'SaleType_New', 'ScreenPorch', 'TotalBsmtSF', 'WoodDeckSF', 
'YearBuilt']

'const', 'BedroomAbvGr', 'BldgType_1Fam', 'BsmtFinSF1', 
'Condition1_Norm', 'Condition2_PosN', 'ExterQual_Ex', 
'Exterior1st_BrkFace', 'Functional_Sev', 'Functional_Typ', 
'GarageArea', 'GrLivArea', 'HeatingQC_Ex', 'KitchenQual_Ex', 
'LandContour_HLS', 'LotArea', 'LotConfig_CulDSac', 'LowQualFinSF', 
'MSSubClass', 'Neighborhood_BrkSide', 'Neighborhood_Crawfor', 
'Neighborhood_NWAmes', 'Neighborhood_NoRidge', 'Neighborhood_NridgHt', 
'Neighborhood_Somerst', 'Neighborhood_StoneBr', 'OverallCond', 
'OverallQual', 'PoolArea', 'RoofMatl_WdShngl', 'SaleCondition_Normal', 
'SaleType_New', 'ScreenPorch', 'TotalBsmtSF', 'WoodDeckSF', 'YearBuilt'

def cma_objective(fs):
    features_use = ['const'] + [
        f for i, f in enumerate(features_select) if fs[i,] == 1
    ]
    lin_mod = sm.OLS(y_cmaes, X_cmaes[features_use], hascnotallow=True).fit()
    return lin_mod.bic


X_cmaes = X.copy()
y_cmaes = y.copy()
features_select = [f for f in X_cmaes.columns if f != 'const']

dim = len(features_select)
bounds = np.tile([0, 1], (dim, 1))
steps = np.ones(dim)
optimizer = CMAwM(
    mean=np.full(dim, 0.5),
    sigma=1 / 6,
    bounds=bounds,
    steps=steps,
    n_max_resampling=10 * dim,
    seed=0,
)

max_gen = 100
best_objective_cmaes_small = np.inf
best_sol_raw_cmaes_small = None
for gen in tqdm(range(max_gen)):
    solutions = []
    for _ in range(optimizer.population_size):
        x_for_eval, x_for_tell = optimizer.ask()
        value = cma_objective(x_for_eval)
        solutions.append((x_for_tell, value))
        if value < best_objective_cmaes_small:
            best_objective_cmaes_small = value
            best_sol_raw_cmaes_small = x_for_eval
    optimizer.tell(solutions)

best_features_cmaes_small = [
    features_select[i]
    for i, val in enumerate(best_sol_raw_cmaes_small.tolist())
    if val == 1.0
]
print(f'best objective: {best_objective_cmaes_small}')
print(f'best features:  {best_features_cmaes_small}')

def cma_objective(fs):
    features_use = ['const'] + [
        f for i, f in enumerate(features_select) if fs[i,] == 1
    ]
    lin_mod = sm.OLS(y_cmaes, X_cmaes[features_use], hascnotallow=True).fit()
    return lin_mod.bic
# copy the original dataframes into local copies, once
X_cmaes = X.copy()
y_cmaes = y.copy()

rs = np.random.RandomState(seed=0)
features_select = [f for f in X_cmaes.columns if f != 'const']
n_features = len(features_select)
cma_bounds = np.tile([0, 1], (n_features, 1))
cma_steps = np.ones(n_features)
n_max_resampling = 10 * n_features
mean = np.full(n_features, 0.5)
sigma = 1 / 6
pop_size = 4 + math.floor(3 * math.log(n_features))
margin = 1 / (n_features * pop_size)

optimizer = CMAwM(
    mean=mean,
    sigma=sigma,
    bounds=cma_bounds,
    steps=cma_steps,
    n_max_resampling=n_max_resampling,
    seed=0,
    population_size=pop_size,
    margin=margin,
)

gen_max = 1000
best_objective_cmaes = np.inf
best_generation_cmaes = 0
best_sol_raw_cmaes = None
history_values_cmaes = np.full((gen_max,), np.nan)
history_values_best_cmaes = np.full((gen_max,), np.nan)
time_to_best_cmaes = np.inf
objective_runs_cmaes = 0
solutions_avg = np.full((gen_max, n_features), np.nan)
n_jobs = os.cpu_count()

iterator = tqdm(range(gen_max), desc='generation')
t_start_cmaes = time.time()
for generation in iterator:
    best_value_gen = np.inf
    # solutions fed back to the optimizer
    solutions_float = []
    # binary-truncated solutions - the yes/no answers we're looking for
    solutions_binary = np.full((pop_size, n_features), np.nan)
    vals = np.full((pop_size,), np.nan)

    for i in range(pop_size):
        # re-seeding the RNG is very important
        # otherwise CMA-ES gets stuck in local extremes
        seed = rs.randint(1, 2**16) + generation
        optimizer._rng.seed(seed)
        fs_for_eval, fs_for_tell = optimizer.ask()
        solutions_binary[i,] = fs_for_eval
        value = cma_objective(fs_for_eval)
        objective_runs_cmaes += 1
        vals[i] = value
        solutions_float.append((fs_for_tell, value))

    optimizer.tell(solutions_float)
    solutions_avg[generation,] = solutions_binary.mean(axis=0)

    best_value_gen = vals.min()
    t_end_loop_cmaes = time.time()
    if best_value_gen < best_objective_cmaes:
        best_objective_cmaes = best_value_gen
        best_generation_cmaes = generation
        best_sol_raw_cmaes = solutions_binary[np.argmin(vals),]
        time_to_best_cmaes = t_end_loop_cmaes - t_start_cmaes
        print(
            f'gen: {best_generation_cmaes:5n}, new best objective: {best_objective_cmaes:.4f}'
        )
    history_values_cmaes[generation] = best_value_gen
    history_values_best_cmaes[generation] = best_objective_cmaes

    if optimizer.should_stop():
        print(f'Optimizer decided to stop.')
        iterator.close()
        break
    if os.path.isfile('break'):
        # to gracefully break the loop, manually create a file called 'break'
        print(f'Found break file, stopping now.')
        iterator.close()
        break
gen_completed_cmaes = generation

best_features_cmaes = [
    features_select[i]
    for i, val in enumerate(best_sol_raw_cmaes.tolist())
    if val == 1.0
]
print(f'best objective:  {best_objective_cmaes}')
print(f'best generation: {best_generation_cmaes}')
print(f'objective runs:  {objective_runs_cmaes}')
print(f'time to best:    {time_to_best_cmaes:.3f} sec')
print(f'best features:   {best_features_cmaes}')

cma_mv_df = pd.DataFrame(solutions_avg, columns=features_select).iloc[
    : gen_completed_cmaes + 1
]

if gen_completed_cmaes > 20:
    x_ticks = list(
        range(0, gen_completed_cmaes, round(gen_completed_cmaes / 20))
    )
else:
    x_ticks = list(range(0, gen_completed_cmaes))
if x_ticks[-1] != gen_completed_cmaes:
    x_ticks.append(gen_completed_cmaes)

fig, ax = plt.subplots(
    2,
    1,
    sharex=True,
    height_ratios=[20, 1],
    figsize=(12, 45),
    layout='constrained',
)
sns.heatmap(
    cma_mv_df.T,
    vmin=0.0,
    vmax=1.0,
    cmap='viridis',
    cbar=True,
    cbar_kws={'fraction': 0.01, 'anchor': (0.0, 1.0)},
    ax=ax[0],
)
ax[0].axvline(x=best_generation_cmaes, color='C7')
ax[0].tick_params(axis='both', reset=True)
ax[0].set_xticks(x_ticks)
ax[0].set_xticklabels(x_ticks)
ax[0].set_title('Population average of feature selector values')
ax[0].set_xlabel('generation')

ax[1].scatter(
    range(gen_completed_cmaes + 1),
    history_values_cmaes[: gen_completed_cmaes + 1],
    s=1,
    label='current value',
)
ax[1].plot(
    range(gen_completed_cmaes + 1),
    history_values_best_cmaes[: gen_completed_cmaes + 1],
    color='C1',
    label='best value',
)
ax[1].vlines(
    x=best_generation_cmaes,
    ymin=history_values_cmaes[: gen_completed_cmaes + 1].min(),
    ymax=history_values_cmaes[: gen_completed_cmaes + 1].max(),
    colors='C7',
)
ax[1].tick_params(axis='both', reset=True)
ax[1].legend()
ax[1].set_title('Objective value')
ax[1].set_xlabel('generation')
fig.suptitle('CMA-ES')
fig.savefig('cmaes-performance.png')
fig.show()

best objective:  33703.070530508514
best generation: 921
objective runs:  20000
time to best:    48.326 sec
best features:   ['BedroomAbvGr', 'BldgType_1Fam', 'BsmtFinSF1', 'Condition1_Norm', 
'Condition2_PosN', 'ExterQual_Ex', 'Exterior1st_BrkFace', 'Functional_Sev', 
'Functional_Typ', 'GarageCars', 'GrLivArea', 'HeatingQC_Ex', 'KitchenQual_Ex', 
'LandContour_HLS', 'LotArea', 'LowQualFinSF', 'Neighborhood_BrkSide', 
'Neighborhood_Crawfor', 'Neighborhood_NWAmes', 'Neighborhood_NoRidge', 
'Neighborhood_NridgHt', 'Neighborhood_Somerst', 'Neighborhood_StoneBr', 'OverallCond', 
'OverallQual', 'PoolArea', 'RoofMatl_WdShngl', 'SaleCondition_Normal', 'SaleType_New', 
'ScreenPorch', 'TotalBsmtSF', 'WoodDeckSF', 'YearBuilt']

上下滑動查看更多CMA-ES 優(yōu)化歷史

best objective:  33703.070530508514
best generation: 921
objective runs:  20000
time to best:    48.326 sec

baseline BIC: 34570.166173470934