工业蒸汽预测-07模型融合

1 基础步骤

1.1 导包

import warnings
warnings.filterwarnings("ignore")
import matplotlib.pyplot as plt
plt.rcParams.update({'figure.max_open_warning': 0})
import seaborn as sns

# modelling
import pandas as pd
import numpy as np
from scipy import stats
from sklearn.model_selection import train_test_split
from sklearn.model_selection import GridSearchCV, RepeatedKFold, cross_val_score,cross_val_predict,KFold
from sklearn.metrics import make_scorer,mean_squared_error
from sklearn.linear_model import LinearRegression, Lasso, Ridge, ElasticNet
from sklearn.svm import LinearSVR, SVR
from sklearn.neighbors import KNeighborsRegressor
from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor,AdaBoostRegressor
from xgboost import XGBRegressor
from sklearn.preprocessing import PolynomialFeatures,MinMaxScaler,StandardScaler

1.2 导入数据

#load_dataset
with open("./data/zhengqi_train.txt")  as fr:
    data_train=pd.read_table(fr,sep="\t")
with open("./data/zhengqi_test.txt") as fr_test:
    data_test=pd.read_table(fr_test,sep="\t")

1.3 合并数据

#merge train_set and test_set
data_train["oringin"]="train"
data_test["oringin"]="test"
data_all=pd.concat([data_train,data_test],axis=0,ignore_index=True)
data_all

	V0	V1	V2	V3	V4	V5	V6	V7	V8	V9	...	V30	V31	V32	V33	V34	V35	V36	V37	target	oringin
0	0.566	0.016	-0.143	0.407	0.452	-0.901	-1.812	-2.360	-0.436	-2.114	...	0.109	-0.615	0.327	-4.627	-4.789	-5.101	-2.608	-3.508	0.175	train
1	0.968	0.437	0.066	0.566	0.194	-0.893	-1.566	-2.360	0.332	-2.114	...	0.124	0.032	0.600	-0.843	0.160	0.364	-0.335	-0.730	0.676	train
2	1.013	0.568	0.235	0.370	0.112	-0.797	-1.367	-2.360	0.396	-2.114	...	0.361	0.277	-0.116	-0.843	0.160	0.364	0.765	-0.589	0.633	train
3	0.733	0.368	0.283	0.165	0.599	-0.679	-1.200	-2.086	0.403	-2.114	...	0.417	0.279	0.603	-0.843	-0.065	0.364	0.333	-0.112	0.206	train
4	0.684	0.638	0.260	0.209	0.337	-0.454	-1.073	-2.086	0.314	-2.114	...	1.078	0.328	0.418	-0.843	-0.215	0.364	-0.280	-0.028	0.384	train
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
4808	-1.362	-1.553	-3.096	-0.444	0.381	1.375	-4.854	-5.331	-4.074	-3.838	...	-4.488	-5.793	-4.050	-1.187	-0.852	-2.131	-2.564	0.597	NaN	test
4809	-2.698	-3.452	-3.620	-1.066	-1.385	1.378	-4.927	-5.103	-4.393	-1.683	...	-0.613	-7.698	-0.674	-1.187	-0.852	-2.131	-2.564	1.215	NaN	test
4810	-2.615	-3.564	-3.402	-0.422	-1.272	1.121	-4.223	-4.315	-5.196	-3.407	...	0.125	-6.111	0.275	-1.851	-1.548	-1.537	-2.544	1.612	NaN	test
4811	-2.661	-3.646	-3.271	-0.699	-1.270	1.116	-3.716	-3.809	-4.735	-2.976	...	1.086	-5.268	0.683	-1.645	-1.471	-1.537	-2.549	1.431	NaN	test
4812	-2.321	-3.037	-3.214	-1.594	-0.910	1.259	-3.616	-3.747	-4.368	-2.976	...	-0.774	-5.211	1.618	-1.703	-1.471	-1.537	-1.123	1.988	NaN	test

4813 rows × 40 columns

1.4 删除相关特征

data_all.drop(["V5","V9","V11","V17","V22","V28"],axis=1,inplace=True)

1.5 数据归一化

# normalise numeric columns
cols_numeric=list(data_all.columns)
cols_numeric.remove("oringin")
def scale_minmax(col):
    return (col-col.min())/(col.max()-col.min())
scale_cols = [col for col in cols_numeric if col!='target']
data_all[scale_cols] = data_all[scale_cols].apply(scale_minmax,axis=0)

1.6 画图：探查特征和标签相关信息

#Check effect of Box-Cox transforms on distributions of continuous variables

fcols = 6
frows = len(cols_numeric)-1
plt.figure(figsize=(4*fcols,4*frows))
i=0

for var in cols_numeric:
    if var!='target':
        dat = data_all[[var, 'target']].dropna() # 获取var列和target列。双重方括号[[  ]]用于选择多个列。
        
        i+=1
        plt.subplot(frows,fcols,i)
        sns.distplot(dat[var] , fit=stats.norm);
        plt.title(var+' Original')
        plt.xlabel('')
        
        i+=1
        plt.subplot(frows,fcols,i)
        _=stats.probplot(dat[var], plot=plt)
        plt.title('skew='+'{:.4f}'.format(stats.skew(dat[var])))
        plt.xlabel('')
        plt.ylabel('')
        
        i+=1
        plt.subplot(frows,fcols,i)
        plt.plot(dat[var], dat['target'],'.',alpha=0.5)
        plt.title('corr='+'{:.2f}'.format(np.corrcoef(dat[var], dat['target'])[0][1]))
 
        i+=1
        plt.subplot(frows,fcols,i)
        trans_var, lambda_var = stats.boxcox(dat[var].dropna()+1)
        trans_var = scale_minmax(trans_var)      
        sns.distplot(trans_var , fit=stats.norm);
        plt.title(var+' Tramsformed')
        plt.xlabel('')
        
        i+=1
        plt.subplot(frows,fcols,i)
        _=stats.probplot(trans_var, plot=plt)
        plt.title('skew='+'{:.4f}'.format(stats.skew(trans_var)))
        plt.xlabel('')
        plt.ylabel('')
        
        i+=1
        plt.subplot(frows,fcols,i)
        plt.plot(trans_var, dat['target'],'.',alpha=0.5)
        plt.title('corr='+'{:.2f}'.format(np.corrcoef(trans_var,dat['target'])[0][1]))

1.7 Box-Cox变换

对特征进行Box-Cox变换，使其满足正态性

cols_transform=data_all.columns[0:-2]
for col in cols_transform:   
    # transform column
    data_all.loc[:,col], _ = stats.boxcox(data_all.loc[:,col]+1)

代码解释

1. cols_transform = data_all.columns[0:-2]：获取data_all数据集中的所有列名，并使用切片操作[0:-2]选择从第一列到倒数第三列的所有列。

2. data_all.loc[:, col], _ = stats.boxcox(data_all.loc[:, col] + 1)：应用了Box-Cox变换来转换指定的列col。data_all.loc[:, col]用于选择DataFrame中的特定列。stats.boxcox()接受一个一维数组作为输入，并返回Box-Cox变换后的结果和变换参数。在这里，原始列的值会先加1（避免零值），然后应用Box-Cox变换。转换后的结果会更新到data_all中的相应列，而下划线_表示不需要使用的变换参数则被忽略。

1.8 分位数计算和绘图

标签数据统计转换后的数据，计算分位数画图展示（基于正态分布）

print(data_all.target.describe())

plt.figure(figsize=(12,4))
plt.subplot(1,2,1)
sns.distplot(data_all.target.dropna() , fit=stats.norm);
plt.subplot(1,2,2)
_=stats.probplot(data_all.target.dropna(), plot=plt)

count    2888.000000
mean        0.126353
std         0.983966
min        -3.044000
25%        -0.350250
50%         0.313000
75%         0.793250
max         2.538000
Name: target, dtype: float64

1.9 标签数据对数变换

标签数据对数变换数据，使数据更符合正态

#Log Transform SalePrice to improve normality
sp = data_train.target
data_train.target1 =np.power(1.5,sp)
print(data_train.target1.describe())

plt.figure(figsize=(12,4))
plt.subplot(1,2,1)
sns.distplot(data_train.target1.dropna(),fit=stats.norm);
plt.subplot(1,2,2)
_=stats.probplot(data_train.target1.dropna(), plot=plt)

count    2888.000000
mean        1.129957
std         0.394110
min         0.291057
25%         0.867609
50%         1.135315
75%         1.379382
max         2.798463
Name: target, dtype: float64

代码解释

1. sp = data_train.target：将data_train中的目标变量列赋值给变量sp。

2. data_train.target1 = np.power(1.5, sp)：将对数变换应用于目标变量sp，使用了np.power()函数来计算1.5的sp次方，也就是对目标变量进行了一个指数转换。转换后的结果存储在data_train数据集中一个名为target1的新列中。

3. print(data_train.target1.describe())：打印转换后的目标变量target1的描述统计信息，包括均值、标准差和分位数等。

4. plt.figure(figsize=(12,4))：创建一个图形窗口，设置大小为（12,4）。

5. plt.subplot(1,2,1)：在图形窗口中创建一个子图，总共有1行2列，当前子图位于第1个位置。

6. sns.distplot(data_train.target1.dropna(), fit=stats.norm)：绘制直方图和核密度估计图，其中data_train.target1.dropna()是转换后的目标变量数据，fit=stats.norm表示通过正态分布拟合曲线来比较数据的分布情况。

7. plt.subplot(1,2,2)：在图形窗口中创建一个子图，总共有1行2列，当前子图位于第2个位置。

8. _=stats.probplot(data_train.target1.dropna(), plot=plt)：使用stats.probplot()函数绘制Q-Q图，用于检验数据是否服从正态分布。data_train.target1.dropna()是转换后的目标变量数据，plot=plt表示将Q-Q图绘制在指定的子图上。

2 获取训练和测试数据

原先将训练集和测试集进行合并处理，现在进行拆分。使用简单交叉验证方法对模型进行验证，划为训练数据为70%，验证数据为30%

# function to get training samples
def get_training_data():
    # extract training samples
    from sklearn.model_selection import train_test_split
    df_train = data_all[data_all["oringin"]=="train"]
    df_train["label"]=data_train.target1
    # split SalePrice and features
    y = df_train.target
    X = df_train.drop(["oringin","target","label"],axis=1)
    X_train,X_valid,y_train,y_valid=train_test_split(X,y,test_size=0.3,random_state=100)
    return X_train,X_valid,y_train,y_valid

# extract test data (without SalePrice)
def get_test_data():
    df_test = data_all[data_all["oringin"]=="test"].reset_index(drop=True)
    return df_test.drop(["oringin","target"],axis=1)

3 模型评价函数

将RMSE和MSE作为模型性能的评价指标。

from sklearn.metrics import make_scorer
# metric for evaluation 自定义
def rmse(y_true, y_pred):
    diff = y_pred - y_true
    sum_sq = sum(diff**2)    
    n = len(y_pred)   
    return np.sqrt(sum_sq/n)

def mse(y_ture,y_pred):
    return mean_squared_error(y_ture,y_pred)

# scorer to be used in sklearn model fitting
rmse_scorer = make_scorer(rmse, greater_is_better=False) # greater_is_better=False 意味着rmse越小越好
mse_scorer = make_scorer(mse, greater_is_better=False)

4 处理异常数据

# function to detect outliers based on the predictions of a model
def find_outliers(model, X, y, sigma=3):

    # predict y values using model
    try:
        y_pred = pd.Series(model.predict(X), index=y.index)
    # if predicting fails, try fitting the model first
    except:
        model.fit(X,y)
        y_pred = pd.Series(model.predict(X), index=y.index)
        
    # calculate residuals between the model prediction and true y values
    resid = y - y_pred
    mean_resid = resid.mean()
    std_resid = resid.std()

    # calculate z statistic, define outliers to be where |z|>sigma
    z = (resid - mean_resid)/std_resid    
    outliers = z[abs(z)>sigma].index
    
    # print and plot the results
    print('R2=',model.score(X,y))
    print('rmse=',rmse(y, y_pred))
    print("mse=",mean_squared_error(y,y_pred))
    print('---------------------------------------')

    print('mean of residuals:',mean_resid)
    print('std of residuals:',std_resid)
    print('---------------------------------------')

    print(len(outliers),'outliers:')
    print(outliers.tolist())

    plt.figure(figsize=(15,5))
    ax_131 = plt.subplot(1,3,1)
    plt.plot(y,y_pred,'.')
    plt.plot(y.loc[outliers],y_pred.loc[outliers],'ro')
    plt.legend(['Accepted','Outlier'])
    plt.xlabel('y')
    plt.ylabel('y_pred');

    ax_132=plt.subplot(1,3,2)
    plt.plot(y,y-y_pred,'.')
    plt.plot(y.loc[outliers],y.loc[outliers]-y_pred.loc[outliers],'ro')
    plt.legend(['Accepted','Outlier'])
    plt.xlabel('y')
    plt.ylabel('y - y_pred');

    ax_133=plt.subplot(1,3,3)
    z.plot.hist(bins=50,ax=ax_133)
    z.loc[outliers].plot.hist(color='r',bins=50,ax=ax_133)
    plt.legend(['Accepted','Outlier'])
    plt.xlabel('z')
    
    plt.savefig('outliers.png')
    
    return outliers

# get training data
from sklearn.linear_model import Ridge
X_train, X_valid,y_train,y_valid = get_training_data()
test=get_test_data()

# find and remove outliers using a Ridge model
outliers = find_outliers(Ridge(), X_train, y_train)

# permanently remove these outliers from the data
#df_train = data_all[data_all["oringin"]=="train"]
#df_train["label"]=data_train.target1
#df_train=df_train.drop(outliers)
X_outliers=X_train.loc[outliers]
y_outliers=y_train.loc[outliers]
X_t=X_train.drop(outliers)
y_t=y_train.drop(outliers)

R2= 0.8794138468263233
rmse= 0.34510338853546346
mse= 0.11909634877865903
---------------------------------------
mean of residuals: -2.96425702171071e-16
std of residuals: 0.3451887995969211
---------------------------------------
23 outliers:
[2655, 2159, 1164, 2863, 1145, 2697, 2528, 2645, 691, 1085, 1874, 2647, 776, 2625, 884, 2696, 2668, 1310, 1901, 2769, 2002, 2669, 1040]

5 网格搜索训练模型

def get_trainning_data_omitoutliers():
    y1=y_t.copy()
    X1=X_t.copy()
    return X1,y1

使用网格搜索训练模型

from sklearn.preprocessing import StandardScaler
def train_model(model, param_grid=[], X=[], y=[], 
                splits=5, repeats=5):

    # get unmodified training data, unless data to use already specified 
    if len(y)==0:
        X,y = get_trainning_data_omitoutliers()
        #poly_trans=PolynomialFeatures(degree=2)
        #X=poly_trans.fit_transform(X)
        #X=MinMaxScaler().fit_transform(X)
    
    # create cross-validation method
    rkfold = RepeatedKFold(n_splits=splits, n_repeats=repeats)
    
    # perform a grid search if param_grid given
    if len(param_grid)>0:
        # setup grid search parameters
        gsearch = GridSearchCV(model, param_grid, cv=rkfold,
                               scoring="neg_mean_squared_error",
                               verbose=1, return_train_score=True)

        # search the grid
        gsearch.fit(X,y)

        # extract best model from the grid
        model = gsearch.best_estimator_        
        best_idx = gsearch.best_index_ # 获取最佳模型的索引
 
        # get cv-scores for best model
        grid_results = pd.DataFrame(gsearch.cv_results_) # .cv_results_ 每个参数组合的得分、标准差等信息       
        cv_mean = abs(grid_results.loc[best_idx,'mean_test_score'])
        cv_std = grid_results.loc[best_idx,'std_test_score']

    # no grid search, just cross-val score for given model    
    else:
        grid_results = []
        cv_results = cross_val_score(model, X, y, scoring="neg_mean_squared_error", cv=rkfold)
        cv_mean = abs(np.mean(cv_results))
        cv_std = np.std(cv_results)
    
    # combine mean and std cv-score in to a pandas series
    cv_score = pd.Series({'mean':cv_mean,'std':cv_std})

    # predict y using the fitted model
    y_pred = model.predict(X)
    
    # print stats on model performance         
    print('----------------------')
    print(model)
    print('----------------------')
    print('score=',model.score(X,y))
    print('rmse=',rmse(y, y_pred))
    print('mse=',mse(y, y_pred))
    print('cross_val: mean=',cv_mean,', std=',cv_std)
    
    # residual plots
    y_pred = pd.Series(y_pred,index=y.index)
    resid = y - y_pred
    mean_resid = resid.mean()
    std_resid = resid.std()
    z = (resid - mean_resid)/std_resid    
    n_outliers = sum(abs(z)>3)
    
    plt.figure(figsize=(15,5))
    ax_131 = plt.subplot(1,3,1)
    plt.plot(y,y_pred,'.')
    plt.xlabel('y')
    plt.ylabel('y_pred');
    plt.title('corr = {:.3f}'.format(np.corrcoef(y,y_pred)[0][1]))
    ax_132=plt.subplot(1,3,2)
    plt.plot(y,y-y_pred,'.')
    plt.xlabel('y')
    plt.ylabel('y - y_pred');
    plt.title('std resid = {:.3f}'.format(std_resid))
    
    ax_133=plt.subplot(1,3,3)
    z.plot.hist(bins=50,ax=ax_133)
    plt.xlabel('z')
    plt.title('{:.0f} samples with z>3'.format(n_outliers))

    return model, cv_score, grid_results

# places to store optimal models and scores
opt_models = dict()
score_models = pd.DataFrame(columns=['mean','std'])

# no. k-fold splits
splits=5
# no. k-fold iterations
repeats=5

6 单一模型训练

6.1 岭回归

使用岭回归模型对数据进行预测，采用RMSE,MSE等指标对模型性能进行评价

model = 'Ridge'

opt_models[model] = Ridge()
alph_range = np.arange(0.25,6,0.25)
param_grid = {'alpha': alph_range}

opt_models[model],cv_score,grid_results = train_model(opt_models[model], param_grid=param_grid, 
                                              splits=splits, repeats=repeats)

cv_score.name = model
score_models = score_models.append(cv_score)

plt.figure()
plt.errorbar(alph_range, abs(grid_results['mean_test_score']),
             abs(grid_results['std_test_score'])/np.sqrt(splits*repeats))
plt.xlabel('alpha')
plt.ylabel('score')

Fitting 25 folds for each of 23 candidates, totalling 575 fits
----------------------
Ridge(alpha=0.25)
----------------------
score= 0.8963563380306225
rmse= 0.31899649233050603
mse= 0.1017587621191668
cross_val: mean= 0.10644934256965052 , std= 0.005302313654798051





Text(0, 0.5, 'score')

​

​

上面的图形主要反映了模型预测的准确度及拟合情况，后面的模型也将进行类似的可视化分析。这里对图形展示的信息进行介绍，从左至右依次如下：

(1)真实值（横轴：y)与模型预测值（竖轴：y pred)的散点图，图形上方显示了相关性数值，其越接近1越好。对于岭回归模型，相关性数值为0.947，预测值与真实值比较一致。

(2)其为在交叉验证训练模型时，真实值（横轴：y)与模型预测值和真实值的残差（竖轴：y-y_pred)的散点图，图形上方显示了方差，其越小说明模型越稳定。可以看到，对于岭回归模型，在真实值y=-3附近的预测值有较大的偏差，同时，方差为0.319，较为稳定。

(3)图是由模型预测值和真实值的残差（横轴：z=(resid-mean resid)/std resid)与落在按z轴划分区间的频率（竖轴：频数）所画的直方图，图形上方显示了预测值与真实值的残差大于三倍标准差的数，其越小越好，越大说明预测中有些样本的偏差很大。对于岭回归模型，预测值与真实值的残差大于三倍标准差的数为5个，模型对偏差大的数据有较好的包容性。

(4)岭回归模型的参数（横轴：alpha)与模型的评价指标MSE(竖轴：score)的误差棒图。
可以看出，对于岭回归模型，随着alpha的增大，评价指标MSE的数值也越来越大，其模型的方差也逐渐增大。

​

代码解释

1.

cv_score.name = model
score_models = score_models.append(cv_score)
````

```python
print(cv_score.name)
print("-------")
print(score_models)

# 结果
Ridge
-------
           mean       std
Ridge  0.106211  0.007446

---

2.

plt.errorbar(alph_range, abs(grid_results['mean_test_score']),
             abs(grid_results['std_test_score'])/np.sqrt(splits*repeats))

代码根据给定的数据点集合在图表上绘制条形图，并添加以标准误差为高度的误差线。

• errorbar() 用于绘制带有误差线的条形图。

参数解释:

• alph_range 是 x 轴上的数据点集合，它表示自变量的取值范围。
• abs(grid_results['mean_test_score']) 是 y 轴上的数据点集合，它表示不同自变量取值下的平均测试分数的绝对值。
• abs(grid_results['std_test_score'])/np.sqrt(splits*repeats) 是误差条的高度。grid_results['std_test_score'] 表示测试分数的标准差，而 splits 和 repeats 是进行交叉验证时划分数据集的折数和重复次数，通过除以sqrt(splits*repeats) 来将标准差转换为标准误差。

​

6.2 Lasso回归

使用Lasso回归模型对数据进行预测，采用RMSE,MSE等指标对模型性能进行评价

model = 'Lasso'

opt_models[model] = Lasso()
alph_range = np.arange(1e-4,1e-3,4e-5)
param_grid = {'alpha': alph_range}

opt_models[model], cv_score, grid_results = train_model(opt_models[model], param_grid=param_grid, 
                                              splits=splits, repeats=repeats)

cv_score.name = model
score_models = score_models.append(cv_score)

plt.figure()
plt.errorbar(alph_range, abs(grid_results['mean_test_score']),abs(grid_results['std_test_score'])/np.sqrt(splits*repeats))
plt.xlabel('alpha')
plt.ylabel('score')

Fitting 25 folds for each of 23 candidates, totalling 575 fits
----------------------
Lasso(alpha=0.0001)
----------------------
score= 0.8965376022911505
rmse= 0.3187174209135977
mse= 0.10158079439381539
cross_val: mean= 0.1062923603647315 , std= 0.007051871077111589



Text(0, 0.5, 'score')

​

​

6.3 ElasticNet 回归

model ='ElasticNet'
opt_models[model] = ElasticNet()

param_grid = {'alpha': np.arange(1e-4,1e-3,1e-4),
              'l1_ratio': np.arange(0.1,1.0,0.1),
              'max_iter':[100000]}

opt_models[model], cv_score, grid_results = train_model(opt_models[model], param_grid=param_grid, 
                                              splits=splits, repeats=1)

cv_score.name = model
score_models = score_models.append(cv_score)

Fitting 5 folds for each of 81 candidates, totalling 405 fits
----------------------
ElasticNet(alpha=0.0001, l1_ratio=0.9, max_iter=100000)
----------------------
score= 0.8965268226984906
rmse= 0.3187340238271355
mse= 0.10159137794503677
cross_val: mean= 0.10558612445735098 , std= 0.006931804612280223

6.4 SVR回归

model='LinearSVR'
opt_models[model] = LinearSVR()

crange = np.arange(0.1,1.0,0.1)
param_grid = {'C':crange,
             'max_iter':[1000]}

opt_models[model], cv_score, grid_results = train_model(opt_models[model], param_grid=param_grid, 
                                              splits=splits, repeats=repeats)


cv_score.name = model
score_models = score_models.append(cv_score)

plt.figure()
plt.errorbar(crange, abs(grid_results['mean_test_score']),abs(grid_results['std_test_score'])/np.sqrt(splits*repeats))
plt.xlabel('C')
plt.ylabel('score')

Fitting 25 folds for each of 9 candidates, totalling 225 fits
----------------------
LinearSVR(C=0.4)
----------------------
score= 0.30645151350483657
rmse= 0.8251880824335185
mse= 0.680935371390307
cross_val: mean= 0.9136545268818084 , std= 0.5904526473477509





Text(0, 0.5, 'score')

​

​

6.5 K近邻

model = 'KNeighbors'
opt_models[model] = KNeighborsRegressor()

param_grid = {'n_neighbors':np.arange(3,11,1)}

opt_models[model], cv_score, grid_results = train_model(opt_models[model], param_grid=param_grid, 
                                              splits=splits, repeats=1)

cv_score.name = model
score_models = score_models.append(cv_score)

plt.figure()
plt.errorbar(np.arange(3,11,1), abs(grid_results['mean_test_score']),abs(grid_results['std_test_score'])/np.sqrt(splits*1))
plt.xlabel('n_neighbors')
plt.ylabel('score')

Fitting 5 folds for each of 8 candidates, totalling 40 fits
----------------------
KNeighborsRegressor(n_neighbors=10)
----------------------
score= 0.7188007679386812
rmse= 0.5254381453364182
mse= 0.27608524457457456
cross_val: mean= 0.35006829225711783 , std= 0.04461643597140863





Text(0, 0.5, 'score')

​

​

7 模型融合 Boosting方法

把n个比较弱的模型，组合成一个比较强的模型。类似于随机森林。

7.1 GBDT 模型

使用 Gradient Boosting模型对数据进行预测，采用RMSE,MSE等指标对模型性能进行评价

model = 'GradientBoosting'
opt_models[model] = GradientBoostingRegressor()

param_grid = {'n_estimators':[150,250,350],
              'max_depth':[1,2,3],
              'min_samples_split':[5,6,7]}

opt_models[model], cv_score, grid_results = train_model(opt_models[model], param_grid=param_grid, 
                                              splits=splits, repeats=1)

cv_score.name = model
score_models = score_models.append(cv_score)

Fitting 5 folds for each of 27 candidates, totalling 135 fits
----------------------
GradientBoostingRegressor(min_samples_split=6, n_estimators=250)
----------------------
score= 0.9677208586746097
rmse= 0.17802275580033733
mse= 0.03169210158274656
cross_val: mean= 0.09849105444854306 , std= 0.014499755846113458

7.2 XGBoost模型

model = 'XGB'
opt_models[model] = XGBRegressor()

param_grid = {'n_estimators':[100,200,300,400,500],
              'max_depth':[1,2,3],
             }

opt_models[model], cv_score,grid_results = train_model(opt_models[model], param_grid=param_grid, 
                                              splits=splits, repeats=1)

cv_score.name = model
score_models = score_models.append(cv_score)

Fitting 5 folds for each of 15 candidates, totalling 75 fits
----------------------
XGBRegressor(base_score=0.5, booster='gbtree', callbacks=None,
             colsample_bylevel=1, colsample_bynode=1, colsample_bytree=1,
             early_stopping_rounds=None, enable_categorical=False,
             eval_metric=None, gamma=0, gpu_id=-1, grow_policy='depthwise',
             importance_type=None, interaction_constraints='',
             learning_rate=0.300000012, max_bin=256, max_cat_to_onehot=4,
             max_delta_step=0, max_depth=3, max_leaves=0, min_child_weight=1,
             missing=nan, monotone_constraints='()', n_estimators=100, n_jobs=0,
             num_parallel_tree=1, predictor='auto', random_state=0, reg_alpha=0,
             reg_lambda=1, ...)
----------------------
score= 0.9702586237542838
rmse= 0.1708815062712164
mse= 0.029200489185519693
cross_val: mean= 0.10472580078722642 , std= 0.006953064644267524

7.3 随机森林模型

model = 'RandomForest'
opt_models[model] = RandomForestRegressor()

param_grid = {'n_estimators':[100,150,200],
              'max_features':[8,12,16,20,24],
              'min_samples_split':[2,4,6]}

opt_models[model], cv_score, grid_results = train_model(opt_models[model], param_grid=param_grid, 
                                              splits=5, repeats=1)

cv_score.name = model
score_models = score_models.append(cv_score)

Fitting 5 folds for each of 45 candidates, totalling 225 fits
----------------------
RandomForestRegressor(max_features=16, n_estimators=150)
----------------------
score= 0.9856162598758819
rmse= 0.11883666269492393
mse= 0.01412215240046713
cross_val: mean= 0.10239417326207083 , std= 0.006602546408163406

8 多模型预测 Bagging方法

Bagging 通常对分类任务使用简单投票法，对回归任务使用简单平均法.

def model_predict(test_data,test_y=[],stack=False):
    #poly_trans=PolynomialFeatures(degree=2)
    #test_data1=poly_trans.fit_transform(test_data)
    #test_data=MinMaxScaler().fit_transform(test_data)
    i=0
    y_predict_total=np.zeros((test_data.shape[0],))
    for model in opt_models.keys():
        if model!="LinearSVR" and model!="KNeighbors":
            y_predict=opt_models[model].predict(test_data)
            y_predict_total+=y_predict
            i+=1
        if len(test_y)>0:
            print("{}_mse:".format(model),mean_squared_error(y_predict,test_y))
    y_predict_mean=np.round(y_predict_total/i,3) # 模型融合的mean
    if len(test_y)>0:
        print("mean_mse:",mean_squared_error(y_predict_mean,test_y))
    else:
        y_predict_mean=pd.Series(y_predict_mean)
        return y_predict_mean

model_predict(X_valid,y_valid)

Ridge_mse: 0.137671879096647
Lasso_mse: 0.13786142881850508
ElasticNet_mse: 0.1378268557408619
LinearSVR_mse: 0.1378268557408619
KNeighbors_mse: 0.1378268557408619
GradientBoosting_mse: 0.1345410925334943
XGB_mse: 0.14010339249433112
RandomForest_mse: 0.13880263283624247
mean_mse: 0.12569497923875433

代码解释

def model_predict(test_data, test_y=[], stack=False):
    i = 0
    y_predict_total = np.zeros((test_data.shape[0],))

这里定义了一个名为 model_predict 的函数，它接受三个参数：test_data（待预测的测试数据），test_y（可选的真实标签），stack（是否进行堆叠）。函数中初始化了变量 i 和 y_predict_total，其中 i 用于计算参与预测的模型数量，y_predict_total 是一个形状与 test_data 行数相同的全零数组，用于累加每个模型的预测结果。

for model in opt_models.keys():
    if model != "LinearSVR" and model != "KNeighbors":
        y_predict = opt_models[model].predict(test_data)
        y_predict_total += y_predict
        i += 1
    if len(test_y) > 0:
        print("{}_mse:".format(model), mean_squared_error(y_predict, test_y))

这部分使用 for 循环遍历 opt_models 字典的键。这个字典存储了优化后的模型。在循环中，根据特定的条件判断，跳过了特定的两个模型（”LinearSVR” 和 “KNeighbors”）。然后，使用对应模型 model 对测试数据 test_data 进行预测，并将预测结果累加到 y_predict_total 中。同时，更新 i 的值。如果提供了真实标签 test_y，则计算并输出该模型的均方误差（MSE）。

y_predict_mean = np.round(y_predict_total / i, 3)

这里计算了平均预测结果 y_predict_mean，通过将 y_predict_total 除以 i 得到。

if len(test_y) > 0:
    print("mean_mse:", mean_squared_error(y_predict_mean, test_y))
else:
    y_predict_mean = pd.Series(y_predict_mean)
    return y_predict_mean

这部分首先判断是否提供了真实标签 test_y。如果提供了，计算并输出平均预测结果 y_predict_mean 和真实标签之间的均方误差（MSE）。如果没有提供真实标签，则将 y_predict_mean 转换为 Pandas Series 对象并返回。

代码的整体思路是：

• 遍历存储了优化后模型的字典 opt_models 的键。
• 对于每个模型，根据特定条件进行预测，然后将预测结果累加到 y_predict_total 中。
• 计算平均预测结果并返回或输出与真实标签的均方误差（如果提供了真实标签）。

目的是对一组模型进行预测，并计算它们的平均预测结果或与真实标签的均方误差。如果没有提供真实标签，它将返回一个包含平均预测结果的 Pandas Series 对象。

可以看到，模型融合预测的MSE最小，预测性能最优。

9 多模型融合 Stacking方法

https://blog.csdn.net/m0_47256162/article/details/119979540

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec # 创建多维子图网格的模块
import itertools
from sklearn.linear_model import LogisticRegression
from sklearn.svm import SVC
from sklearn.ensemble import RandomForestClassifier

##主要使用pip install mlxtend安装mlxtend
from mlxtend.classifier import EnsembleVoteClassifier
from mlxtend.data import iris_data
from mlxtend.plotting import plot_decision_regions
%matplotlib inline

# Initializing Classifiers
clf1 = LogisticRegression(random_state=0)
clf2 = RandomForestClassifier(random_state=0)
clf3 = SVC(random_state=0, probability=True)
eclf = EnsembleVoteClassifier(clfs=[clf1, clf2, clf3], weights=[2, 1, 1], voting='soft')

# Loading some example data
X, y = iris_data()
X = X[:,[0, 2]]

# Plotting Decision Regions
gs = gridspec.GridSpec(2, 2)
fig = plt.figure(figsize=(10, 8))

for clf, lab, grd in zip([clf1, clf2, clf3, eclf],
                         ['Logistic Regression', 'Random Forest', 'RBF kernel SVM', 'Ensemble'],
                         itertools.product([0, 1], repeat=2)):
    clf.fit(X, y)
    ax = plt.subplot(gs[grd[0], grd[1]])
    fig = plot_decision_regions(X=X, y=y, clf=clf, legend=2)
    plt.title(lab)
plt.show()

​

​

代码解释

EnsembleVoteClassifier()

用于创建集成投票分类器。它可以将多个基分类器的预测结果进行投票，从而得到集成分类器的最终预测结果。

参数包括：

• clfs：一个列表，包含要集成的基分类器对象。
• voting：设置投票方式，可选的取值为'hard'或'soft'，分别表示硬投票和软投票。硬投票是指根据少数服从多数的原则选择最多数量的类别作为预测结果，而软投票则是根据分类器预测的概率进行加权计算，并选择概率加权总和最大的类别作为预测结果。
• weights：一个可选的列表，用于指定每个基分类器的权重。如果未指定，则默认所有分类器的权重相等。

---

1. for clf, lab, grd in zip([clf1, clf2, clf3, eclf], ['Logistic Regression', 'Random Forest', 'RBF kernel SVM', 'Ensemble'], itertools.product([0, 1], repeat=2)):：

   • zip([clf1, clf2, clf3, eclf], ['Logistic Regression', 'Random Forest', 'RBF kernel SVM', 'Ensemble'], itertools.product([0, 1], repeat=2))将分类器列表、标签列表和网格坐标的笛卡尔积打包在一起，返回一个迭代器。每次迭代返回一个元组(clf, lab, grd)，其中clf表示分类器对象，lab表示标签，grd表示网格坐标。
   • for clf, lab, grd in ...:遍历迭代器中的每个元组，依次赋值给变量clf、lab和grd。

2. clf.fit(X, y)：使用当前分类器clf拟合数据集X和目标变量y，即进行训练。

3. ax = plt.subplot(gs[grd[0], grd[1]])：在网格规则对象gs上创建一个子图对象，并将其赋值给变量ax。通过索引grd[0]和grd[1]，选择网格规则中的坐标位置。

4. fig = plot_decision_regions(X=X, y=y, clf=clf, legend=2)：使用mlxtend.plotting模块的plot_decision_regions函数，在当前子图上绘制决策边界。传入参数包括输入数据集X、目标变量y、当前分类器clf，以及将图例放置的位置legend=2。绘制结果会返回一个图形对象，将其赋值给变量fig。

5. plt.title(lab)：设置当前子图的标题为相应的标签lab。

6. 最后一行plt.show()：显示整个图形，包含所有子图和决策边界。

通过这段代码，可以实现对每个分类器的训练和决策边界的绘制，并在同一个图形中进行比较和展示。每个分类器都会有一个独立的子图，显示相应的决策边界和标签。最后，通过plt.show()将整个图形显示出来。

10工业蒸汽赛题多模型融合 stacking方法

10.1 基础代码

from sklearn.model_selection import train_test_split
import pandas as pd
import numpy as np
from scipy import sparse
import xgboost
import lightgbm

from sklearn.ensemble import RandomForestRegressor,AdaBoostRegressor,GradientBoostingRegressor,ExtraTreesRegressor
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error


# 针对于一个模型进行K折交叉验证的函数
def stacking_reg(clf,train_x,train_y,test_x,clf_name,kf,label_split=None):
    train=np.zeros((train_x.shape[0],1))
    test=np.zeros((test_x.shape[0],1))
    test_pre=np.empty((folds,test_x.shape[0],1)) # 对话详解
    cv_scores=[]
    for i,(train_index,test_index) in enumerate(kf.split(train_x,label_split)):       
        tr_x=train_x[train_index]
        tr_y=train_y[train_index]
        te_x=train_x[test_index]
        te_y = train_y[test_index]
        if clf_name in ["rf","ada","gb","et","lr","lsvc","knn"]:
            clf.fit(tr_x,tr_y)
            pre=clf.predict(te_x).reshape(-1,1) # 对话详解
            train[test_index]=pre
            test_pre[i,:]=clf.predict(test_x).reshape(-1,1)
            cv_scores.append(mean_squared_error(te_y, pre))
        elif clf_name in ["xgb"]:
            train_matrix = clf.DMatrix(tr_x, label=tr_y, missing=-1) # 对话详解
            test_matrix = clf.DMatrix(te_x, label=te_y, missing=-1)
            z = clf.DMatrix(test_x, label=te_y, missing=-1)
            params = {'booster': 'gbtree',
                      'eval_metric': 'rmse',
                      'gamma': 1,
                      'min_child_weight': 1.5,
                      'max_depth': 5,
                      'lambda': 10,
                      'subsample': 0.7,
                      'colsample_bytree': 0.7,
                      'colsample_bylevel': 0.7,
                      'eta': 0.03,
                      'tree_method': 'exact',
                      'seed': 2017,
                      'nthread': 12
                      }
            num_round = 10000
            early_stopping_rounds = 100
            watchlist = [(train_matrix, 'train'),
                         (test_matrix, 'eval')
                         ]
            if test_matrix:
                model = clf.train(params, train_matrix, num_boost_round=num_round,evals=watchlist,
                                  early_stopping_rounds=early_stopping_rounds
                                  )
                pre= model.predict(test_matrix,ntree_limit=model.best_ntree_limit).reshape(-1,1)
                train[test_index]=pre
                test_pre[i, :]= model.predict(z, ntree_limit=model.best_ntree_limit).reshape(-1,1)
                cv_scores.append(mean_squared_error(te_y, pre))

        elif clf_name in ["lgb"]:
            train_matrix = clf.Dataset(tr_x, label=tr_y)
            test_matrix = clf.Dataset(te_x, label=te_y)
            #z = clf.Dataset(test_x, label=te_y)
            #z=test_x
            params = {
                      'boosting_type': 'gbdt',
                      'objective': 'regression_l2',
                      'metric': 'mse',
                      'min_child_weight': 1.5,
                      'num_leaves': 2**5,
                      'lambda_l2': 10,
                      'subsample': 0.7,
                      'colsample_bytree': 0.7,
                      'colsample_bylevel': 0.7,
                      'learning_rate': 0.03,
                      'tree_method': 'exact',
                      'seed': 2017,
                      'nthread': 12,
                      'silent': True,
                      }
            num_round = 10000
            callbacks=[lightgbm.log_evaluation(period=100), lightgbm.early_stopping(stopping_rounds=100)]
            if test_matrix:
                model = clf.train(params, train_matrix,num_round,valid_sets=test_matrix,
                                  callbacks = callbacks
                                  )
                pre= model.predict(te_x,num_iteration=model.best_iteration).reshape(-1,1)
                train[test_index]=pre
                test_pre[i, :]= model.predict(test_x, num_iteration=model.best_iteration).reshape(-1,1)
                cv_scores.append(mean_squared_error(te_y, pre))
        else:
            raise IOError("Please add new clf.")
        print("%s now score is:"%clf_name,cv_scores)
    test[:]=test_pre.mean(axis=0)
    print("%s_score_list:"%clf_name,cv_scores)
    print("%s_score_mean:"%clf_name,np.mean(cv_scores))
    return train.reshape(-1,1),test.reshape(-1,1)

代码详解

这段代码是一个用于堆叠（stacking）回归模型的函数。下面是对每一行代码的逐行解释：

1. def stacking_reg(clf,train_x,train_y,test_x,clf_name,kf,label_split=None):
   这个函数定义了一个堆叠回归模型的函数，它接受以下参数：

   • clf: 使用的基础回归模型
   • train_x: 训练集的特征数据
   • train_y: 训练集的目标数据
   • test_x: 测试集的特征数据
   • clf_name: 基础回归模型的名称
   • kf: 交叉验证的迭代器
   • label_split（可选）: 标签的分割方式

2. train=np.zeros((train_x.shape[0],1))
   创建一个形状为(训练集样本数, 1)的全零数组，并将其赋值给变量train。这个数组用于存储训练集在基础模型上的预测结果。

3. test=np.zeros((test_x.shape[0],1))
   创建一个形状为(测试集样本数, 1)的全零数组，并将其赋值给变量test。这个数组用于存储测试集在基础模型上的预测结果。

4. test_pre=np.empty((folds,test_x.shape[0],1))
   创建一个形状为(folds, 测试集样本数, 1)的空数组，并将其赋值给变量test_pre。这个数组用于存储每次交叉验证中测试集在基础模型上的预测结果。

5. cv_scores=[]
   创建一个空列表，并将其赋值给变量cv_scores。这个列表用于存储每次交叉验证的均方误差（MSE）。

6. for i,(train_index,test_index) in enumerate(kf.split(train_x,label_split)):
   遍历交叉验证迭代器kf生成的训练集和验证集的索引。每次迭代中，会得到一个训练集的索引和一个验证集的索引。

7. tr_x=train_x[train_index]
   根据索引从训练集中获取相应的特征数据，并赋值给变量tr_x。

8. tr_y=train_y[train_index]
   根据索引从训练集中获取相应的目标数据，并赋值给变量tr_y。

9. te_x=train_x[test_index]
   根据索引从训练集中获取相应的特征数据，并赋值给变量te_x。

10. te_y=train_y[test_index]
    根据索引从训练集中获取相应的目标数据，并赋值给变量te_y。

15-62. 基于不同的基础回归模型进行训练和预测，并计算均方误差（MSE）：

• 如果clf_name在["rf","ada","gb","et","lr","lsvc","knn"]中，使用clf.fit()方法进行训练，然后使用clf.predict()方法进行预测，并将预测结果存储在相应的变量中。
• 如果clf_name是"xgb"，则使用XGBoost库进行训练和预测。具体的参数设置和训练过程与XGBoost相关。
• 如果clf_name是"lgb"，则使用LightGBM库进行训练和预测。具体的参数设置和训练过程与LightGBM相关。

64. raise IOError("Please add new clf.")
    如果clf_name不属于任何已定义的模型名称，则抛出一个错误。

65. print("%s now score is:"%clf_name,cv_scores)
    打印当前基础回归模型的均方误差。

66. test[:]=test_pre.mean(axis=0)
    计算所有交叉验证模型的预测结果的平均值，并将结果赋值给test变量。

70-71. 分别打印基础回归模型的均方误差列表和均值。

73. return train.reshape(-1,1),test.reshape(-1,1)
    返回训练集和测试集的预测结果，其中预测结果的形状为(样本数, 1)。

---

15-62.代码是一个用于训练和评估XGBoost或LightGBM模型的逻辑。我将逐行解释代码的功能：

1. train_matrix = clf.DMatrix(tr_x, label=tr_y, missing=-1): 使用XGBoost库的DMatrix方法创建训练数据矩阵，其中tr_x是特征数据，tr_y是目标数据，missing=-1表示缺失值的表示方式。

2. test_matrix = clf.DMatrix(te_x, label=te_y, missing=-1): 使用XGBoost库的DMatrix方法创建测试数据矩阵，其中te_x是特征数据，te_y是目标数据，missing=-1表示缺失值的表示方式。

3. z = clf.DMatrix(test_x, label=te_y, missing=-1): 使用XGBoost库的DMatrix方法创建用于预测的数据矩阵，其中test_x是特征数据，te_y是目标数据，missing=-1表示缺失值的表示方式。

4. params = {...}: 设置XGBoost模型的参数，包括树的深度、学习率、正则化等。

5. num_round = 10000: 设置迭代轮数。

6. early_stopping_rounds = 100: 设置在验证集上早停的轮数。

7. watchlist = [(train_matrix, 'train'), (test_matrix, 'eval')]: 创建监控列表用于跟踪训练和测试集上的性能。

8. model = clf.train(params, train_matrix, num_boost_round=num_round,evals=watchlist,early_stopping_rounds=early_stopping_rounds): 使用XGBoost库的train方法训练模型，其中params是模型参数，train_matrix是训练数据矩阵，num_boost_round是迭代轮数，evals是监控列表，early_stopping_rounds是早停的轮数。

9. pre = model.predict(test_matrix,ntree_limit=model.best_ntree_limit).reshape(-1,1): 对测试数据进行预测，并将结果重塑为一列。

10. train[test_index] = pre: 将对测试集的预测结果保存在训练集的相应位置。

11. test_pre[i, :] = model.predict(z, ntree_limit=model.best_ntree_limit).reshape(-1,1): 对预测数据进行预测，并将结果保存在test_pre数组的第i行。

12. cv_scores.append(mean_squared_error(te_y, pre)): 将平均均方误差（MSE）添加到交叉验证分数列表中。

13. train_matrix = clf.Dataset(tr_x, label=tr_y): 使用LightGBM库的Dataset方法创建训练数据集。

14. test_matrix = clf.Dataset(te_x, label=te_y): 使用LightGBM库的Dataset方法创建测试数据集。

15. params = {...}: 设置LightGBM模型的参数，包括树的深度、学习率、正则化等。

16. model = clf.train(params, train_matrix,num_round,valid_sets=test_matrix,early_stopping_rounds=early_stopping_rounds): 使用LightGBM库的train方法训练模型，其中params是模型参数，train_matrix是训练数据集，num_round是迭代轮数，valid_sets是用于验证的数据集，early_stopping_rounds是早停的轮数。

17. pre = model.predict(te_x,num_iteration=model.best_iteration).reshape(-1,1): 对测试数据进行预测，并将结果重塑为一列。

18. train[test_index] = pre: 将对测试集的预测结果保存在训练集的相应位置。

19. test_pre[i, :] = model.predict(test_x, num_iteration=model.best_iteration).reshape(-1,1): 对预测数据进行预测，并将结果保存在test_pre数组的第i行。

20. cv_scores.append(mean_squared_error(te_y, pre)): 将平均均方误差（MSE）添加到交叉验证分数列表中。

代码整体思路

1. 创建一些空数组和变量，包括train（训练集预测结果）、test（测试集预测结果）、test_pre（测试集在各折交叉验证下的预测结果）和cv_scores（交叉验证得分列表）。

2. 对于交叉验证的每一折，进行以下操作：

   • 根据当前折的索引，获取训练集的索引和测试集的索引。
   • 根据索引从训练集和标签中抽取对应的样本和标签。
   • 如果分类器名称在[‘rf’、’ada’、’gb’、’et’、’lr’、’lsvc’、’knn’]中，则使用分类器训练模型，进行预测，并将预测结果存储在相应的数组中。
   • 如果分类器名称为’xgb’，则使用XGBoost训练模型，进行预测，并将预测结果存储在相应的数组中。
   • 如果分类器名称为’lgb’，则使用LightGBM训练模型，进行预测，并将预测结果存储在相应的数组中。
   • 计算当前折的均方误差得分，并将其添加到cv_scores列表中。

3. 计算所有折的平均预测结果，并将其存储在test数组中。

4. 打印分类器名称、交叉验证得分列表和平均得分。

5. 将train和test数组重新整形为列向量，并返回它们作为函数的输出。

10.2 模型融合stacking基学习器

def rf_reg(x_train, y_train, x_valid, kf, label_split=None):
    randomforest = RandomForestRegressor(n_estimators=600, max_depth=20, n_jobs=-1, random_state=2017, max_features="auto",verbose=1)
    rf_train, rf_test = stacking_reg(randomforest, x_train, y_train, x_valid, "rf", kf, label_split=label_split)
    return rf_train, rf_test,"rf_reg"

def ada_reg(x_train, y_train, x_valid, kf, label_split=None):
    adaboost = AdaBoostRegressor(n_estimators=30, random_state=2017, learning_rate=0.01)
    ada_train, ada_test = stacking_reg(adaboost, x_train, y_train, x_valid, "ada", kf, label_split=label_split)
    return ada_train, ada_test,"ada_reg"

def gb_reg(x_train, y_train, x_valid, kf, label_split=None):
    gbdt = GradientBoostingRegressor(learning_rate=0.04, n_estimators=100, subsample=0.8, random_state=2017,max_depth=5,verbose=1)
    gbdt_train, gbdt_test = stacking_reg(gbdt, x_train, y_train, x_valid, "gb", kf, label_split=label_split)
    return gbdt_train, gbdt_test,"gb_reg"

def et_reg(x_train, y_train, x_valid, kf, label_split=None):
    extratree = ExtraTreesRegressor(n_estimators=600, max_depth=35, max_features="auto", n_jobs=-1, random_state=2017,verbose=1)
    et_train, et_test = stacking_reg(extratree, x_train, y_train, x_valid, "et", kf, label_split=label_split)
    return et_train, et_test,"et_reg"

def lr_reg(x_train, y_train, x_valid, kf, label_split=None):
    lr_reg=LinearRegression(n_jobs=-1)
    lr_train, lr_test = stacking_reg(lr_reg, x_train, y_train, x_valid, "lr", kf, label_split=label_split)
    return lr_train, lr_test, "lr_reg"

def xgb_reg(x_train, y_train, x_valid, kf, label_split=None):
    xgb_train, xgb_test = stacking_reg(xgboost, x_train, y_train, x_valid, "xgb", kf, label_split=label_split)
    return xgb_train, xgb_test,"xgb_reg"

def lgb_reg(x_train, y_train, x_valid, kf, label_split=None):
    lgb_train, lgb_test = stacking_reg(lightgbm, x_train, y_train, x_valid, "lgb", kf, label_split=label_split)
    return lgb_train, lgb_test,"lgb_reg"

代码解释

1. rf_reg 函数：

   • 首先，使用随机森林回归器 RandomForestRegressor 创建一个随机森林模型。设置了一些参数，如树的数量(n_estimators)、最大深度(max_depth)、并行运算的CPU核数(n_jobs)等。
   • 接下来，调用 stacking_reg 函数，将创建的随机森林模型作为基学习器进行堆叠集成学习。该函数的作用是对训练数据 x_train 进行 K 折交叉验证，并返回基学习器在训练集和验证集上的预测结果。
   • 最后，返回训练集的预测结果 rf_train、验证集的预测结果 rf_test 和标识符字符串 "rf_reg"。

2. ada_reg 函数：

   • 首先，使用 AdaBoost 回归器 AdaBoostRegressor 创建一个 AdaBoost 模型。设置了一些参数，如基学习器的数量(n_estimators)、学习率(learning_rate)等。
   • 接下来，调用 stacking_reg 函数，将创建的 AdaBoost 模型作为基学习器进行堆叠集成学习。该函数的作用和上述相同，对训练数据 x_train 进行 K 折交叉验证，并返回基学习器在训练集和验证集上的预测结果。
   • 最后，返回训练集的预测结果 ada_train、验证集的预测结果 ada_test 和标识符字符串 "ada_reg"。

······

代码的目的是创建不同的基学习器，并将其作为基础模型进行堆叠集成学习。通过堆叠集成学习，可以利用各个基学习器的优势，提高整体模型的预测性能。最后，返回训练集和验证集上的预测结果以及标识符字符串。

10.3 定义模型融合stacking预测函数

def stacking_pred(x_train, y_train, x_valid, kf, clf_list, label_split=None, clf_fin="lgb", if_concat_origin=True):
    # 第一层-结果为train test
    for k, clf_list in enumerate(clf_list):
        clf_list = [clf_list]
        column_list = []
        train_data_list=[]
        test_data_list=[]
        for clf in clf_list:
            train_data,test_data,clf_name=clf(x_train, y_train, x_valid, kf, label_split=label_split)
            train_data_list.append(train_data)
            test_data_list.append(test_data)
            column_list.append("clf_%s" % (clf_name))
    train = np.concatenate(train_data_list, axis=1)
    test = np.concatenate(test_data_list, axis=1)
     
    # 是否选择原始特征拼接
    if if_concat_origin:
        train = np.concatenate([x_train, train], axis=1)
        test = np.concatenate([x_valid, test], axis=1)
    print(x_train.shape)
    print(train.shape)
    print(clf_name)
    print(clf_name in ["lgb"])
    # 第二层
    if clf_fin in ["rf","ada","gb","et","lr","lsvc","knn"]:
        if clf_fin in ["rf"]:
            clf = RandomForestRegressor(n_estimators=600, max_depth=20, n_jobs=-1, random_state=2017, max_features="auto",verbose=1)
        elif clf_fin in ["ada"]:
            clf = AdaBoostRegressor(n_estimators=30, random_state=2017, learning_rate=0.01)
        elif clf_fin in ["gb"]:
            clf = GradientBoostingRegressor(learning_rate=0.04, n_estimators=100, subsample=0.8, random_state=2017,max_depth=5,verbose=1)
        elif clf_fin in ["et"]:
            clf = ExtraTreesRegressor(n_estimators=600, max_depth=35, max_features="auto", n_jobs=-1, random_state=2017,verbose=1)
        elif clf_fin in ["lr"]:
            clf = LinearRegression(n_jobs=-1)
        clf.fit(train, y_train)
        pre = clf.predict(test).reshape(-1,1)
        return pred
    elif clf_fin in ["xgb"]:
        clf = xgboost
        train_matrix = clf.DMatrix(train, label=y_train, missing=-1)
        test_matrix = clf.DMatrix(train, label=y_train, missing=-1)
        params = {'booster': 'gbtree',
                  'eval_metric': 'rmse',
                  'gamma': 1,
                  'min_child_weight': 1.5,
                  'max_depth': 5,
                  'lambda': 10,
                  'subsample': 0.7,
                  'colsample_bytree': 0.7,
                  'colsample_bylevel': 0.7,
                  'eta': 0.03,
                  'tree_method': 'exact',
                  'seed': 2017,
                  'nthread': 12
                  }
        num_round = 10000
        #early_stopping_rounds = 100
        callbacks=[lightgbm.log_evaluation(period=100), lightgbm.early_stopping(stopping_rounds=100)]
        watchlist = [(train_matrix, 'train'),
                     (test_matrix, 'eval')
                     ]
        model = clf.train(params, train_matrix, num_boost_round=num_round,evals=watchlist,
                          #early_stopping_rounds=early_stopping_rounds
                          callbacks=callbacks
                          )
        pre = model.predict(test,ntree_limit=model.best_ntree_limit).reshape(-1,1)
        return pre
    elif clf_fin in ["lgb"]:
        print(clf_name)
        clf = lightgbm
        train_matrix = clf.Dataset(train, label=y_train)
        test_matrix = clf.Dataset(train, label=y_train)
        params = {
                  'boosting_type': 'gbdt',
                  'objective': 'regression_l2',
                  'metric': 'mse',
                  'min_child_weight': 1.5,
                  'num_leaves': 2**5,
                  'lambda_l2': 10,
                  'subsample': 0.7,
                  'colsample_bytree': 0.7,
                  'colsample_bylevel': 0.7,
                  'learning_rate': 0.03,
                  'tree_method': 'exact',
                  'seed': 2017,
                  'nthread': 12,
                  'silent': True,
                  }
        num_round = 10000
        #early_stopping_rounds = 100
        callbacks=[lightgbm.log_evaluation(period=100), lightgbm.early_stopping(stopping_rounds=100)]
        model = clf.train(params, train_matrix,num_round,valid_sets=test_matrix,
                          #early_stopping_rounds=early_stopping_rounds
                          callbacks=callbacks
                          )
        print('pred')
        pre = model.predict(test,num_iteration=model.best_iteration).reshape(-1,1)
        print(pre)
        return pre

10.4 使用lr_reg和lgb_reg进行融合预测

加载数据

# load_dataset
with open("./data/zhengqi_train.txt")  as fr:
    data_train=pd.read_table(fr,sep="\t")
with open("./data/zhengqi_test.txt") as fr_test:
    data_test=pd.read_table(fr_test,sep="\t")
    

K折交叉验证

# K折交叉验证
from sklearn.model_selection import StratifiedKFold, KFold

folds = 5
seed = 1
kf = KFold(n_splits=5, shuffle=True, random_state=0)

训练集和测试集数据

# 训练集
x_train = data_train[data_test.columns].values
y_train = data_train['target'].values
# 测试集
x_valid = data_test[data_test.columns].values

使用lr_reg和lgb_reg进行融合预测

clf_list = [lr_reg, lgb_reg]

##很容易过拟合
pred = stacking_pred(x_train, y_train, x_valid, kf, clf_list, label_split=None, clf_fin="lgb", if_concat_origin=True)

预测结果

pred

array([[ 0.41700189],
       [ 0.37101204],
       [ 0.1479981 ],
       ...,
       [-2.48570445],
       [-2.4264629 ],
       [-2.48367665]])

# 将pred转换为一维数组
pred_flat = pred.flatten()

# 保存为文本文件
np.savetxt('pred.txt', pred_flat)