AQI分析與預(yù)測(cè)-2

一.項(xiàng)目背景

AQI（air Quality Index）指空氣質(zhì)量指數(shù)，用來(lái)衡量空氣清潔或者污染程度。值
越小，表示空氣質(zhì)量越好。近年來(lái)因?yàn)榄h(huán)境問(wèn)題，空氣質(zhì)量越來(lái)越受到人們重視。

上篇文章[AQI分析與預(yù)測(cè)（一）](http://mp.weixin.qq.com/s?__biz=MzIzODI4ODM2MA==&mid=2247486525&idx=1&sn=a92d40e0d715f9ab63563b4e8b054a0d&chksm=e93ae0bade4d69ac57407216484cdbb6d7a8a812baef1c0e4ab0d5aeae0c7f0aff9fd274127b&scene=21#wechat_redirect)我們進(jìn)行了初步分析，主要分析了空氣質(zhì)量最好/差城市和臨海城市是否空氣質(zhì)量?jī)?yōu)于內(nèi)陸城市這兩個(gè)問(wèn)題，本篇我們?cè)谥盎A(chǔ)上繼續(xù)研究如下問(wèn)題。

1.空氣質(zhì)量受那些因素影響
2.關(guān)于空氣質(zhì)量驗(yàn)證
3.構(gòu)建空氣質(zhì)量預(yù)測(cè)模型

二.實(shí)現(xiàn)過(guò)程

1.空氣質(zhì)量受那些因素影響
  》》指標(biāo)：協(xié)方差和相關(guān)系數(shù)
  》》圖形：熱力圖

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import warnings


sns.set(style="darkgrid", font_scale=1.2)
plt.rcParams["font.family"] = "SimHei"
plt.rcParams["axes.unicode_minus"] = False
warnings.filterwarnings("ignore")


#讀取文件
data = pd.read_csv("data.csv")
#查看數(shù)據(jù)
data.head()

#kind：繪制圖像的類(lèi)型。可選值：
#scatter：散點(diǎn)圖（默認(rèn)值）。
#reg：帶有回歸線的散點(diǎn)圖。
#vars：顯示哪些變量之間的兩兩關(guān)系，默認(rèn)為顯示所有變量。
sns.pairplot(data, vars=["AQI", "PopulationDensity", "GreenCoverageRate"])

#計(jì)算相關(guān)數(shù)據(jù)
x = data["AQI"]
y = data["Precipitation"]
#計(jì)算AQI與Precipitation的協(xié)方差。
a = (x - x.mean()) * (y - y.mean())
#計(jì)算協(xié)方差
cov = np.sum(a) / (len(a) - 1)
print("協(xié)方差：", cov)

#計(jì)算AQI與Precipitation的相關(guān)系數(shù)。
corr = cov / np.sqrt(x.var() * y.var())
print("相關(guān)系數(shù)：", corr)

#pandas封裝了相關(guān)方法
print("協(xié)方差：", x.cov(y))
print("相關(guān)系數(shù)：", x.corr(y))

#初始化畫(huà)布
plt.figure(figsize=(15, 10))
#繪制熱力圖
ax = sns.heatmap(data.corr(),
                 cmap=plt.cm.RdYlGn, 
                 annot=True,
                 fmt=".2f")
  
 結(jié)論：降雨量越多，空氣質(zhì)量越好；維度越低，空氣質(zhì)量越好

2.關(guān)于空氣質(zhì)量驗(yàn)證
  》》問(wèn)題：全國(guó)所有城市的空氣質(zhì)量指數(shù)均值在71左右，請(qǐng)問(wèn)此結(jié)論準(zhǔn)確嗎？
  》》方法：假設(shè)檢驗(yàn)

該需求是驗(yàn)證樣本均值是否等于總體均值，根據(jù)條件，我們可以使用單樣本t檢驗(yàn)，
置信度為95%。

#進(jìn)行單樣本t檢驗(yàn)
r = stats.ttest_1samp(data["AQI"], 71)
#輸出檢驗(yàn)統(tǒng)計(jì)量
print("t值：", r.statistic)
#輸出p值
print("p值：", r.pvalue)


結(jié)論：我們可以看到P值是大于0.05的，y因此我們無(wú)法拒絕原假設(shè)，因此接受原假設(shè)

#計(jì)算均值
mean = data["AQI"].mean()
#計(jì)算標(biāo)準(zhǔn)差
std = data["AQI"].std()
#計(jì)算置信區(qū)間
stats.t.interval(0.95, df=len(data) - 1, loc=mean, scale=std / np.sqrt(len(data)))


結(jié)論：我們就計(jì)算出全國(guó)城市平均空氣質(zhì)量指數(shù)，95%的可能大致在70.63~80.04之間

3.構(gòu)建空氣質(zhì)量預(yù)測(cè)模型
  》》對(duì)于一些城市，如果能夠已知降雨量，溫度，經(jīng)緯度等指標(biāo)，我們是否能夠預(yù)測(cè)該
城市的空氣質(zhì)量指數(shù)呢？因此我們需要構(gòu)建模型，預(yù)測(cè)新的數(shù)據(jù)。
  》》過(guò)程：基模型構(gòu)建，異常值處理后構(gòu)建模型和特征選擇后構(gòu)建模型，將結(jié)果與基模
型進(jìn)行對(duì)比，看看是否進(jìn)行優(yōu)化。

#進(jìn)行類(lèi)別轉(zhuǎn)換
data["Coastal"] = data["Coastal"].map({"是": 1, "否": 0})
#統(tǒng)計(jì)類(lèi)別數(shù)目
data["Coastal"].value_counts()

from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split


#City(城市名稱(chēng))對(duì)預(yù)測(cè)毫無(wú)用處，刪掉。
X = data.drop(["City","AQI"], axis=1)
#目標(biāo)值
y = data["AQI"]
#分離測(cè)試集和訓(xùn)練集
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=0)
#構(gòu)建基礎(chǔ)線性回歸模型
lr = LinearRegression()
#訓(xùn)練模型
lr.fit(X_train, y_train)
#輸出訓(xùn)練集模型評(píng)分
print(lr.score(X_train, y_train))
#輸出測(cè)試集模型評(píng)分
print(lr.score(X_test, y_test))

#繪制圖形觀察


#預(yù)測(cè)y值
y_hat = lr.predict(X_test)
#初始化畫(huà)布
plt.figure(figsize=(15, 5))
#繪制真實(shí)值折線圖
plt.plot(y_test.values,
         "-r", 
         label="真實(shí)值",
         marker="o")
#繪制預(yù)測(cè)值折線圖
plt.plot(y_hat,
         "-g",
         label="預(yù)測(cè)值",
         marker="D")
#設(shè)置圖例
plt.legend(loc="upper left")
#設(shè)置標(biāo)題
plt.title("線性回歸預(yù)測(cè)結(jié)果", fontsize=20)

# Coastal是類(lèi)別變量，映射為離散變量，不會(huì)有異常值。


#遍歷列
for col in X.columns.drop("Coastal"):
    #對(duì)數(shù)值型數(shù)據(jù)進(jìn)行判斷
    if pd.api.types.is_numeric_dtype(X_train[col]):
        #獲取分位數(shù)
        quartile = np.quantile(X_train[col], [0.25, 0.75])
        #計(jì)算IQR
        IQR = quartile[1] - quartile[0]
        #計(jì)算正常數(shù)值下限
        lower = quartile[0] - 1.5 * IQR
        #計(jì)算正常數(shù)值上限
        upper = quartile[1] + 1.5 * IQR
        #用邊界值進(jìn)行填充異常值
        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

#訓(xùn)練模型
lr.fit(X_train, y_train)
#去除異常值后評(píng)估模型效果
print(lr.score(X_train, y_train))
print(lr.score(X_test, y_test))


結(jié)論：去除異常值后所構(gòu)建模型效果相比之前有所改進(jìn)

#對(duì)數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化處理，


from sklearn.preprocessing import StandardScaler
#數(shù)據(jù)標(biāo)準(zhǔn)化
s = StandardScaler()
#對(duì)訓(xùn)練集進(jìn)行標(biāo)準(zhǔn)化
X_train_scale = s.fit_transform(X_train)
#對(duì)測(cè)試集進(jìn)行標(biāo)準(zhǔn)化
X_test_scale = s.transform(X_test)

#對(duì)數(shù)據(jù)進(jìn)行特征選擇，目的是提高模型準(zhǔn)確率和訓(xùn)練速度


REFCV方法
  》》RFE（REcursive feature elimination）：遞歸特征消除，用來(lái)對(duì)特征進(jìn)行重要性評(píng)級(jí)
  》》CV（Cross Validation）：交叉驗(yàn)證，通過(guò)交叉驗(yàn)證，選擇最佳數(shù)量特征
具體過(guò)程如下：
RFE階段：
  1.初始特征集為所有可用特征
  2.使用當(dāng)前特征集進(jìn)行建模，然后計(jì)算每個(gè)特征的重要性
  3.刪除最不重要的一個(gè)或者多個(gè)特征，更新特征集
  4.跳到步驟2，直到完成所有特征集重要性評(píng)級(jí)
CV階段
  1.根據(jù)REF階段確定的特征重要性，依次選擇不同數(shù)量特征
  2.對(duì)選定的特征集進(jìn)行交叉驗(yàn)證
  3.確定平均分最高的特征數(shù)量，完成特征選擇。

from sklearn.feature_selection import RFECV
#estimator：要操作的模型。
#step：每次刪除的變量數(shù)。
#cv：使用的交叉驗(yàn)證折數(shù)。
#n_jobs：并發(fā)的數(shù)量。
#scoring: 評(píng)估的方式。
rfecv = RFECV(estimator=lr, 
              step=1,
              cv=5, 
              n_jobs=-1, 
              scoring="r2")
rfecv.fit(X_train_scale, y_train)
#返回經(jīng)過(guò)選擇之后，剩余的特征數(shù)量。
print(rfecv.n_features_)
#返回經(jīng)過(guò)特征選擇后，使用縮減特征訓(xùn)練后的模型。
print(rfecv.estimator_)
#返回每個(gè)特征的等級(jí)，數(shù)值越小，特征越重要。
print(rfecv.ranking_)
#返回布爾數(shù)組，用來(lái)表示特征是否被選擇。
print(rfecv.support_)
#返回對(duì)應(yīng)數(shù)量特征時(shí)，模型交叉驗(yàn)證的評(píng)分。
print(rfecv.grid_scores_)

#繪制圖形
plt.plot(range(1, len(rfecv.grid_scores_) + 1), 
         rfecv.grid_scores_,
         marker="o")
#設(shè)置x軸標(biāo)簽
plt.xlabel("特征數(shù)量")
#設(shè)置y軸標(biāo)簽
plt.ylabel("交叉驗(yàn)證$R^2$值")

print("剔除的變量：", X.columns.values[~rfecv.support_])
#應(yīng)用到訓(xùn)練集
X_train_eli = rfecv.transform(X_train_scale)
#應(yīng)用到測(cè)試集
X_test_eli = rfecv.transform(X_test_scale)
#輸出模型評(píng)分
print(rfecv.estimator_.score(X_train_eli, y_train))
print(rfecv.estimator_.score(X_test_eli, y_test))

#獲取列名與對(duì)應(yīng)的權(quán)重，構(gòu)成一個(gè)元組，作為列表的元素。
li = list(zip(X.columns.values[rfecv.support_], rfecv.estimator_.coef_))
#根據(jù)權(quán)重的絕對(duì)值，對(duì)列表進(jìn)行降序排列。
li.sort(key=lambda x: abs(x[1]), reverse=True)
#轉(zhuǎn)換為Series
s = pd.Series(dict(li))
#初始化畫(huà)布
plt.figure(figsize=(15, 5))
#繪制柱狀圖
ax = sns.barplot(y=s.index, x=s.values)
for y, x in enumerate(s):
    #繪制標(biāo)注
    t = ax.text(x / 2, y, round(x, 3))
    #設(shè)置居中對(duì)齊
    t.set_ha("center")
#顯示圖形
plt.show()

from sklearn.preprocessing import KBinsDiscretizer
#KBinsDiscretizer K個(gè)分箱的離散器。用于將數(shù)值（通常是連續(xù)變量）變量進(jìn)行區(qū)間離散化操作。
#n_bins：分箱（區(qū)間）的個(gè)數(shù)。
#encode：離散化編碼方式。分為：onehot，onehot-dense與ordinal。
#   onehot：使用獨(dú)熱編碼，返回稀疏矩陣。
#   onehot-dense：使用獨(dú)熱編碼，返回稠密矩陣。
#   ordinal：使用序數(shù)編碼（0,1,2……）。
#strategy：分箱的方式。分為：uniform，quantile，kmeans。
#uniform：每個(gè)區(qū)間的長(zhǎng)度范圍大致相同。
#quantile：每個(gè)區(qū)間包含的元素個(gè)數(shù)大致相同。
#kmeans：使用一維kmeans方式進(jìn)行分箱。


#對(duì)數(shù)據(jù)進(jìn)行分箱操作
k=KBinsDiscretizer(n_bins=[4, 5, 10, 6],
                     encode="onehot-dense", 
                     strategy="uniform")
#定義離散化的特征。
discretize=["Longitude", "Temperature", "Precipitation", "Latitude"]
#訓(xùn)練集數(shù)據(jù)轉(zhuǎn)換為DataFrame
X_train_eli=pd.DataFrame(data=X_train_eli, columns=X.columns[rfecv.support_])
#測(cè)試集數(shù)據(jù)轉(zhuǎn)換為DataFrame
X_test_eli=pd.DataFrame(data=X_test_eli, columns=X.columns[rfecv.support_])
#應(yīng)用到訓(xùn)練集
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)
#將離散化后的特征與其他特征進(jìn)行重新組合。
X_train_dis=pd.concat([X_train_dis, r], axis=1)
#對(duì)測(cè)試集進(jìn)行同樣的離散化操作。
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)
#查看轉(zhuǎn)換之后的格式。
display(X_train_dis.head()

#訓(xùn)練模型
lr.fit(X_train_dis, y_train)
#去除異常值后評(píng)估模型效果
print(lr.score(X_train_dis, y_train))
print(lr.score(X_test_dis, y_test))


結(jié)論：離散化后模型效果進(jìn)一步提升