AQI指的是空气质量指数，用来衡量一个城市的空气清洁或污染的程度，数值越小则空气质量越好。
怎么来预测一个城市的空气质量？【线性回归】
现在获取了2015年的空气质量数据集，该数据集包含了全国主要城市的相关数据及空气质量指数，数据情况如下： 列名含义City城市名AQI空气质量指数Precipitation降雨量GDP城市生产总值Temperature温度Longitude纬度Latitude经度Altitude海拔高度PopulationDensity人口密度Coastal是否沿海GreenCoverageRate绿化覆盖率Incineration(10000ton)焚烧量(10000吨)
data = pd.read_csv("data.csv") data.drop_duplicates(inplace=True) 最好的5个城市 t = data[["City", "AQI"]].sort_values("AQI") display(t.iloc[:5]) plt.xticks(rotation=45) # x 轴刻度旋转度数（与水平方向的逆时针夹角度数） sns.barplot(x="City", y="AQI", data=t.iloc[:5]) 最差的5个城市 t = data[["City", "AQI"]].sort_values("AQI", ascending=False) display(t.iloc[:5]) plt.xticks(rotation=45) sns.barplot(x="City", y="AQI", data=t.iloc[:5])
对于AQI，可以对空气质量进行等级划分，划分标准如下： AQI指数等级描述0-50一级优51-100二级良101-150三级轻度污染151-200四级中度污染201-300五级重度污染>300六级严重污染 根据该标准，来统计下全国空气质量每个等级的数量。
## 将AQI转换为对应的等级 def value_to_level(AQI): if AQI >= 0 and AQI <= 50: return "一级" elif AQI >= 51 and AQI <= 100: return "二级" elif AQI >= 101 and AQI <= 150: return "三级" elif AQI >= 151 and AQI <= 200: return "四级" elif AQI >= 201 and AQI <= 300: return "五级" else: return "六级" level = data["AQI"].apply(value_to_level) display(level.value_counts()) sns.countplot(x=level, order=["一级", "二级", "三级", "四级", "五级"])
可见，城市空气质量主要以一级(优)与二级(良)为主，三级(轻度污染)占一部分，更高污染的城市占少数。
接下来绘制下全国各城市的空气质量分布图。
## hue：分组依据的类别变量 palette：在对数据进行分组时，设置不同组数据的显示颜色 RdYlGn_r：红黄绿的渐变颜色，r是reverse的意思 sns.scatterplot(x="Longitude", y="Latitude", hue="AQI", palette=plt.cm.RdYlGn_r, data=data)
从结果上可以发现，从大致的地理位置，西部城市要好于东部城市，南部城市要好于北部城市。
先来统计下临海城市与内陆城市的数量：
display(data["Coastal"].value_counts()) sns.countplot(x="Coastal", data=data)
然后来观察临海城市与内陆城市的散点分布：
sns.stripplot(x="Coastal", y="AQI", data=data) # 分类散点图 sns.swarmplot(x="Coastal", y="AQI", data=data) # 分簇散点图（蜂群图）
然后分组计算空气质量的均值：
display(data.groupby("Coastal")["AQI"].mean()) sns.barplot(x="Coastal", y="AQI", data=data) # 柱状图
在柱状图（图中均值的竖线是均值的置信区间）中，仅显示了内陆城市与沿海城市的空气质量指数（AQI）的均值对比，可以用箱线图来显示更多的信息。
sns.boxplot(x="Coastal", y="AQI", data=data)
也可以绘制小提琴图，除了能够展示箱线图的信息外，还能呈现出分布的密度。
sns.violinplot(x="Coastal", y="AQI", data=data)
也可以同时绘制蜂群图和小提琴图：
sns.violinplot(x="Coastal", y="AQI", data=data, inner=None) sns.swarmplot(x="Coastal", y="AQI", color="g", data=data)
为了进一步验证结果，可以使用两独立样本t检验，来查看临海城市与内陆城市的均值差异是否显著。
在统计学中，两独立样本t检验的前提条件为:
同时，还需要考虑方差是否相同，方差是否相同，检验结果也会有所差异。
想判断一组数据是否来自于正态分布的总体，可以采用以下的几个方式进行验证：
首先，绘制两个样本的分布图：
from scipy import stats coastal = data[data["Coastal"] == 1]["AQI"] inland = data[data["Coastal"] == 0]["AQI"] fig, ax = plt.subplots(1, 2) fig.set_size_inches(15, 5) sns.distplot(coastal, ax=ax[0]) sns.distplot(inland, ax=ax[1])
PP图(Probability-Probability plot)与QQ图(Quantle-Quantile plot)本质上基本是相同的。用于检验样本数据的分布是否符合某个分布(默认为正态分布)。
import statsmodels.api as sm def plot_pp_qq(d): '''绘制PP图与QQ图 d: 要绘制的数值，array-like ''' fig,ax=plt.subplots(1,2) fig.set_size_inches(15,5) scale_data=(d-d.mean())/d.std() p=sm.ProbPlot(data=scale_data,dist=stats.norm) p.ppplot(line='45',ax=ax[0]) ax[0].set_title('PP图') p.qqplot(line='45',ax=ax[1]) ax[1].set_title('QQ图') plt.show() plot_pp_qq(coastal) plot_pp_qq(inland)
从结果上可以发现，不太可能服从正态分布。
## 原假设，观测值来自于正态分布的总体 print(stats.normaltest(coastal)) print(stats.normaltest(inland))
P<0.05，拒绝原假设，说明两个样本都不是来自于正态分布的总体。
经过一系列的周折，结果令人失望，两样本同时服从正态分布的可能性几乎为0。这似乎意味着之前所作出的一切努力付诸东流。
在假设检验中，很多假设都是以总体为正态分布为前提的，例如Z检验与t检验。当总体并非正态分布时，可以使用非参数检验。
可以使用非参数检验，来实现两样本是否来自同一分布。
## 曼-惠特尼检验。应该仅在每个停本答量 > 20时使用 原假设，两个样本服从相同的分布 print(stats.mannwhitneyu(coastal, inland, alternative="two-sided")) # 默认双边检验 威尔科克森秩和检验 原候设，两个样本服从相同的分布 print(stats.ranksums(coastal, inland))
P<0.05，拒绝原假设，说明两个样本是服从不同的分布的，意味着它们相应的均值中位数是不同的，也就意味着临海城市与内陆城市它们的空气质量并非是相同的。
可以绘制散点图矩阵，该矩阵会显示两两变量之间的关系。因为变量较多，图像会变得很小，故这里仅选择部分变量进行绘制。
sns.pairplot(data[["AQI", "PopulationDensity", "GreenCoverageRate"]])
协方差体现的是两个变量之间的离散程度以及两个变量变化步调是否一致。
相关系数，用来体现两个数值变量之间的相关性，是协方差的标准化版本，最为常用的为皮尔逊相关系数，衡量两个数值变量的线性相关性。
相关系数的取值范围为[-1,1]，可以根据相关系数的取值(绝对值)来衡量两个变量的相关性：
0.0-0.2：极弱相关或无相关
x = data["AQI"] y = data["Precipitation"] 计算协方差 a = (x - x.mean()) * (y - y.mean()) cov = np.sum(a) / (len(a) - 1) print("协方差：", cov) 计算相关系数 corr = cov / np.sqrt(x.var() * y.var()) print("相关系数：", corr) series对象也提供了cov和corr方法 print("协方差：", x.cov(y)) print("相关系数：", x.corr(y)) dataframe对象也提供了corr方法，可计算任意两个数值变量的相关系数 data.corr() 用热图来展示相关系数 plt.figure(figsize=(15, 10)) annot：是否标记数值 ax = sns.heatmap(data.corr(), cmap=plt.cm.RdYlGn, annot=True, fmt=".2f")
AQI跟纬度有0.55的正相关，也就是说，随着纬度的增大，AQI也是增大的，即越北空气质量越差，说明南方城市的空气质量好一些。
从结果中可知，空气质量指数主要受降雨量(-0.40)与纬度(0.55)影响。
此外，还能够发现其他一些明显的细节:
Altitude (海拔)与Precipitation (降雨量)负相关(-0.32)
网上有传闻说，全国所有城市的空气质量指数均值在71左右，请问，这个消息可靠嘛？
data["AQI"].mean() r = stats.ttest_1samp(data["AQI"], 71) # 单样本t检验 print("t值：", r.statistic) print("p值：", r.pvalue)
P>0.05，故无法拒绝原假设，因此接受原假设。
要清楚，+-1.96倍的标准差，是正态分布在置信区间为95%下的临界值，严格来说，对 t 分布不是如此，只不过，当样本容量较大时，t分布近似正态分布。但当样本容量较小时，二者会有较大的差异。可以获取更准确的置信区间：
mean = data["AQI"].mean() std = data["AQI"].std() stats.t.interval(0.95, df=len(data) - 1, loc=mean, scale=std / np.sqrt(len(data)))
由此，计算出全国所有城市平均空气质量指数，有95%的可能在70.63-80.04之间。
是否能够预测该城市的空气质量指数呢？答案是肯定的。可以通过对以往的数据，去建立一种模式，然后将这种模式去应用于未知的数据，进而预测结果。
对于模型来说，内部进行的都是数学上的运算，在进行建模之前，我们需要首先进行转换，将类别变量转换为离散变量。
data["Coastal"] = data["Coastal"].map({"是": 1, "否": 0}) data["Coastal"].value_counts()
Coastal变量是一个二分类变量，在对变量进行转换时，可以将“是”与“否”转换为其他离散整数值，例如1和0。
假设某变量是一个多分类变量(三个或是三个以上的分类)，在对变量进行转换时，不可以将n个类别映射为0 1 2 ... n-1。
这时，可以采用独热编码即 One-Hot 编码，又称一位有效编码，其方法是把一个特征拆分为多个特征。 x→总经理部门经理普通员工总经理100部门经理010普通员工001
首先，不进行任何处理，建立一个基模型，作为基线（baseline）。后续的操作，可以在此基础上进行改进。
from sklearn.linear_model import LinearRegression from sklearn.model_selection import train_test_split city（城市名称）对预测毫无用处，删掉 x = data.drop(["City", "AQI"], axis=1) y = data["AQI"] x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.3, random_state=0) lr = LinearRegression() lr.fit(x_train, y_train) print(lr.score(x_train, y_train)) print(lr.score(x_test, y_test)) y_hat = lr.predict(x_test) plt.figure(figsize=(15, 5)) plt.plot(y_test.values, "-r", label="真实值", marker="o") plt.plot(y_hat, "-g", label="测试值", marker="D") plt.legend(loc="upper left") plt.title("线性回归预测结果", fontsize=20)
可以依据箱线图判断离群点的原则去探索异常值，然后使用临界值替换掉异常值。
## Coastal是类别变量，映射为离散变量，不会有异常值 for col in x.columns.drop("Coastal"): if pd.api.types.is_numeric_dtype(x_train[col]): quartile = np.quantile(x_train[col], [0.25, 0.75]) IQR = quartile[1] - quartile[0] lower = quartile[0] - 1.5 * IQR upper = quartile[1] + 1.5 * IQR x_train[col][x_train[col] < lower] = lower x_train[col][x_train[col] > upper] = upper x_test[col][x_test[col] < lower] = lower x_test[col][x_test[col] > upper] = upper
去除异常值后，用新的训练集与测试集来评估模型的效果。
lr.fit(x_train, y_train) print(lr.score(x_train, y_train)) print(lr.score(x_test, y_test))
结果跟之前相比有了一定的少许提升。
数据标准化不会对线性回归的评估效果造成影响，不过，会影响模型的权重系数。如果需要考虑线性回归的权重，则需要对数据首先进行标准化处理。
from sklearn.preprocessing import StandardScaler s = StandardScaler() x_train_scale = s.fit_transform(x_train) x_test_scale = s.transform(x_test)
上面我们使用所有可能的原始数据作为特征，建立模型，然而，特征并非越多越好，有些对模型质量没有什么改善的特征，可以进行删除，以提高模型质量。
特征选择的方式有很多，这里使用RFECV方法来实现特征选择。RFECV分成两个部分，分别是:
具体过程如下：
RFE阶段
CV阶段
确定平均分最高的特征数量，完成特征选择
from sklearn.feature_selection import RFECV estimator：要操作的模型 step：每次删除的变量数 cv：使用的交叉验证折数（就是把数据分成几份） scoring：评估方式 rfecv = RFECV(estimator=lr, step=1, cv=5, scoring="r2") rfecv.fit(x_train, y_train) 返回经过选择之后，剩余的特征数量 print(rfecv.n_features_) 返回经过特征选择后，使用缩减特征训练后的模型 print(rfecv.estimator_) 返回每个特征的等级，数值越小，特征越重要 print(rfecv.ranking_) 返回布尔数组，用来表示特征是否被选择 print(rfecv.support_) 返回对应数量特征时，模型交叉验证的评分 print(rfecv.cv_results_) 绘制特征选择过程中，使用交叉验证获取的R^2 plt.plot(range(1, len(rfecv.cv_results_["mean_test_score"]) + 1), rfecv.cv_results_["mean_test_score"], marker="o") plt.xlabel("特征数量") plt.ylabel("交叉验证R^2值")
然后，对测试集应用这种特征选择（变换），进行测试，获取测试结果。
print("剔除的变量：", x_train.columns[~rfecv.support_]) x_train_eli = rfecv.transform(x_train) x_test_eli = rfecv.transform(x_test) print(rfecv.estimator_.score(x_train_eli, y_train)) print(rfecv.estimator_.score(x_test_eli, y_test))
可以发现，经过特征选择后，使用剩余特征训练的模型，与之前模型表现上几乎没有变化，也就是说，我们成功删除了一些特征。
## 获取列名与对应的权重，构成一个元组，作为列表的元素 li = list(zip(x.columns.values[rfecv.support_], rfecv.estimator_.coef_)) 根据权重的绝对值，对列表进行降序排列 li.sort(key=lambda x: abs(x[1]), reverse=True) s = pd.Series(dict(li)) plt.figure(figsize=(15, 5)) ax = sns.barplot(y=s.index, x=s.values) for y, x in enumerate(s): 在条形图上绘制数值 t = ax.text(x / 2, y, round(x, 3)) 教值居中对齐 t.set_ha("center") plt.show()
分箱后，我们不能将每个区间映射为离散数值，而是应当使用One-Hot编码。
from sklearn.preprocessing import KBinsDiscretizer KBinsDiscretizer K个分箱的离散器，用于将数值（通常是连续变量）变量进行区间离散化操作 n_bin：分箱（区间）个数 encode：离散化编码方式，分为：onehot、onehot-dense与ordinal onehot：使用独热编码，返回稀疏矩阵 onehot-dense：使用独热编码，返回稠密矩阵 ordinal：使用序数编码（0,1,2……） strategy：分箱的方式，分为：uniform、quantile、kmeans uniform：每个区间的长度范围大致相同 quantile：每个区间包含的元素个数大致相同 kmeans：使用一维kmeans方式进行分箱 k = KBinsDiscretizer(n_bins=[4, 5, 10, 6], encode="onehot-dense", strategy="uniform") 定义离散化的特征 discretize = ["Longitude", "Temperature", "Precipitation", "Latitude"] x_train_eli = pd.DataFrame(data=x_train_eli, columns=x.columns[rfecv.support_]) x_test_eli = pd.DataFrame(data=x_test_eli, columns=x.columns[rfecv.support_]) r = k.fit_transform(x_train_eli[discretize]) r = pd.DataFrame(r, index=x_train_eli.index) 获取除离散化特征之外的其他特征 x_train_dis = x_train_eli.drop(discretize, axis=1) 将离散化后的特征与其他特征进行重新组合 x_train_dis = pd.concat([x_train_dis, r], axis=1) 对测试集进行同样的离散化操作 r = pd.DataFrame(k.transform(x_test_eli[discretize]), index=x_test_eli.index) x_test_dis = x_test_eli.drop(discretize, axis=1) x_test_dis = pd.concat([x_test_dis, r], axis=1) 查看转换之后的格式 display(x_train_dis.head())
本文作者:a
本文链接:
版权声明:本博客所有文章除特别声明外，均采用 BY-NC-SA 许可协议。转载请注明出处！