空气质量预测¶

问题陈述：¶

我们给出了 2015 年至 2020 年印度 AQI 数据集。任务是建立一个模型，可以根据数据集中可用的独立特征准确预测 AQI。

关于数据集¶

image.png

主要空气污染物¶

颗粒物

颗粒物是固体和液体的混合物，包括碳、复杂有机化学品、硫酸盐、硝酸盐、矿物粉尘和悬浮在空气中的水。根据其尺寸（以微米为单位），它们分为 PM 2.5 和 PM 10。

氮氧化物

氮氧化物由氮和氧组成，归类为NOx气体，最常见的是一氧化氮（NO）和二氧化氮（NO2）。它们主要是化石燃料燃烧的结果。

氨

氨是大气中最丰富的碱性气体。 NH3 在大气颗粒物的形成中起着重要作用。氨被认为是农业和工业的副产品

二氧化硫

二氧化硫是一种无色、有强烈气味的气体，当煤炭和石油等含硫燃料燃烧时会形成二氧化硫，造成空气污染。当它们与大气中的物质发生反应形成酸雨时，会影响环境。

臭氧

地面臭氧是一种无色且高度刺激性的气体，形成于地球表面上方。

一氧化碳

含碳燃料的不完全燃烧会产生一氧化碳。一氧化碳是一种无色、剧毒的气体。

苯

室内的苯来自胶水、油漆、家具蜡、清洁剂等含苯产品。

甲苯

机动车和工业排放是污染物的主要来源。

二甲苯

机动车排放是城市空气环境中二甲苯的主要来源。

导入库和数据¶

In [1]:
# 导入重要的库
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
In [2]:
df_city_day = pd.read_csv('city_day.csv')
In [3]:
df_city_day.head()
Out[3]:
City Date PM2.5 PM10 NO NO2 NOx NH3 CO SO2 O3 Benzene Toluene Xylene AQI AQI_Bucket
0 Ahmedabad 2015-01-01 NaN NaN 0.92 18.22 17.15 NaN 0.92 27.64 133.36 0.00 0.02 0.00 NaN NaN
1 Ahmedabad 2015-01-02 NaN NaN 0.97 15.69 16.46 NaN 0.97 24.55 34.06 3.68 5.50 3.77 NaN NaN
2 Ahmedabad 2015-01-03 NaN NaN 17.40 19.30 29.70 NaN 17.40 29.07 30.70 6.80 16.40 2.25 NaN NaN
3 Ahmedabad 2015-01-04 NaN NaN 1.70 18.48 17.97 NaN 1.70 18.59 36.08 4.43 10.14 1.00 NaN NaN
4 Ahmedabad 2015-01-05 NaN NaN 22.10 21.42 37.76 NaN 22.10 39.33 39.31 7.01 18.89 2.78 NaN NaN

理解数据¶

In [4]:
df_city_day.shape
Out[4]:
(29531, 16)
In [5]:
df_city_day.describe()
Out[5]:
PM2.5 PM10 NO NO2 NOx NH3 CO SO2 O3 Benzene Toluene Xylene AQI
count 24933.000000 18391.000000 25949.000000 25946.000000 25346.000000 19203.000000 27472.000000 25677.000000 25509.000000 23908.000000 21490.000000 11422.000000 24850.000000
mean 67.450578 118.127103 17.574730 28.560659 32.309123 23.483476 2.248598 14.531977 34.491430 3.280840 8.700972 3.070128 166.463581
std 64.661449 90.605110 22.785846 24.474746 31.646011 25.684275 6.962884 18.133775 21.694928 15.811136 19.969164 6.323247 140.696585
min 0.040000 0.010000 0.020000 0.010000 0.000000 0.010000 0.000000 0.010000 0.010000 0.000000 0.000000 0.000000 13.000000
25% 28.820000 56.255000 5.630000 11.750000 12.820000 8.580000 0.510000 5.670000 18.860000 0.120000 0.600000 0.140000 81.000000
50% 48.570000 95.680000 9.890000 21.690000 23.520000 15.850000 0.890000 9.160000 30.840000 1.070000 2.970000 0.980000 118.000000
75% 80.590000 149.745000 19.950000 37.620000 40.127500 30.020000 1.450000 15.220000 45.570000 3.080000 9.150000 3.350000 208.000000
max 949.990000 1000.000000 390.680000 362.210000 467.630000 352.890000 175.810000 193.860000 257.730000 455.030000 454.850000 170.370000 2049.000000

数据似乎没有任何异常值

In [6]:
#数据摘要
df_city_day.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 29531 entries, 0 to 29530
Data columns (total 16 columns):
 #   Column      Non-Null Count  Dtype  
---  ------      --------------  -----  
 0   City        29531 non-null  object 
 1   Date        29531 non-null  object 
 2   PM2.5       24933 non-null  float64
 3   PM10        18391 non-null  float64
 4   NO          25949 non-null  float64
 5   NO2         25946 non-null  float64
 6   NOx         25346 non-null  float64
 7   NH3         19203 non-null  float64
 8   CO          27472 non-null  float64
 9   SO2         25677 non-null  float64
 10  O3          25509 non-null  float64
 11  Benzene     23908 non-null  float64
 12  Toluene     21490 non-null  float64
 13  Xylene      11422 non-null  float64
 14  AQI         24850 non-null  float64
 15  AQI_Bucket  24850 non-null  object 
dtypes: float64(13), object(3)
memory usage: 3.6+ MB

由于我们的日期列是对象类型，因此我们需要将其转换为日期时间格式。

In [7]:
#将日期转换为日期时间格式
df_city_day['Date']=pd.to_datetime(df_city_day['Date'])

检查空值

In [8]:
#检查空值
df_city_day.isna().sum()
Out[8]:
City              0
Date              0
PM2.5          4598
PM10          11140
NO             3582
NO2            3585
NOx            4185
NH3           10328
CO             2059
SO2            3854
O3             4022
Benzene        5623
Toluene        8041
Xylene        18109
AQI            4681
AQI_Bucket     4681
dtype: int64

因为，我们有很多空值。放弃它们并不是一个好的选择，而是我们可以用中值代替它，因为对于特定城市，特征值将遵循类似的趋势。

In [9]:
# 替换空值
null_col = ['PM2.5',	'PM10',	'NO',	'NO2',	'NOx',	'NH3',	'CO',	'SO2',	'O3',	'Benzene',	'Toluene',	'Xylene',	'AQI']
for col in null_col:
  df_city_day[col] = df_city_day[col].fillna(df_city_day[col].median())
In [10]:
df_city_day.isna().sum()
Out[10]:
City             0
Date             0
PM2.5            0
PM10             0
NO               0
NO2              0
NOx              0
NH3              0
CO               0
SO2              0
O3               0
Benzene          0
Toluene          0
Xylene           0
AQI              0
AQI_Bucket    4681
dtype: int64
In [11]:
df_city_day.head()
Out[11]:
City Date PM2.5 PM10 NO NO2 NOx NH3 CO SO2 O3 Benzene Toluene Xylene AQI AQI_Bucket
0 Ahmedabad 2015-01-01 48.57 95.68 0.92 18.22 17.15 15.85 0.92 27.64 133.36 0.00 0.02 0.00 118.0 NaN
1 Ahmedabad 2015-01-02 48.57 95.68 0.97 15.69 16.46 15.85 0.97 24.55 34.06 3.68 5.50 3.77 118.0 NaN
2 Ahmedabad 2015-01-03 48.57 95.68 17.40 19.30 29.70 15.85 17.40 29.07 30.70 6.80 16.40 2.25 118.0 NaN
3 Ahmedabad 2015-01-04 48.57 95.68 1.70 18.48 17.97 15.85 1.70 18.59 36.08 4.43 10.14 1.00 118.0 NaN
4 Ahmedabad 2015-01-05 48.57 95.68 22.10 21.42 37.76 15.85 22.10 39.33 39.31 7.01 18.89 2.78 118.0 NaN

我们仍然有 AQI_Bucket 的 NaN 值

In [12]:
# AQI_Bucket 的唯一计数
df_city_day['AQI_Bucket'].value_counts()
Out[12]:
Moderate        8829
Satisfactory    8224
Poor            2781
Very Poor       2337
Good            1341
Severe          1338
Name: AQI_Bucket, dtype: int64
In [13]:
df_city_day.groupby('AQI_Bucket')['AQI'].mean()
Out[13]:
AQI_Bucket
Good             40.268456
Moderate        136.671877
Poor            245.663430
Satisfactory     76.765686
Severe          567.886398
Very Poor       343.004279
Name: AQI, dtype: float64

以上表中的平均 AQI 为参考，我们将替换 AQI_Bucket 的 NaN 值。

In [14]:
#使用条件替换AQI_Bucket的Null值
d1 =  df_city_day[df_city_day.AQI<=50 ].AQI_Bucket.fillna('Good')
d2 =  df_city_day[(df_city_day.AQI>50) & (df_city_day.AQI<=100)].AQI_Bucket.fillna('Satisfactory')
d3 =  df_city_day[(df_city_day.AQI>100) & (df_city_day.AQI<=180)].AQI_Bucket.fillna('Moderate')
d4 =  df_city_day[(df_city_day.AQI>180) & (df_city_day.AQI<=300)].AQI_Bucket.fillna('Poor')
d5 =  df_city_day[(df_city_day.AQI>300) & (df_city_day.AQI<=450)].AQI_Bucket.fillna('Very Poor')
d6 =  df_city_day[df_city_day.AQI>450].AQI_Bucket.fillna('Severe')
df_city_day.AQI_Bucket = pd.concat([d1,d2,d3,d4,d5,d6]).sort_index()
In [15]:
df_city_day.isna().sum()
Out[15]:
City          0
Date          0
PM2.5         0
PM10          0
NO            0
NO2           0
NOx           0
NH3           0
CO            0
SO2           0
O3            0
Benzene       0
Toluene       0
Xylene        0
AQI           0
AQI_Bucket    0
dtype: int64

现在我们的数据没有任何空值。

In [16]:
df_city_day.head()
Out[16]:
City Date PM2.5 PM10 NO NO2 NOx NH3 CO SO2 O3 Benzene Toluene Xylene AQI AQI_Bucket
0 Ahmedabad 2015-01-01 48.57 95.68 0.92 18.22 17.15 15.85 0.92 27.64 133.36 0.00 0.02 0.00 118.0 Moderate
1 Ahmedabad 2015-01-02 48.57 95.68 0.97 15.69 16.46 15.85 0.97 24.55 34.06 3.68 5.50 3.77 118.0 Moderate
2 Ahmedabad 2015-01-03 48.57 95.68 17.40 19.30 29.70 15.85 17.40 29.07 30.70 6.80 16.40 2.25 118.0 Moderate
3 Ahmedabad 2015-01-04 48.57 95.68 1.70 18.48 17.97 15.85 1.70 18.59 36.08 4.43 10.14 1.00 118.0 Moderate
4 Ahmedabad 2015-01-05 48.57 95.68 22.10 21.42 37.76 15.85 22.10 39.33 39.31 7.01 18.89 2.78 118.0 Moderate

让我们做一些 EDA

探索性数据分析 (EDA)¶

In [17]:
#AQI_Bucket的分布
import seaborn as sns
sns.displot(data=df_city_day,x='AQI_Bucket')
Out[17]:
<seaborn.axisgrid.FacetGrid at 0x14a03f9e910>
  • 我们的数据不平衡。
  • 最常见的 AQI 标签是“中等”和“满意”。

让我们看看AQI_Bucket与城市的分布。

In [18]:
# AQI_Bucket 相对于城市的分布
df_city_day.groupby(['City', 'AQI_Bucket']).size().unstack().plot(kind='bar', stacked=True,figsize=(12,7))
plt.title(f"Total count of AQI_Bucket grouped by {'City'} ")
plt.ylabel('Count')
plt.show()

注意：- 此数据仅代表数据集中每个城市的数据计数。

对于大多数城市来说，与 AQI_Bucket 中的其他 AQI 标签相比，中等 AQI 的出现频率最高。

比较每个城市的平均 AQI。

In [19]:
# 比较每个城市的平均 AQI
df_city_day.groupby('City')['AQI'].mean().sort_values(ascending=False).plot(kind='bar',figsize=(12,7))
plt.ylabel('Average AQI')
plt.show()
  • 艾哈迈达巴德的AQI最差（超过300），可以归类为“非常差”。
  • 平均 AQI 超过 200 的其他城市包括德里、古尔冈、勒克瑙和巴特那。这些城市可被归类为空气质量较差的城市。

现在让我们看看影响每个城市 AQI 的污染物浓度。

In [20]:
# 各城市污染物浓度
pollutant_list = ['PM2.5', 'PM10', 'NO', 'NO2', 'NOx', 'NH3', 'CO', 'SO2', 'O3', 'Benzene', 'Toluene', 'Xylene']

# 计算需要的行列数
n = int(len(pollutant_list)**0.5) 
if n * n < len(pollutant_list):
    n += 1

# 创建一个 n x n 的画布
fig, axes = plt.subplots(n, n, figsize=(12, 12))
fig.tight_layout(pad=3.0)

# 遍历污染物列表并创建条形图
for i, col in enumerate(pollutant_list):
    ax = axes[i//n, i%n]
    df_city_day.groupby('City')[col].mean().sort_values(ascending=False).plot(kind='bar', ax=ax, fontsize=6)
    ax.set_ylabel(f'Average {col}', fontsize=12)
    ax.set_xlabel('City', fontsize=12)
    ax.set_title(f"Concentration of {col}", fontsize=12)

# 隐藏多余的图表
for j in range(i+1, n*n):
    axes[j//n, j%n].axis('off')
plt.subplots_adjust(wspace=0.5, hspace=0.8)  # 调整子图间距
plt.show()

探索历年来所有城​​市AQI的趋势。

In [21]:
df_city_day['Date'].min()
Out[21]:
Timestamp('2015-01-01 00:00:00')
In [22]:
df_city_day['Date'].max()
Out[22]:
Timestamp('2020-07-01 00:00:00')

我们的数据是从2015年1月1日到2020年7月1日

In [23]:
pd.DatetimeIndex(df_city_day['Date']).year.value_counts()
Out[23]:
2019    7446
2018    6471
2017    4689
2020    4646
2016    3478
2015    2801
Name: Date, dtype: int64
In [24]:
# 创建“年份”列
df_city_day['Year'] = pd.DatetimeIndex(df_city_day['Date']).year
In [25]:
df_city_day.head()
Out[25]:
City Date PM2.5 PM10 NO NO2 NOx NH3 CO SO2 O3 Benzene Toluene Xylene AQI AQI_Bucket Year
0 Ahmedabad 2015-01-01 48.57 95.68 0.92 18.22 17.15 15.85 0.92 27.64 133.36 0.00 0.02 0.00 118.0 Moderate 2015
1 Ahmedabad 2015-01-02 48.57 95.68 0.97 15.69 16.46 15.85 0.97 24.55 34.06 3.68 5.50 3.77 118.0 Moderate 2015
2 Ahmedabad 2015-01-03 48.57 95.68 17.40 19.30 29.70 15.85 17.40 29.07 30.70 6.80 16.40 2.25 118.0 Moderate 2015
3 Ahmedabad 2015-01-04 48.57 95.68 1.70 18.48 17.97 15.85 1.70 18.59 36.08 4.43 10.14 1.00 118.0 Moderate 2015
4 Ahmedabad 2015-01-05 48.57 95.68 22.10 21.42 37.76 15.85 22.10 39.33 39.31 7.01 18.89 2.78 118.0 Moderate 2015
In [26]:
import seaborn as sns
import matplotlib.pyplot as plt

# 假设df_city_day是您的数据集
# df_city_day = ...

# 设置绘图尺寸
sns.set(rc={'figure.figsize':(25, 15)})

# 绘制折线图
sns.lineplot(data=df_city_day, x='Year', y='AQI', hue='City', err_style=None, linewidth=3)

# 设置标题和坐标轴标签
plt.title('Trend in Average AQI on year basis', fontsize=30)
plt.xlabel('Year', fontsize=20)
plt.ylabel('Avg. AQI', fontsize=20)

# 将图例放置在图表外部的右侧
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0,prop={'size': 20})

# 显示图表
plt.show()
  • 2017 年至 2018 年，艾哈迈达巴德市的 AQI 快速上升，这与其他城市的趋势不同。
  • 造成艾哈迈达巴德如此高 AQI 的最主要污染物是高二氧化碳、二氧化硫、二氧化氮、甲苯和二甲苯。

在新年或排灯节等特定日子查看 AQI 水平将会很有趣，这些日子污染增加的可能性最高。

现在为了控制这些污染物，我们需要了解这些污染物的来源。

为了简单起见，我们将考虑两类污染物 -

  1. 工业污染物（SO2、O3、苯、甲苯、二甲苯）
  2. 车辆污染物(PM2.5,PM10,NO,NO2,NOx,NH3,CO)
In [27]:
#创建新的派生特征
Pollutants_df = df_city_day.copy()
Pollutants_df['Industrial_pollutants'] = Pollutants_df['SO2']+Pollutants_df['O3']+Pollutants_df['Benzene']+Pollutants_df['Toluene']+Pollutants_df['Xylene']
Pollutants_df['Vehical_pollutants'] = Pollutants_df['PM2.5']+Pollutants_df['PM10']+Pollutants_df['NO']+Pollutants_df['NO2']+Pollutants_df['NOx']+Pollutants_df['NH3']+Pollutants_df['CO']
Pollutants_df.drop(columns=['PM2.5','PM10','NO','NO2','NOx','NH3','CO','SO2','O3','Benzene','Toluene','Xylene'],axis=1,inplace=True)
In [28]:
Pollutants_df.head()
Out[28]:
City Date AQI AQI_Bucket Year Industrial_pollutants Vehical_pollutants
0 Ahmedabad 2015-01-01 118.0 Moderate 2015 161.02 197.31
1 Ahmedabad 2015-01-02 118.0 Moderate 2015 71.56 194.19
2 Ahmedabad 2015-01-03 118.0 Moderate 2015 85.22 243.90
3 Ahmedabad 2015-01-04 118.0 Moderate 2015 70.24 199.95
4 Ahmedabad 2015-01-05 118.0 Moderate 2015 107.32 263.48
In [29]:
#创建新的数据框
Pollutants_df_new = Pollutants_df.groupby('City')[['Industrial_pollutants','Vehical_pollutants']].mean()
In [30]:
Pollutants_df_new.head()
Out[30]:
Industrial_pollutants Vehical_pollutants
City
Ahmedabad 106.505620 299.591354
Aizawl 32.157257 86.678850
Amaravati 55.136877 171.539064
Amritsar 44.451360 253.407789
Bengaluru 46.840806 202.521229
In [31]:
#创建情节
Pollutants_df_new.plot.bar()
plt.title('Comparison for Industrial and Vehicle pollutants based on City',fontsize=30)
plt.xlabel('City',fontsize=20)
plt.ylabel('Pollutant Level',fontsize=20)
plt.show()

  • 可以看出，机动车污染物是所有城市污染的主要贡献者。

  • 德里在车辆污染方面名列前茅，而艾哈迈达巴德在工业污染方面名列前茅。

由于我们的数据集包含锁定期间的数据，因此可视化每个城市锁定之前和之后的污染图将非常有趣。

In [32]:
Pollutants_df.head()
Out[32]:
City Date AQI AQI_Bucket Year Industrial_pollutants Vehical_pollutants
0 Ahmedabad 2015-01-01 118.0 Moderate 2015 161.02 197.31
1 Ahmedabad 2015-01-02 118.0 Moderate 2015 71.56 194.19
2 Ahmedabad 2015-01-03 118.0 Moderate 2015 85.22 243.90
3 Ahmedabad 2015-01-04 118.0 Moderate 2015 70.24 199.95
4 Ahmedabad 2015-01-05 118.0 Moderate 2015 107.32 263.48
In [33]:
#数据范围
print(Pollutants_df['Date'].min())
print(Pollutants_df['Date'].max())

2015-01-01 00:00:00
2020-07-01 00:00:00
In [34]:
#创建 precovid 数据集
precovid_df = Pollutants_df[Pollutants_df['Date']<='2020-03-24']
In [35]:
precovid_df.shape
Out[35]:
(26957, 7)
In [36]:
#创建 postcovid 数据集
postcovid_df = Pollutants_df[Pollutants_df['Date']>'2020-03-24']
postcovid_df.shape
Out[36]:
(2574, 7)
In [37]:
# 过滤数据以创建新数据
precovid_df_new = precovid_df.groupby('City')[['Industrial_pollutants','Vehical_pollutants']].mean()
In [38]:
#创建情节
precovid_df_new.plot.bar()
plt.title('Pre-Covid Comparison for Industrial and Vehicle pollutants based on City',fontsize=30)
plt.xlabel('City',fontsize=20)
plt.ylabel('Pollutant Level',fontsize=20)
plt.show()
In [39]:
# 过滤数据以创建新数据
postcovid_df_new = postcovid_df.groupby('City')[['Industrial_pollutants','Vehical_pollutants']].mean()
In [40]:
#创建情节
postcovid_df_new.plot.bar()
plt.title('Post-Covid Comparison for Industrial and Vehicle pollutants based on City',fontsize=30)
plt.xlabel('City',fontsize=20)
plt.ylabel('Pollutant Level',fontsize=20)
plt.show()

有趣的是，看到后疫情场景如何因封锁而将车辆污染物降低到显着水平。

  • 下降最显着的是德里。
  • Brajrajnagar 表现出异常，新冠疫情后车辆污染物水平有所上升。

相关性¶

In [41]:
# 变量之间的相关性
plt.figure(figsize=(15, 8))
# 显式设置 numeric_only 参数为 True
correlation = df_city_day.corr(numeric_only=True)
sns.heatmap(correlation, annot=True, cmap='coolwarm')
Out[41]:
<Axes: >

没有观察到显着相关性

数据准备¶

In [42]:
#创建数据副本
model_data = df_city_day.copy()
model_data.head()
Out[42]:
City Date PM2.5 PM10 NO NO2 NOx NH3 CO SO2 O3 Benzene Toluene Xylene AQI AQI_Bucket Year
0 Ahmedabad 2015-01-01 48.57 95.68 0.92 18.22 17.15 15.85 0.92 27.64 133.36 0.00 0.02 0.00 118.0 Moderate 2015
1 Ahmedabad 2015-01-02 48.57 95.68 0.97 15.69 16.46 15.85 0.97 24.55 34.06 3.68 5.50 3.77 118.0 Moderate 2015
2 Ahmedabad 2015-01-03 48.57 95.68 17.40 19.30 29.70 15.85 17.40 29.07 30.70 6.80 16.40 2.25 118.0 Moderate 2015
3 Ahmedabad 2015-01-04 48.57 95.68 1.70 18.48 17.97 15.85 1.70 18.59 36.08 4.43 10.14 1.00 118.0 Moderate 2015
4 Ahmedabad 2015-01-05 48.57 95.68 22.10 21.42 37.76 15.85 22.10 39.33 39.31 7.01 18.89 2.78 118.0 Moderate 2015

由于 AQI_Bucket 是 AQI 的派生特征，并且 AQI 是我们的预测目标变量，因此我们将删除 AQI_Bucket 列。

Year 是我们从 Date 派生的一个功能，因此我们还将删除 Date 列。

In [43]:
# 删除日期列
model_data.drop(columns=['Date','AQI_Bucket'],axis=1,inplace=True)
In [44]:
#名义变量的一种热编码
model_data = pd.get_dummies(model_data, columns=["City","Year"])
In [45]:
#标准化连续特征
from scipy.stats import zscore
model_data[['PM2.5','PM10','NO','NO2','NOx','NH3','CO','SO2','O3','Benzene','Toluene','Xylene']]=model_data[['PM2.5','PM10','NO','NO2','NOx','NH3','CO','SO2','O3','Benzene','Toluene','Xylene']].apply(zscore)
In [46]:
model_data.head()
Out[46]:
PM2.5 PM10 NO NO2 NOx NH3 CO SO2 O3 Benzene ... City_Shillong City_Talcher City_Thiruvananthapuram City_Visakhapatnam Year_2015 Year_2016 Year_2017 Year_2018 Year_2019 Year_2020
0 -0.26654 -0.193291 -0.731090 -0.412430 -0.472010 -0.23605 -0.183488 0.812044 4.918625 -0.200657 ... 0 0 0 0 1 0 0 0 0 0
1 -0.26654 -0.193291 -0.728765 -0.522191 -0.495418 -0.23605 -0.176052 0.630337 0.003261 0.057542 ... 0 0 0 0 1 0 0 0 0 0
2 -0.26654 -0.193291 0.035219 -0.365576 -0.046258 -0.23605 2.267235 0.896135 -0.163059 0.276450 ... 0 0 0 0 1 0 0 0 0 0
3 -0.26654 -0.193291 -0.694820 -0.401151 -0.444192 -0.23605 -0.067495 0.279859 0.103251 0.110164 ... 0 0 0 0 1 0 0 0 0 0
4 -0.26654 -0.193291 0.253765 -0.273602 0.227173 -0.23605 2.966167 1.499475 0.263137 0.291184 ... 0 0 0 0 1 0 0 0 0 0

5 rows × 45 columns

In [47]:
#检查列
model_data.columns
Out[47]:
Index(['PM2.5', 'PM10', 'NO', 'NO2', 'NOx', 'NH3', 'CO', 'SO2', 'O3',
       'Benzene', 'Toluene', 'Xylene', 'AQI', 'City_Ahmedabad', 'City_Aizawl',
       'City_Amaravati', 'City_Amritsar', 'City_Bengaluru', 'City_Bhopal',
       'City_Brajrajnagar', 'City_Chandigarh', 'City_Chennai',
       'City_Coimbatore', 'City_Delhi', 'City_Ernakulam', 'City_Gurugram',
       'City_Guwahati', 'City_Hyderabad', 'City_Jaipur', 'City_Jorapokhar',
       'City_Kochi', 'City_Kolkata', 'City_Lucknow', 'City_Mumbai',
       'City_Patna', 'City_Shillong', 'City_Talcher',
       'City_Thiruvananthapuram', 'City_Visakhapatnam', 'Year_2015',
       'Year_2016', 'Year_2017', 'Year_2018', 'Year_2019', 'Year_2020'],
      dtype='object')
In [48]:
#目标变量“AQI”的分布
import seaborn as sns
plt.figure(figsize=(7,7))
sns.histplot(model_data['AQI'],color="grey")
Out[48]:
<Axes: xlabel='AQI', ylabel='Count'>
In [49]:
# 因变量的偏度
model_data['AQI'].skew()
Out[49]:
3.76975148929373

我们的目标变量是正偏态的，我们将尝试使用变换对其进行标准化。

In [50]:
# 使用不同的变换检查偏度
print(np.log10(model_data['AQI']).skew())
print(np.sqrt(model_data['AQI']).skew())
print(1/(model_data['AQI']).skew())

0.4742432997234769
1.6788266180736655
0.2652694754123837

由于倒数变换的偏度值在三种变换中最低。我们认为这种转变是最好的。

In [51]:
# 将因变量重新定位到最后一个索引
last_column = model_data.pop('AQI')
model_data.insert(44, 'AQI', last_column)
model_data.head()
Out[51]:
PM2.5 PM10 NO NO2 NOx NH3 CO SO2 O3 Benzene ... City_Talcher City_Thiruvananthapuram City_Visakhapatnam Year_2015 Year_2016 Year_2017 Year_2018 Year_2019 Year_2020 AQI
0 -0.26654 -0.193291 -0.731090 -0.412430 -0.472010 -0.23605 -0.183488 0.812044 4.918625 -0.200657 ... 0 0 0 1 0 0 0 0 0 118.0
1 -0.26654 -0.193291 -0.728765 -0.522191 -0.495418 -0.23605 -0.176052 0.630337 0.003261 0.057542 ... 0 0 0 1 0 0 0 0 0 118.0
2 -0.26654 -0.193291 0.035219 -0.365576 -0.046258 -0.23605 2.267235 0.896135 -0.163059 0.276450 ... 0 0 0 1 0 0 0 0 0 118.0
3 -0.26654 -0.193291 -0.694820 -0.401151 -0.444192 -0.23605 -0.067495 0.279859 0.103251 0.110164 ... 0 0 0 1 0 0 0 0 0 118.0
4 -0.26654 -0.193291 0.253765 -0.273602 0.227173 -0.23605 2.966167 1.499475 0.263137 0.291184 ... 0 0 0 1 0 0 0 0 0 118.0

5 rows × 45 columns

因变量和自变量分割

In [52]:
#分割数据
X= model_data.iloc[:,:-1]
y= 1/model_data['AQI']

我们的数据已准备好用于模型

模型选择¶

这里我们的任务是根据可用的污染物数据来预测 AQI。该任务可以使用回归模型来完成。 现在的问题是我们应该使用哪种回归模型？

这个问题的答案取决于业务需求。因为我们知道，当我们使用复杂模型提高预测的准确性时，模型的可解释性就会降低，这可能会给利益相关者带来问题。

考虑到上述情况，我们可以有三种选择——

  • 模型简单，精度低，但可解释性高
  • 模型复杂，准确率高，但可解释性低
  • 上述两种情况的折衷模型

Accuracy_vs_interpretability.png

上图最能说明模型的准确性和可解释性之间的关系

现在让我们尝试实现回归模型并观察结果。

＃ 预测函数

In [53]:
#创建预测函数来测试不同模型
from sklearn.model_selection import train_test_split 
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import r2_score
from sklearn.metrics import mean_squared_error
def predict(ml_model,X,y):
    X_train,X_test,y_train,y_test=train_test_split(X,y,test_size=0.25,random_state=12)
    sc=StandardScaler()
    X_train=sc.fit_transform(X_train)
    X_test=sc.transform(X_test)
    model=ml_model.fit(X_train,y_train)
    y_pred=model.predict(X_test)
    plt.scatter(y_pred,y_test,color='b')
    plt.xlabel('Predicted')
    plt.ylabel('Actual')
# 评估指标
    evaluation_metric=(f'R^2 is {r2_score(10**(y_test), 10**(y_pred))}\n Adj R^2 is {1-(1-r2_score(10**(y_test), 10**(y_pred)))*((X_test.shape[0]-1)/(X_test.shape[0]-X_test.shape[1]-1))}\n RMSE is: {np.sqrt(mean_squared_error(10**(y_test), 10**(y_pred)))}')
    print(evaluation_metric)
In [54]:
sns.set(rc={'figure.figsize':(10,7)})

1 线性回归¶

In [55]:
from sklearn.linear_model import LinearRegression
predict(LinearRegression(),X,y)

R^2 is 0.5727955326783882
 Adj R^2 is 0.5702339359814474
 RMSE is: 0.00890246727125143

2 多项式回归¶

In [56]:
from sklearn.preprocessing import PolynomialFeatures  
poly_regs= PolynomialFeatures(degree= 2)  
x_poly= poly_regs.fit_transform(X)  
predict(LinearRegression(),x_poly,y)

R^2 is 0.743583398924464
 Adj R^2 is 0.7017697575012436
 RMSE is: 0.0068970847496440295

3 决策树回归器¶

In [57]:
from sklearn.tree import DecisionTreeRegressor
predict(DecisionTreeRegressor(),X,y)

R^2 is 0.6937112748902108
 Adj R^2 is 0.6918747112618611
 RMSE is: 0.007538032369465511

4 随机森林回归器¶

In [58]:
from sklearn.ensemble import RandomForestRegressor
predict(RandomForestRegressor(),X,y)

R^2 is 0.8385357477899551
 Adj R^2 is 0.8375675783844982
 RMSE is: 0.005473066005065331

5 XGBoost 回归器¶

In [59]:
import xgboost as xgb
predict(xgb.XGBRegressor(),X,y)

R^2 is 0.8300554466563952
 Adj R^2 is 0.8290364278028767
 RMSE is: 0.005614952905727235

结论¶

我们可以清楚地观察到，随机森林回归器最适合我们的数据，具有最佳 R2 分数和最小均方根误差 (RMSE)