In [1]:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.model_selection import train_test_split
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import OneHotEncoder
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import GridSearchCV
from sklearn.linear_model import LinearRegression
from sklearn.tree import DecisionTreeRegressor
from sklearn.ensemble import RandomForestRegressor
import warnings
warnings.filterwarnings('ignore')
In [2]:
housing = pd.read_csv('./datasets/housing/housing.csv')
In [3]:
housing.head()
Out[3]:
- 数据之间的量级相差很大。
- ocean_proximity 是 categorical value
In [4]:
housing.info()
- total_bedrooms 存在一些缺失值。
In [5]:
housing.describe()
Out[5]:
In [6]:
housing.hist(bins=50, figsize=(20,15))
plt.show()
- households, population, total_bedrooms, total_rooms 都有比较长的拖尾。可以考虑做log转换
- housing_median_age，median_house_value 最大值数据存在winsorization，也即把更高的数据统一设置成同一个值
Exploratory data analysis¶
- 所有数据项目，包括计量经济学项目，数据探索(exploratory data analysis, EDA)都是重要的一步。
- 在预测问题里，EDA 更重要，因为不必过分拘泥于前人的数据处理方法。在计量经济学项目中，如果前人的文献里包括了某个变量，那么很可能你也必须要包括这个变量。在预测问题里，这不是必然的。因此，充分探索数据，找到能使你的模型表现更好的变量或者变量组合更重要
In [7]:
numeric_features = ["housing_median_age","total_rooms","total_bedrooms",
"population","households","median_income","median_house_value"]
In [8]:
sns.pairplot(housing[numeric_features])
Out[8]:
- median_income 和 median_house_value 有正向关系。
- 因为数据的 winsorization，散点图中出现了很多点聚集在同一条水平线上的情况
- 考虑“家庭平均”,“人均”的变量组合
In [9]:
housing['rooms_per_household'] = housing['total_rooms'] / housing['households']
housing['bedrooms_per_household'] = housing['total_bedrooms'] / housing['households']
housing['rooms_per_population'] = housing['total_rooms'] / housing['population']
housing['bedrooms_per_population'] = housing['total_bedrooms'] / housing['population']
housing['population_per_household'] = housing['population'] / housing['households']
In [10]:
numeric_features = numeric_features + ['rooms_per_household','bedrooms_per_household',
'rooms_per_population','bedrooms_per_population',
'population_per_household']
In [11]:
numeric_features
Out[11]:
查看和median_house_value的相关性
In [12]:
housing.corr()['median_house_value'].sort_values(ascending=False)
Out[12]:
- median_income 和 median_house_value 相关性比较强
- 纬度也和 median_house_value 有较强相关性
- 上面只考虑了数值变量之间的联系。下面考察一下ocean_proximity，也即离海远近和房价之间的关系
In [13]:
housing.groupby('ocean_proximity').median_house_value.describe()
Out[13]:
很明显，靠海的房子更值钱
train test split¶
- 简单的策略：把数据8-2划分，80%的随机样本用作训练，剩下的20%来衡量模型好坏。
- 也可以做更复杂的划分，比如继续把训练样本划分为train, valiation。或者考虑数据的分布情况做抽样划分。
- 简单起见，这里做8-2划分
- 可以自己做随机抽样，把数据按比例划出来。也可以用
sklearn中的train_test_split
In [14]:
train, test = train_test_split(housing, test_size=0.2, random_state=42)
In [15]:
train_X = train.drop('median_house_value', axis=1)
train_y = train['median_house_value'].copy()
test_X = test.drop('median_house_value', axis=1)
test_y = test['median_house_value'].copy()
Data cleaning¶
根据以上的考察，有以下几个需要注意的地方：
- 缺失值处理
- categorical feature，也即 ocean_proximity 处理
- 量级统一
Missing values¶
对于缺失值，这里采取用median填充
In [16]:
imputer = SimpleImputer(strategy='median')
train_X_num = train_X.drop('ocean_proximity', axis=1)
train_X_num = imputer.fit_transform(train_X_num)
In [17]:
test_X_num = test_X.drop('ocean_proximity', axis=1)
test_X_num = imputer.transform(test_X_num)
Categorical feature¶
一般有两种处理分类变量的方法。
- 第一种是给每一类编号，比如'<1H OCEAN'编为0，'INLAND'编成1，以此类推。
- 第二种是虚拟变量。
- 第一类的缺点是编号之间实际上没有数量关系，但算法不会这样认为。例如，'<1H OCEAN'是0，'INLAND'是1，但这里的0和1并没有数量上的关系。
- 虚拟变量的缺点是分类太多时可能会有很多的变量被做出来。
这里用虚拟变量的方法处理ocean_proximity
In [18]:
train_X_cat = train_X[['ocean_proximity']]
one_hot = OneHotEncoder()
train_X_cat = one_hot.fit_transform(train_X_cat)
In [19]:
train_X_cat = train_X_cat.toarray()
In [20]:
test_X_cat = test_X[['ocean_proximity']]
test_X_cat = one_hot.transform(test_X_cat)
test_X_cat = test_X_cat.toarray()
Feature scaling¶
对于量级，这里采取Standard Scaler的处理办法，也即对每一个变量，减去其均值，除以其标准差
In [21]:
std_scaler = StandardScaler()
train_X_num = std_scaler.fit_transform(train_X_num)
In [22]:
test_X_num = std_scaler.transform(test_X_num)
Combine numerical and categorical columns¶
In [23]:
train_X_num.shape
Out[23]:
In [24]:
train_X_cat.shape
Out[24]:
In [25]:
train_X_prepared = np.hstack((train_X_num, train_X_cat))
train_X_prepared.shape
Out[25]:
In [26]:
test_X_prepared = np.hstack((test_X_num, test_X_cat))
In [27]:
test_X_num.shape
Out[27]:
In [28]:
test_X_cat.shape
Out[28]:
In [29]:
test_X_prepared.shape
Out[29]:
Fit models¶
Linear regression¶
In [30]:
lin_reg = LinearRegression()
lin_reg.fit(train_X_prepared, train_y)
Out[30]:
考察模型表现
In [31]:
lin_reg_predict = lin_reg.predict(train_X_prepared)
lin_reg_rmse = np.sqrt(mean_squared_error(train_y, lin_reg_predict))
In [32]:
lin_reg_rmse
Out[32]:
Decision tree¶
In [33]:
tree_reg = DecisionTreeRegressor()
tree_reg.fit(train_X_prepared, train_y)
Out[33]:
In [34]:
tree_predict = tree_reg.predict(train_X_prepared)
tree_rmse = np.sqrt(mean_squared_error(train_y, tree_predict))
In [35]:
tree_rmse
Out[35]:
很显然，决策树给出的模型产生了严重的过拟合。原因在于，sklearn默认的DecisionTreeRegressor对于很多参数不做任何限制。
如果我们对决策树的深度做出一定的限制，看看结果如何：
In [36]:
tree_reg = DecisionTreeRegressor(max_depth=5)
tree_reg.fit(train_X_prepared, train_y)
Out[36]:
In [37]:
tree_predict = tree_reg.predict(train_X_prepared)
tree_rmse = np.sqrt(mean_squared_error(train_y, tree_predict))
print(tree_rmse)
如何确定参数？可以把训练数据进一步分组进行 cross validation.
Decision tree with cross validation¶
我们考虑 $3*3 = 9$ 种参数组合。对于每一种参数组合，我们把训练数据分成5组，在其中4组种训练模型，剩下的1组上评估模型表现，并取5组表现的均值。然后选出表现最优的一组参数作为决策树的参数。
In [38]:
param_grid = [
{'max_depth': [5, 10, 20], 'max_features': [2,5,8]}
]
In [39]:
tree_reg = DecisionTreeRegressor()
grid_search = GridSearchCV(tree_reg, param_grid, cv=5,
scoring='neg_mean_squared_error',
return_train_score=True)
In [40]:
grid_search.fit(train_X_prepared, train_y)
Out[40]:
In [41]:
grid_search.best_params_
Out[41]:
由此可以确定，在我们选定的9组候选参数中，max_depth设为10，max_features设为8，决策树的表现是最好的。如果我们把备选参数增多，则可以选出更优的模型，但程序运行的时间也会相应变长。
In [42]:
cv_results = grid_search.cv_results_
for mean_score, params in zip(cv_results['mean_test_score'],
cv_results['params']):
print(np.sqrt(-mean_score), params)
In [43]:
tree_reg = DecisionTreeRegressor(max_depth=10,max_features=8)
tree_reg.fit(train_X_prepared, train_y)
Out[43]:
Ensemble methods¶
这里我们考察最简单的ensemble模型，Random Forest, 随机森林
In [44]:
forest_reg = RandomForestRegressor()
forest_reg.fit(train_X_prepared, train_y)
Out[44]:
In [45]:
forest_predict = forest_reg.predict(train_X_prepared)
forest_rmse = np.sqrt(mean_squared_error(train_y,forest_predict))
In [46]:
print(forest_rmse)
可以看到，随机森林给出的均方误差是最小的，同时，这个误差也没有像第一次的决策树那样不合理（为0）。因此，直观上，随机森林是最好的模型。
Evaluation on the test set¶
在测试集上最终看一下模型的表现优劣
In [47]:
test_lin_predict = lin_reg.predict(test_X_prepared)
test_tree_predict = tree_reg.predict(test_X_prepared)
test_forest_predict = forest_reg.predict(test_X_prepared)
test_lin_rmse = np.sqrt(mean_squared_error(test_y, test_lin_predict))
test_tree_rmse = np.sqrt(mean_squared_error(test_y, test_tree_predict))
test_forest_rmse = np.sqrt(mean_squared_error(test_y, test_forest_predict))
print('Linear Regression rmse on test set:', test_lin_rmse)
print('Decision Tree regression rmse on test set:', test_tree_rmse)
print('Random Forest regression rmse on test set:', test_forest_rmse)
可以看到，在测试集上，Random Forest 的表现是最好的。