应用回归分析作业

第一步：导入数据

将数据的Features和Values分开，以便后续步骤可以分开使用（训练时只使用Features，验证时二者都要使用）。

import pandas as pd
import numpy as np
import visuals as vs

data = pd.read_csv('housing.csv')
features = data.drop('MEDV',axis=1)
prices = data['MEDV']

数据解释如下：
CRIM：城镇人均犯罪率。
ZN：住宅用地超过 25000 sq.ft. 的比例。
INDUS：城镇非零售商用土地的比例。
CHAS：查理斯河空变量（如果边界是河流，则为1；否则为0）。
NOX：一氧化氮浓度。
RM：住宅平均房间数。
AGE：1940 年之前建成的自用房屋比例。
DIS：到波士顿五个中心区域的加权距离。
RAD：辐射性公路的接近指数。
TAX：每 10000 美元的全值财产税率。
PTRATIO：城镇师生比例。
B：1000（Bk-0.63）^ 2，其中 Bk 指代城镇中黑人的比例。
LSTAT：人口中地位低下者的比例。
MEDV：自住房的平均房价，以千美元计。

data.info()

#查看数据集有没有缺失值
<class ‘pandas.core.frame.DataFrame’>
RangeIndex: 506 entries, 0 to 505
Data columns (total 14 columns):

Column Non-Null Count Dtype

---

0 CRIM 506 non-null float64
1 ZN 506 non-null float64
2 INDUS 506 non-null float64
3 CHAS 506 non-null int64
4 NOX 506 non-null float64
5 RM 506 non-null float64
6 AGE 506 non-null float64
7 DIS 506 non-null float64
8 RAD 506 non-null int64
9 TAX 506 non-null float64
10 PIRATIO 506 non-null float64
11 B 506 non-null float64
12 LSTAT 506 non-null float64
13 MEDV 506 non-null float64
dtypes: float64(12), int64(2)
memory usage: 55.5 KB

def find_outliers(model,X , y ,sigma=3):
    try :
        y_pred = pd.Series(model.predict(X), index=y.index)
    except:
        model.fit(X,y)
        y_pred = pd.Series(model.predict(X), index=y.index)
        
    resid = y - y_pred
    mean_resid = resid.mean()
    std_resid = resid.std()
    
    z = (resid - mean_resid) / std_resid
    outliers = z[abs(z)>sigma].index #相当于标准化残差（ZRE）大于3
    
    print('R2=',model.score(X,y))
    print('mse=',mean_squared_error(y,y_pred))
    
    print('mean of residuals:',mean_resid)
    print('std of residuals:',std_resid)
    print('-------------------------------------------')
    
    print('mean of residuals:',mean_resid)
    print('std of reiduals:',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.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')
    
    return outliers

from sklearn.linear_model import Ridge
from sklearn.metrics import mean_squared_error
X_train = data.iloc[:,0:-1]
y_train = data.iloc[:,-1]
outliers = find_outliers(Ridge(),X_train,y_train)

#通过岭回归查找异常值，红色点为异常值

train_cols = 6
train_rows = len(train_data.columns)
plt.figure(figsize=(4*train_cols,4*train_rows))

i = 0
for col in train_data.columns:
    i += 1
    ax=plt.subplot(train_rows,train_cols,i)
    sns.distplot(train_data[col],fit=stats.norm)
    
    i+=1
    ax=plt.subplot(train_rows,train_cols,i)
    res = stats.probplot(train_data[col],plot=plt)
    
plt.tight_layout()
plt.show()

通过画直方图和Q-Q图，如果数据满足正态分布，则所有点都会落在直线上。

from sklearn.model_selection import train_test_split

X_train, X_test, y_train, y_test = train_test_split(features, prices, test_size=0.2,random_state=1)

将数据划分为测试机和训练集

模型衡量标准
房价预测属于回归问题，要判断模型预测值与实际值的拟合程度，可以使用R2分数进行评价。R2是决定系数（coefficient of determination）,在统计学中，它表征回归方程多大程度上解释了因变量的变化，可用来表示表示对模型拟合程度的好坏，其取值范围是0~1，R2等于1表示模型可以100%解释目标值变化，R2等于0表示模型完全不能解释目标值的变化。

from sklearn.metrics import r2_score
def performance_metric(y_true,y_predict):
    '''计算实际值与预测值的R2分数'''
    score = r2_score(y_true,y_predict)
    return score

选择最优参数
网格搜索类似于穷举法，会将全部参数计算一遍，并根据评分方法（R2分数）计算得分，最后找出最优的参数。

from sklearn.model_selection import GridSearchCV
from sklearn.model_selection import KFold
from sklearn.tree import DecisionTreeRegressor
from sklearn.metrics import make_scorer

def fit_model(X,y):
    
    cross_validator = KFold(n_splits=10,shuffle=False)    #使用10折交叉验证
    
    regressor = DecisionTreeRegressor()
    
    params = {'max_depth':[1,2,3,4,5,6,7,8,9,10]}
    
    scoring_fnc = make_scorer(score_func = performance_metric)
    
    grid = GridSearchCV(estimator = regressor,
                        param_grid = params,
                        scoring = scoring_fnc,
                        cv = cross_validator)
    
    grid = grid.fit(X,y)    #基于数据X,y进行网格搜索
    
    return grid.best_estimator_    #返回网格搜索后的最优模型

optimal_reg = fit_model(X_train,y_train)    #使用数据对模型进行训练，找出最优参数

print('The best parameter is max_depth={}'.format(optimal_reg.get_params()['max_depth']))

The best parameter is max_depth=5

GridSearchCV参数

cross_validator参数表示使用的交叉验证方法，本例中定义为10折交叉验证，默认为3折交叉验证；
regressor参数表示使用的回归函数，本例使用决策树回归；
params参数表示要使用交叉验证的参数，本例中的参数为决策树的最大深度（max_depth）使用字典表示；
scoring_fnc表示评分函数，本例使用R2分数
从结果来看，使用网格搜索所得出的最优参数是max_depth=5。

我们已经训练好一个模型，现在我们想知道模型的预测能力怎样，使用在第二步分割好的X_test与y_test对模型的预测能力评分。

predicted_test = optimal_reg.predict(X_test)  #将X_test放入模型进行预测，并将预测结果放在变量predicted_test中

test_set_scoring = performance_metric(predicted_test,y_test)   #将预测值predicted_test与实际值y_test计算R2分数

print('The test set score is {}'.format(test_set_scoring))

The test set score is 0.8219648568947842
最终测试集R2分数为0.82。这个分数表示，房价82%的价格模型可以解释，还有18%的房价模型无法解释。
决策树模型还可以通过显示每一个变量对模型的重要性进行自变量的选择，在spss中我们是通过前进法（每次增加一个，知道没有可引进的变量为止，通过观测p值来看该变量是否显著，最后通过R方来判断选择的模型）和后退法（首先将所有变量都引进来再逐一剔除，最后通过R方进行选择）