线性回归、Lasso回归、岭回归预测北京PM2.5浓度
一、项目背景
1.数据
数据主要包括2010年1月1日至2014年12月31日间北京pm2.5指数以及相关天气指数数据。
数据分为训练数据和测试数据,分别保存在pm25_train.csv和pm25_test.csv两个文件中。
其中训练数据主要包括35746条记录,13个字段,主要字段说明如下:
- date:观测数据发生的日期(年-月-日)
- hour:观测数据发生的时间点(时)
- pm2.5:观测时间点对应的pm2.5指数((ug/m^3)
- DEWP:露点,空气中水气含量达到饱和的气温(℃)
- TEMP:温度,观测时间点对应的温度(℃)
- PRES:压强,观测时间点对应的压强(hPa)
- Iws:累积风速,观测时间点对应的累积风速(m/s)
- Is:累计降雪,到观测时间点为止累计降雪的时长(小时)
- Ir:累计降雨,到观测时间点为止累计降雨的时长(小时)cbwd_NE:观测时间点对应的风向为东北风(m/s)
- cbwd_NW:观测时间点对应的风向为西北风(m/s)
- cbwd_SE:观测时间点对应的风向为东南风(m/s)
- cbwd_cv:观测时间点对应的风向为静风(m/s)
测试数据主要包括6011条记录,12个字段,测试数据的字段信息和训练数据相比,除了不包括pm2.5字段以外其他完全相同。学员需要通过所学的知识,利用训练数据建立回归模型,并用于预测测试数据相应的pm2.5指数。
2.评分算法
算法通过计算平均预测误差来衡量回归模型的优劣。平均预测误差越小,说明回归模型越好。平均预测误差计算公式如下:
其中,mse是平均预测误差,m是测试数据的记录数(即6011), 是参赛者提交的PM2.5预测指数,y是对应的真实PM2.5指数。
二、数据预处理
1.读取数据
[1]:import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
[2]:train_data = pd.read_csv("pm25_train.csv")
test_data = pd.read_csv('pm25_test.csv')
[3]:train_data.info()
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 35746 entries, 0 to 35745
Data columns (total 13 columns):
date 35746 non-null object
hour 35746 non-null int64
pm2.5 35746 non-null float64
DEWP 35746 non-null int64
TEMP 35746 non-null float64
PRES 35746 non-null float64
Iws 35746 non-null float64
Is 35746 non-null int64
Ir 35746 non-null int64
cbwd_NE 35746 non-null int64
cbwd_NW 35746 non-null int64
cbwd_SE 35746 non-null int64
cbwd_cv 35746 non-null int64
dtypes: float64(4), int64(8), object(1)
memory usage: 3.5+ MB
[4]:train_data.head(3)
date hour pm2.5 DEWP TEMP PRES Iws Is Ir cbwd_NE cbwd_NW cbwd_SE cbwd_cv
0 2010-01-02 0 129.0 -16 -4.0 1020.0 1.79 0 0 0 0 1 0
1 2010-01-02 1 148.0 -15 -4.0 1020.0 2.68 0 0 0 0 1 0
2 2010-01-02 2 159.0 -11 -5.0 1021.0 3.57 0 0 0 0 1 0
可以看出数据是经过处理的,没有空值,非常干净。但其中date是object,所以需要将其转化为数值型数据。这里将其分解为年、月、日、星期。
2.数据处理
将object型数据转化为数值型数据
[5]:train_data["year"] = pd.to_datetime(train_data["date"]).dt.year # 年
train_data["month"] = pd.to_datetime(train_data["date"]).dt.month # 月
train_data["day"] = pd.to_datetime(train_data["date"]).dt.day # 日
train_data["weekday"] = pd.to_datetime(train_data["date"]).dt.weekday # 星期
test_data["year"] = pd.to_datetime(test_data["date"]).dt.year
test_data["month"] = pd.to_datetime(test_data["date"]).dt.month
test_data["day"] = pd.to_datetime(test_data["date"]).dt.day
test_data["weekday"] = pd.to_datetime(test_data["date"]).dt.weekday
3.正态检验
在做多元线性回归前,需要对因变量进行正态检验。对于某些非正态的连续变量,需要进行处理,比如对数化等。
此处借鉴机器学习实践系列(三)----达观杯–北京PM2.5浓度回归分析训练赛
[6]:from scipy.stats import norm
import seaborn as sns
sns.distplot(train_data["pm2.5"],fit=norm)
[7]:from scipy import stats
stats.probplot(train_data['pm2.5'],plot=plt)
可以看出因变量并不符合正态分布,有严重的右偏。probplot()可以绘制数据的概率图,直观的检查数据是否满足某种特性的分布。可以看出并不满足线性分布。于是尝试将因变量对数化:
[8]:train_data.drop(train_data[train_data["pm2.5"] == 0].index,inplace=True)
[9]:sns.distplot(np.log(train_data["pm2.5"]),fit=norm)
[10]:stats.probplot(np.log(train_data['pm2.5']),plot=plt)
对图像进行比较我们可以发现,对数化之后的数据更符合线性分布。所以我们添加一列对数化pm2.5.
[11]:train_data["log_pm2.5"] = np.log(train_data['pm2.5'])
三、建模
1.普通线性回归
#引包
[12]:from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error,r2_score
# 特征标签分离
[13]:features = [ 'hour', 'DEWP', 'TEMP', 'PRES', 'Iws', 'Is', 'Ir',
'cbwd_NE', 'cbwd_NW', 'cbwd_SE', 'cbwd_cv', 'year', 'month', 'day',
'weekday']
labels = "log_pm2.5"
# 划分训练集与测试集
[14]:X_train,X_test,Y_train,Y_test = train_test_split(train_data[features],train_data[labels],test_size=0.2,random_state = 10)
# 训练模型
[15]:model = LinearRegression()
model.fit(X_train,Y_train)
pre = model.predict(X_test)
# 检测误差
[16]:mean_squared_error(pre,Y_test) #MSE
0.6169782885037972
[17]:r2_score(Y_test,pre) #拟合度
0.4214605730380331
初步来看拟合度并不高,均方误差虽然在对数化情况下很小,但是还原到真实数值依旧很大。所以考虑使用L1正则化和L2正则化试一下。
2.Lasso回归(L1正则化)
[18]:from sklearn.linear_model import Lasso,LassoCV
# 通过交叉验证找到最优λ
[19]:Lambdas=np.logspace(-10,2,2000)
lasso_cv=LassoCV(alphas=Lambdas,normalize=True,cv=10,max_iter=10000)
lasso_cv.fit(X_train,Y_train)
lasso_cv.alpha_
1.86911411820748e-07
[20]:lasso=Lasso(alpha=lasso_cv.alpha_,normalize=True,max_iter=10000)
lasso.fit(X_train,Y_train)
lasso_pre=lasso.predict(X_test)
[21]:mean_squared_error(Y_test,lasso_pre)
0.616976900661324
[22]:r2_score(Y_test,lasso_pre)
0.42146187441541394
Lasso回归对模型影响并不大,MSE只是减少了0.000002,基本可以忽略。但是Lasso回归可以帮我们找出对预测结果影响不大的特征,如cbwd_cv为0,所以基本对预测值影响不大。(因为Lasso回归会下降过快到0)
[23]:print(pd.Series(index=['Intercept']+X_train.columns.tolist(),
data=[lasso.intercept_]+lasso.coef_.tolist()))
Intercept 55.253721
hour 0.012310
DEWP 0.055657
TEMP -0.076257
PRES -0.020764
Iws -0.003627
Is -0.026637
Ir -0.079902
cbwd_NE -0.394258
cbwd_NW -0.447996
cbwd_SE 0.136994
cbwd_cv 0.000000
year -0.014454
month -0.016871
day 0.006179
weekday 0.013773
dtype: float64
3.岭回归(L2正则化)
# 引包
[24]:from sklearn.linear_model import Ridge,RidgeCV
# 通过交叉验证找到最优λ
[25]:lambdas = np.logspace(-5,2,200)
[26]:ridge_cv = RidgeCV(alphas=lambdas,normalize=True,scoring="neg_mean_squared_error",cv=10)
ridge_cv.fit(X_train,Y_train)
ridge_cv.alpha_
5.478901179593945e-05
# 使用最优
[27]:model=Ridge(alpha=ridge_cv.alpha_,normalize=True)
model.fit(X_train,Y_train)
ridge_pre = model.predict(X_test)
[28]:mean_squared_error(Y_test,ridge_pre)
0.6169759952335254
[29]:r2_score(Y_test,ridge_pre)
0.4214627234334256
误差依旧没有显著减少,所以可能需要使用其他方法如PCA降维,或者其他模型会更好一些。再或者根据Lasso回归给出的回归系数,有选择的剔除一些特征,也可以,待验证。
四、总结
本文分别使用线性回归,Lasso回归和岭回归对北京PM2.5进行回归预测。虽然预测并能说是准确,但是对于回归模型的具体使用有了更深的了解。希望以后会使用其他模型解决这个问题!