回归分析--python应用篇（scikit-learn）

最新推荐文章于 2024-07-29 18:31:04 发布

一瑞111

最新推荐文章于 2024-07-29 18:31:04 发布

阅读量2.7k

点赞数 1

分类专栏： algorithm

本文链接：https://blog.csdn.net/weixin_41712808/article/details/85323417

版权

algorithm 专栏收录该内容

6 篇文章 2 订阅

订阅专栏

实例一：目前有汽车数据，找到相关规律并做汽车价格预测。

1.调取相关工具包、读取数据、观察数据

2.清理数据：数据格式化、缺失值处理、异常值处理、预处理、特征相关性等

4.预测及检验

实例一：目前有汽车数据，找到相关规律并做汽车价格预测。

数据集简介
主要包括3类指标:	类别属性	连续指标
汽车的各种特性. 保险风险评级：(-3, -2, -1, 0, 1, 2, 3). 每辆保险车辆年平均相对损失支付.	make: 汽车的商标（奥迪，宝马。。。） fuel-type: 汽油还是天然气 aspiration: 涡轮 num-of-doors: 两门还是四门 body-style: 硬顶车、轿车、掀背车、敞篷车 drive-wheels: 驱动轮 engine-location: 发动机位置 engine-type: 发动机类型 num-of-cylinders: 几个气缸 fuel-system: 燃油系统	bore: continuous from 2.54 to 3.94. stroke: continuous from 2.07 to 4.17. compression-ratio: continuous from 7 to 23. horsepower: continuous from 48 to 288. peak-rpm: continuous from 4150 to 6600. city-mpg: continuous from 13 to 49. highway-mpg: continuous from 16 to 54. price: continuous from 5118 to 45400.

1.调取相关工具包、读取数据、观察数据

import numpy as np
import pandas as pd
from pandas import datetime
#pandas文档：http://pandas.pydata.org/pandas-docs/stable/generated/pandas.DataFrame.html#pandas.DataFrame

#作图库
import matplotlib.pyplot as plt
#matplotlib文档：https://matplotlib.org/gallery/index.html
import seaborn as sns 
import missingno as msno # missing values 缺失值可视化展示
#missingno文档 https://github.com/ResidentMario/missingno

#只适用jupyter notebook,在console直接生成图像，不适用spyder等编辑器
%matplotlib inline 

# machine learning 机器学习工具包
from statsmodels.distributions.empirical_distribution import ECDF  #ECDF：累积分布函数
from sklearn.metrics import mean_squared_error, r2_score
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import Lasso, LassoCV 
from sklearn.model_selection import train_test_split, cross_val_score 
from sklearn.ensemble import RandomForestRegressor
seed = 123 #做交叉验证，随机策略，制定随机种子，把随机因素控制住

f=open(r'C:\Users\Administrator\Desktop\算法案例\数学基础\回归分析\Auto-Data.csv')
#csv数据中的‘？’填充为缺失值
data = pd.read_csv(f,na_values = '?')
#观察数据
print(data.dtypes,data.shape)

#describe():对pandas进行数据统计，展示各列的均值、最大小值等
data.describe()

2.清理数据：数据格式化、缺失值处理、异常值处理、预处理、特征相关性等

2.1 缺失值处理

方法：1、缺失值较少是，1%以下，可以直接去掉nan；2、用已有的值取平均值或众值；3、用已有的数做回归模型，再用其他特征数据预测缺失值。

msno.matrix(data)

msno库可用作观察缺失值，matrix作图是黑白风格，易于观察，白色代表是缺失值，观上图，可知'normalized-losses'列缺失较为严重，不能单纯删除处理。我们先将缺失这列数据提出来看。

#pd.isnull返回的是bool类型，再传入pd，返回需要的值
print(data[pd.isnull(data['normalized-losses'])].head())

#指定seaborn画图风格ticks
sns.set(style = "ticks") #https://blog.csdn.net/llh_1178/article/details/77923033
plt.figure(figsize = (12, 5)) #创建画板，并设置画板长宽
c = '#366DE8'

plt.subplot(121) #在第一个画板的第一子图上作图
cdf = ECDF(data['normalized-losses']) # ECDF：累计连续分布函数
plt.plot(cdf.x, cdf.y, label = "statmodels", color = c)
plt.xlabel('normalized losses')
plt.ylabel('ECDF')

plt.subplot(122)  #在第一个画板的第二子图上作图
plt.hist(data['normalized-losses'].dropna(), 
         bins = int(np.sqrt(len(data['normalized-losses']))), #np.sqrt(x):对x数组进行平方根计算
         color = c)

由图可见，大部分数据集中在75到200中，数据分布不均衡，如果用均值填充缺失值是不适合的。这个特征跟哪些因素可能有关呢？应该是和保险风险等级有关吧，所以我们可以按'symboling'分组来，用相应等级的平均值进行填充这样会更精确一些。

data.groupby('symboling')['normalized-losses'].describe()#对'symboling'分组，计算不同组的'normalized-losses'平均值

# 对于其他列的缺失值是比较少的，可用dropna去掉
data = data.dropna(subset = ['price', 'bore', 'stroke', 'peak-rpm', 'horsepower', 'num-of-doors'])
#按保险风险等级分组得到平均值并对normalized-losses列的缺失值进行填充
data['normalized-losses'] = data.groupby('symboling')['normalized-losses'].transform(lambda x:x.fillna(x.mean()))
data.head()

2.2 特征相关性

缺失值处理后，我们还需要分析特征之间存在的相关性，如果相关性趋向于1，说明这两个特征是一回事，可以去掉其中一个特征，避免出现多重共线，除此我们还需要分析特征之间是否存在一些逻辑关系，比如：数据中有长、宽、高三个特征，我们可以将其合并相乘得一个体积变量。

cormatrix = data.corr()#corr:得到相关矩阵，每个值是各个列的相关系数，如：symboling与normalized-losses列的相关系数是0.593658。
cormatrix *= np.tri(*cormatrix.values.shape, k=-1).T  #返回函数的上三角矩阵，把对角线上的置0，让他们不是最高的。
cormatrix = cormatrix.stack()#在用pandas进行数据重排时，经常用到stack和unstack两个函数。stack:以列为索引堆积，unstack:以行为索引展开。
cormatrix = cormatrix.reindex(cormatrix.abs().sort_values(ascending=False).index).reset_index()
'''reindex(新索引)：按新索引排序；
abs():返回绝对值；
sort_values()：排序，ascending=False：升序，默认true：降序；
reset_index()：将行索引转为新列的值，并命名level_'''
cormatrix.columns = ["FirstVariable", "SecondVariable", "Correlation"]
cormatrix.head(10)

从左图：'city_mpg' 和 'highway-mpg'可保留一个特征，避免产生多重共线性。

data['volume'] = data.length * data.width * data.height
data.drop(['width', 'length', 'height', 'curb-weight', 'city-mpg'], axis = 1,inplace = True) # axis=1 删除columns

我们也可作图看特征之间的线性关系：

sns.pairplot(data, hue = 'fuel-type', palette = 'plasma')
#hue:分类体现，例子分了两类，紫色烧汽油，橙色烧天然气，palette：设置颜色

让我们仔细看看价格和马力变量之间的关系

lmplot:是专门做回归分析作图。公式：sns.lmplot(横轴列名，纵轴列名，数据变量，hue=类型列，col=列，row=列，palette=作图风格， fit_reg=是否画阴影及回归线) hue：子图，col：横轴子图对比，row：纵轴子图对比

sns.lmplot('price', 'horsepower', data, 
           hue = 'fuel-type', col = 'fuel-type',  row = 'num-of-doors', 
           palette = 'plasma', 
           fit_reg = True);

2.3 预处理（标准化，分类处理）

对连续值进行标准化（或归一化）

target = data.price #目标标签
regressors = [x for x in data.columns if x not in ['price']]  #得到剃掉目标变量的其他特征变量名
features = data[regressors] 
cols=list(features.loc[:,features.dtypes!=object].columns)  #筛选出数据类型不是字符型的列
features[cols]=StandardScaler().fit_transform(features[cols]) #进行标准化处理
features.head()

对分类属性就行one-hot encoding（或称哑变量，独热编码）即对非数值类型的字符进行分类转换成数字。

classes=[x for x in features.columns if x not in cols]  #筛选出分类变量列名
dummies = pd.get_dummies(features[classes])  #将数据转化成独热编码
#将分类处理后的数据添加到数据表中，并删除处理前的数据
features = features.join(dummies).drop(classes, axis = 1)
features.head()

3.Lasso 回归

多加了一个绝对值项来惩罚过大的系数（(预测值Y-实际值Y)/N-入|w|），alphas参数即是入|w|，惩罚力度越大，alphas值越大，若alphas=0，即没有惩罚项，就等于线性回归，即最小二乘法。

#将数据划分出训练集和测试集
X_train, X_test, y_train, y_test = train_test_split(features, target, test_size = 0.3,random_state = seed)  
#features：处理后数据集  target：目标标签  test_size：测试集占数据集的3成  random_state：划分训练集和测试集机器会随机从数据集中提取，而random_state=seed，指定每次随机方式都是一样的

#选择alphas值的方法有两个：
#方法一：建一个for循环，测试alphas等于多少时方程最优
#设定不同的参数值
alphas = 2.** np.arange(2, 12)# A**B =A^B
scores = np.empty_like(alphas)#empty_like:
#用不同的alphas值做模型
#enumerate() 函数用于将一个可遍历的数据对象(如列表、元组或字符串)组合为一个索引序列，同时列出数据和数据下标，一般用在 for 循环当中
for i, a in enumerate(alphas):#i:alphas索引序列，a：alphas数据
    lasso = Lasso(random_state = seed)  #指定lasso模型
    lasso.set_params(alpha = a) #确定alpha值
    lasso.fit(X_train, y_train)  #.fit(x,y) 执行
    scores[i] = lasso.score(X_test, y_test)  #用测试值测算公式

#方法二：常用交叉验证
#lassocv：交叉验证模型，cv：交叉验证切分个数 
lassocv = LassoCV(cv = 10, random_state = seed)#制定模型，将训练集平均切10分，9份用来做训练，1份用来做验证，可设置alphas=[]是多少（序列格式），默认不设置则找适合训练集最优alpha
lassocv.fit(X_train, y_train) #用训练集训练模型
lassocv_score = lassocv.score(X_test, y_test) #用测试集测试模型，返回R^2
lassocv_alpha = lassocv.alpha_  
print(lassocv_alpha,lassocv_score)

#Out[1]: 17.429955242800645 0.827600976449054

#作图
plt.figure(figsize = (10, 4))
plt.plot(alphas, scores, '-ko')
plt.axhline(lassocv_score, color = c)
plt.xlabel(r'$\alpha$')
plt.ylabel('CV Score')
plt.xscale('log', basex = 2)
sns.despine(offset = 15)

coefs = pd.Series(lassocv.coef_, index = features.columns)# 返回变量系数（lassocv）
#lasso挑选了24个特征，排除了其他42个特征。 
print("Lasso picked " + str(sum(coefs != 0)) + " features and eliminated the other "+ str(sum(coefs == 0)) + " features.")

coefs = pd.concat([coefs.sort_values().head(5), coefs.sort_values().tail(5)]) #筛选出变量系数值前5和后5 #concat()：连接函数
#作图，将前5个值和后5个值画出来
plt.figure(figsize = (10, 4))
coefs.plot(kind = "barh", color = c)
plt.title("Coefficients in the Lasso Model")
plt.show()

4.预测及检验

model=LassoCV(alphas=alphas,cv=10,random_state=seed).fit(X_train,y_train) #调用lassocv函数，并用测试集训练模型
y_pred=model.predict(X_test)#预测值
model.score(X_test,y_test) #返回判定系数R^2

#Out[1]:0.8302697453896324

plt.rcParams['figure.figsize'] = (6.0, 6.0) #设置绘图大小
#建dataframe两列，用模型预测出的预测值及测试集的观察值
preds = pd.DataFrame({"preds": model_l1.predict(X_test), "true": y_test})
preds["residuals"] = preds["true"] - preds["preds"] #增加一列，观察值-预测值（差值）
preds.plot(x = "preds", y = "residuals", kind = "scatter", color = c) #画图

#用MSE和R^2检验预测值和观察值
def MSE(y_true,y_pred):
    #mean_squared_error：均方误差MSE，反映估计量与被估计量之间差异程度的一种度量，我们希望MSE越小越好
    mse = mean_squared_error(y_true, y_pred)
    print('MSE: %2.3f' % mse)
    return mse

def R2(y_true,y_pred): 
    #r2_score:R^2
    r2 = r2_score(y_true, y_pred)
    print('R2: %2.3f' % r2)     
    return r2

MSE(y_test, y_pred); R2(y_test, y_pred);

#Out[1]:
#MSE: 3935577.218
#R2: 0.830

用MES和R方检验预测值和观察值，（个人觉得MSE偏大）不过返回 $R^2$ =0.83结果还是喜人的。这个实例就到此咯，后续如果发现什么问题再更。