线性回归、Lasso回归、岭回归预测北京PM2.5浓度

一、项目背景

北京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进行回归预测。虽然预测并能说是准确，但是对于回归模型的具体使用有了更深的了解。希望以后会使用其他模型解决这个问题！
数学推导：深入理解线性回归与岭回归(数学推导)