AQI全称是Air Quality Index，指空气质量指数，用来衡量空气清洁或者污染的程度，值越小，表示空气质量越好。
本文的分析目标是：
一、描述性统计
那些城市的空气质量较好/较差？ 空气质量在地理位置分布上，是否具有一定的规律？ 二、推断统计
三、相关系数分析
四、区间估计
五、统计建模
导包并读取数据：
import numpy as np
import pandas as pd
import matplotlib. pyplot as plt
import seaborn as sns
sns. set ( style= "darkgrid" )
plt. rcParams[ "font.family" ] = "SimHei"
plt. rcParams[ "axes.unicode_minus" ] = False
data = pd. read_csv( "data/data.csv" )
print ( data. shape)
data. head( )
数据集描述：
检查缺失值：
data. isnull( ) . sum ( axis= 0 )
查看含缺失值列数据的分布：
#print(data["Precipitation"].skew())#偏度
sns. distplot( data[ "Precipitation" ] . dropna( ) ) #要删除NA值才能做分布密度图
plt. title( "分布密度图" )
0.27360760671177387
数值型变量，数据呈现右偏分布，所以使用中位数填充。
对缺失值进行中位数填充
data. fillna( { "Precipitation" : data[ "Precipitation" ] . median( ) } , inplace= True )
data.describe() 查看数据的描述：分位数、均值与标准差 基于正太分布 ±三个标准差涵盖99.7%的数据 箱线图（四分位距IQR=Q3-Q1，上下边界：Q3/Q1 ±1.5IQR） 查看数据集的偏度：
data. skew( )
AQI 1.198754
Precipitation 0.273608
GDP 3.761428
Temperature -0.597343
Longitude -1.407505
Latitude 0.253563
Altitude 3.067242
PopulationDensity 3.125853
GreenCoverageRate -0.381786
Incineration(10,000ton) 4.342614
dtype: float64
可以看到GDP和人口密度等都出现了严重的右偏分布，意味着存在很多极大的异常值。
下面我们查看以下GDP的异常值：
mean, std = data. GDP. mean( ) , data. GDP. std( )
lower, upper = mean - 3 * std, mean + 3 * std
print ( "均值：" , mean)
print ( "标准差：" , std)
print ( "下限：" , lower)
print ( "上限：" , upper)
data. loc[ ( data. GDP < lower) | ( data. GDP > upper) , "GDP" ]
均值： 2390.901815384616
标准差： 3254.876921271434
下限： -7373.728948429687
上限： 12155.532579198918
16 22968.60
63 18100.41
202 24964.99
207 17502.99
215 14504.07
230 16538.19
256 17900.00
314 15719.72
Name: GDP, dtype: float64
可以通过以下代码查看每个数值列异常值的个数：
def func ( s) :
mean, std = s. mean( ) , s. std( )
lower, upper = mean - 3 * std, mean + 3 * std
a, b = ( s < lower) . sum ( ) , ( s > upper) . sum ( )
return pd. Series( ( a, b, a+ b) , index= ( "极小异常值" , "极大异常值" , "总和" ) )
data. select_dtypes( include= 'number' ) . agg( func) . T
箱线图
我们可以直接通过箱线图 来查看数据集整体的异常情况：
plt. figure( figsize= ( 16 , 7 ) )
plt. xticks( rotation= 45 )
sns. boxplot( data= data)
也可以单独查看某个列的异常情况：
sns. boxplot( data= data[ "GDP" ] )
处理方式有：
删除异常数据 视为缺失值处理 对数转换 使用边界值替换 分箱离散化 由于右偏分布在取对数后，往往呈现正态分布，所以我们可以对右偏分布进行取对数转换。
例如，GDP变量呈现严重的右偏分布，我们可以进行取对数转换：
fig, ax = plt. subplots( 1 , 2 )
fig. set_size_inches( 15 , 5 )
sns. distplot( data[ "GDP" ] , ax= ax[ 0 ] )
sns. distplot( np. log( data[ "GDP" ] ) , ax= ax[ 1 ] )
查看处理前的数据：
data. select_dtypes( "number" ) . agg( [ "max" , "min" ] )
用边界值替换处理：
data_c = data. copy( )
for column in data. columns:
s = data_c[ column]
if pd. api. types. is_numeric_dtype( s) :
quantile = np. quantile( s, [ 0.25 , 0.75 ] )
IQR = quantile[ 1 ] - quantile[ 0 ]
lower = quantile[ 0 ] - 1.5 * IQR
upper = quantile[ 1 ] + 1.5 * IQR
s. clip( lower, upper, inplace= True )
data_c. describe( )
再查看：
data_c. describe( )
结果：
再次查看箱线图：
plt. figure( figsize= ( 16 , 7 ) )
plt. xticks( rotation= 45 )
sns. boxplot( data= data_c)
# 发现重复值。
print ( data. duplicated( ) . sum ( ) )
# 查看哪些记录出现了重复值。
data[ data. duplicated( keep= False ) ]
结果：
2
重复值对数据分析通常没有什么用，直接删除即可：
data. drop_duplicates( inplace= True )
print ( data. duplicated( ) . sum ( ) ) #去重后检查
结果：
0
最好的5个城市
t = data[ [ "City" , "AQI" ] ]
t. nsmallest( 5 , "AQI" )
可以发现：空气质量最好的5个城市为 韶关、南平、梅州、基隆、三明
最差的8个城市
t. nlargest( 5 , "AQI" )
可以发现：空气质量最差的5个城市为 北京、朝阳、保定、锦州、东营
对于AQI，可以对空气质量进行等级划分，划分标准为：
分级如函数：
import math
levels = [ "一级" , "二级" , "三级" , "四级" , "五级" , "六级" ]
def value_to_level ( AQI) :
i = math. ceil( AQI/ 50 ) - 1
if i == 5 :
i = 4
elif i > 5 :
i = 5
return levels[ i]
level = data[ "AQI" ] . apply ( value_to_level)
display( level. value_counts( ) )
sns. countplot( x= level, order= levels)
结果：
sns. scatterplot( x= "Longitude" , y= "Latitude" , hue= "AQI" , palette= plt. cm. RdYlGn_r, data= data)
plt. title( "空气质量指数分布图" , size= 16 )
plt. show( )
结果：
从上图大致的地理位置来看，西部城市普遍好于东部城市，南部城市普遍好于北部城市。
首先看下临海城市与内陆城市的数量：
data[ "Coastal" ] . value_counts( )
否 243
是 80
Name: Coastal, dtype: int64
分类散点图：
sns. stripplot( x= "Coastal" , y= "AQI" , data= data)
分簇散点图：
sns. swarmplot( 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= "r" , data= data)
均值
display( data. groupby( "Coastal" ) [ "AQI" ] . mean( ) )
sns. barplot( x= "Coastal" , y= "AQI" , data= data) #只能看均值数据
结果：
Coastal
否 79.045267
是 64.062500
Name: AQI, dtype: float64
注：柱形图上方的竖线表示置信区间。
还可以查看箱线图：
sns. boxplot( x= "Coastal" , y= "AQI" , data= data)
而小提琴图除了展示箱线图的信息外，还能呈现分布的密度。
就样本数据的分布图而言，我们可以看到沿海城市的AQI分布普遍低于内陆城市。但我们仅仅只有样本数据，对于总体的支持力度还需假设检验。
假设检验： 检验内陆AQI和沿海城市AQI的均值是否相等（两样本t检验）
from scipy import stats
coastal = data. loc[ data[ "Coastal" ] == "是" , "AQI" ]
inland = data. loc[ data[ "Coastal" ] == "否" , "AQI" ]
首先进行方差齐性检验（检验两样本的方差是否一致）：
stats. levene( coastal, inland)
结果：
LeveneResult(statistic=0.08825036641952543, pvalue=0.7666054880248168)
方差齐性检验，原假设H0为方差相等（齐性），显然p-value超过显著性水平，说明两样本的方差是相等的。
下面进行两样本t检验：
#equal_var=True表示两样本方差一致
r = stats. ttest_ind( coastal, inland, equal_var= True )
r
结果：
Ttest_indResult(statistic=-2.7303827520948905, pvalue=0.006675422541012958)
P值小于显著性水平，所以拒绝原假设，接受备择假设，即内陆AQI和沿海城市AQI的均值不相等。
由于统计量由左样本减右样本得到，统计量小于0说明，临海的均值小于内陆。
下面我们假设临海空气质量的均值小于内陆空气质量的均值，则这是一个右边假设检验，可以通过以下方法得到P值：
p = stats. t. sf( r. statistic, df= len ( coastal) + len ( inland) - 2 )
p
结果：
0.9966622887294936
说明我们有99.7%的信心认为临海的空气质量整体好于内陆空气质量（均值小于内陆）。
一般我们会有如下疑问：
人口密度大，是否会对空气质量造成负面影响 绿化率高，是否能提高空气质量 首先查看散点图矩阵：
sns. pairplot( data[ [ "AQI" , "PopulationDensity" , "GreenCoverageRate" ] ] )
#参数 kind="reg"给散点图绘制一条回归线
# sns.pairplot(data[["AQI", "PopulationDensity", "GreenCoverageRate"]], kind="reg")
相关系数：
变量X与变量Y的协方差为：
Cov
(
X
,
Y
)
=
∑
i
=
1
n
(
x
i
−
x
ˉ
)
(
y
i
−
y
ˉ
)
n
−
1
\large{\operatorname{Cov}(X, Y)=\frac{\sum_{i=1}^{n}\left(x_{i}-\bar{x}\right)\left(y_{i}-\bar{y}\right)}{n-1}}
C o v ( X , Y ) = n − 1 ∑ i = 1 n ( x i − x ˉ ) ( y i − y ˉ )
r
(
X
,
Y
)
=
Cov
(
X
,
Y
)
σ
X
σ
Y
\large{r(X, Y)=\frac{\operatorname{Cov}(X, Y)}{\sigma_{X} \sigma_{Y}}}
r ( X , Y ) = σ X σ Y C o v ( X , Y )
=
∑
i
=
1
n
(
x
i
−
x
ˉ
)
(
y
i
−
y
ˉ
)
∑
i
=
1
n
(
x
i
−
x
ˉ
)
2
∑
i
=
1
n
(
y
i
−
y
ˉ
)
2
\large{=\frac{\sum_{i=1}^{n}\left(x_{i}-\bar{x}\right)\left(y_{i}-\bar{y}\right)}{\sqrt{\sum_{i=1}^{n}\left(x_{i}-\bar{x}\right)^{2} \sum_{i=1}^{n}\left(y_{i}-\bar{y}\right)^{2}}}}
= ∑ i = 1 n ( x i − x ˉ ) 2 ∑ i = 1 n ( y i − y ˉ ) 2
∑ i = 1 n ( x i − x ˉ ) ( y i − y ˉ )
相关系数的绝对值的取值范围是[0,1]，可以根据相关系数衡量两个变量的相关性，正数代表正相关，负数代码负相关。绝对值的大小体现相关的程度：
0.8-1：极强相关
0.6-0.8：强相关
0.4-0.6：中等程度相关
0.2-0.4：弱相关
0-0.2：极弱相关或无相关
x = data[ "AQI" ]
y = data[ "Precipitation" ]
# 计算AQI与Precipitation的协方差。
a = ( x - x. mean( ) ) * ( y - y. mean( ) )
cov = np. sum ( a) / ( len ( a) - 1 )
print ( "协方差：" , cov)
# 计算AQI与Precipitation的相关系数。
corr = cov / np. sqrt( x. var( ) * y. var( ) )
print ( "相关系数：" , corr)
结果：
协方差： -10098.209013903044
相关系数： -0.40184407003013883
其实可以直接使用内置方法计算：
print ( "协方差：" , x. cov( y) )
print ( "相关系数：" , x. corr( y) )
结果：
协方差： -10098.209013903042
相关系数： -0.4018440700301393
调用DataFrame中的方法计算相关系数：
data. corr( )
使用热力图呈现相关性数值更佳：
plt. figure( figsize= ( 15 , 10 ) )
ax = sns. heatmap( data. corr( ) , cmap= plt. cm. RdYlGn, annot= True , fmt= ".2f" )
结果：
从上述相关性矩阵我们可以看到，AQI与人口密度绿化面积几乎不相关。
AQI与纬度（0.55）和降雨量（-0.4）的相关性最高。说明：
纬度越低，AQI越低，空气质量越好。 降雨量越多，AQI越低，空气质量越好。 传言说，全国所有城市的空气质量指数均值为71左右，这个消息可靠吗？
下面作为原假设全部城市的AQI均值为71，并进行假设检验：
print ( "样本均值 ：" , data[ "AQI" ] . mean( ) )
print ( stats. ttest_1samp( data[ "AQI" ] , 71 ) )
结果：
样本均值 ： 75.3343653250774
Ttest_1sampResult(statistic=1.8117630617496872, pvalue=0.07095431526986647)
可以看到，P值大于显著性水平0.05，故我们没有充足的证据拒绝原假设，于是接受原假设全部城市的AQI均值为71。
下面计算一下置信区间：
stats. t. interval( 0.95 , df= len ( data) - 1 , loc= data. AQI. mean( ) , scale= stats. sem( data. AQI) )
结果： (70.6277615675309, 80.0409690826239)
因此，我们可以认为全国所有城市的平均空气质量，95%的可能在(70.63, 80.04)范围内。
问题：已知某市的降雨量、温度、经纬度等指标，如何预测其空气质量？
为了进行模型计算，我们首先需要将一些文本型变量转换为离散数值变量。
data[ "Coastal" ] = data[ "Coastal" ] . map ( { "是" : 1 , "否" : 0 } )
data[ "Coastal" ] . value_counts( )
0 243
1 80
Name: Coastal, dtype: int64
注意：只有两个类别时，转换为任意两个数值都可以。如果存在三个以上的分类，可以使用独热编码对变量进行转换。
基础模型：
#线性回归模型
from sklearn. linear_model import LinearRegression
#数据集的处理。训练集和测试集
from sklearn. model_selection import train_test_split
#构建X、y变量
X = data. drop( [ "City" , "AQI" ] , axis= 1 ) #删除多余变量 ，Y变量不能在此
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)
#R方值 模型评价
print ( lr. score( X_train, y_train) )
print ( lr. score( X_test, y_test) )
结果：
0.4538897765064037
0.4040770562383651
绘制预测结果图
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 )
plt. show( )
下面我们将尝试优化模型提高模型效果。
首先我们使用临界值替换异常值：
for column in X. columns. drop( "Coastal" ) :
s = X_train[ column]
if pd. api. types. is_numeric_dtype( s) :
quantile = np. quantile( s, [ 0.25 , 0.75 ] )
IQR = quantile[ 1 ] - quantile[ 0 ]
lower = quantile[ 0 ] - 1.5 * IQR
upper = quantile[ 1 ] + 1.5 * IQR
s. clip( lower, upper, inplace= True )
X_test[ column] . clip( lower, upper, inplace= True )
再次查看模型效果：
lr. fit( X_train, y_train)
print ( lr. score( X_train, y_train) )
print ( lr. score( X_test, y_test) )
结果：
0.4631142291492417
0.44614202658395546
可以看到模型相对之前有一定的提高。
下面使用RFE(Recursive feature elimination 递归特征消除) 方法实现特征选择 ，因为删除一些不重要的特征，反而有助于模型效果的提高。
from sklearn. feature_selection import RFECV
# estimator： 要操作的模型。
# step： 每次删除的变量数。
# cv： 使用的交叉验证折数。
# n_jobs： 并发的数量。
# scoring: 评估的方式。
rfecv = RFECV( estimator= lr, step= 1 , cv= 5 , n_jobs= - 1 , scoring= "r2" )
rfecv. fit( X_train, y_train)
print ( "剩余的特征数量:" , rfecv. n_features_)
# 返回经过特征选择后，使用缩减特征训练后的模型。
print ( rfecv. estimator_)
# 返回每个特征的等级，数值越小，特征越重要。
print ( rfecv. ranking_)
print ( "被选择的特征布尔数组：" , rfecv. support_)
print ( "交叉验证的评分:" , rfecv. grid_scores_)
结果：
9
LinearRegression(copy_X=True, fit_intercept=True, n_jobs=1, normalize=False)
[1 1 1 1 1 1 2 1 1 1]
[ True True True True True True False True True True]
[-0.06091362 0.1397744 0.19933237 0.16183209 0.18281661 0.20636585
0.29772708 0.307344 0.30877162 0.30022701]
绘图表示，在特征选择的过程中，使用交叉验证获取R平方的值
plt. plot( range ( 1 , len ( rfecv. grid_scores_) + 1 ) , rfecv. grid_scores_, 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)
X_train_eli = X_train[ X_train. columns[ rfecv. support_] ]
X_test_eli = X_test[ X_test. columns[ rfecv. support_] ]
print ( rfecv. estimator_. score( X_train_eli, y_train) )
print ( rfecv. estimator_. score( X_test_eli, y_test) )
结果：
剔除的变量： Index(['PopulationDensity'], dtype='object')
0.46306656191488593
0.44502255894081927
下面我们对部分数据使用分箱操作，并进行one_Hot编码：
from sklearn. preprocessing import KBinsDiscretizer
# KBinsDiscretizer K个分箱的离散器。用于将数值（通常是连续变量）变量进行区间离散化操作。
# n_bins：分箱（区间）的个数。
# 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 , 14 , 6 ] ,
encode= "onehot-dense" , strategy= "uniform" )
# 定义离散化的特征。
discretize = [ "Longitude" , "Temperature" , "Precipitation" , "Latitude" ]
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( ) )
再重新查看模型效果：
lr. fit( X_train_dis, y_train)
print ( lr. score( X_train_dis, y_train) )
print ( lr. score( X_test_dis, y_test) )
结果：
0.6892388692774563
0.6546062348355671
可以看到模型在经过分箱离散化操作后，预测效果大幅度提高。
残差就是模型的预测值与真实值之间的差异，可以绘制残差图对模型进行评估，横坐标为预测值，纵坐标为真实值
对于一个模型，误差应该的随机性的，而不是有规律的。残差也随机分布于中心线附近，如果残差图中残差是有规律的，则说明模型遗漏了某些能够影响残差的解释信息。
异方差性 ，是指残差具有明显的方差不一致性。
绘制残差图：
fig, ax = plt. subplots( 1 , 2 )
fig. set_size_inches( 15 , 5 )
data = [ X_train, X_train_dis]
title = [ "原始数据" , "处理后数据" ]
for d, a, t in zip ( data, ax, title) :
model = LinearRegression( )
model. fit( d, y_train)
y_hat_train = model. predict( d)
residual = y_hat_train - y_train. values
a. set_xlabel( "预测值" )
a. set_ylabel( " 残差" )
a. axhline( y= 0 , color= "red" )
a. set_title( t)
sns. scatterplot( x= y_hat_train, y= residual, ax= a)
在左图中，可以看出随着预测值的增大，残差也有增大的趋势。此时我们可以使用对y值取对数的方式处理。
model = LinearRegression( )
y_train_log = np. log( y_train)
y_test_log = np. log( y_test)
model. fit( X_train, y_train_log)
y_hat_train = model. predict( X_train)
residual = y_hat_train - y_train_log. values
plt. xlabel( "预测值" )
plt. ylabel( " 残差" )
plt. axhline( y= 0 , color= "red" )
sns. scatterplot( x= y_hat_train, y= residual)
此时异方差性得到解决。模型效果也可能会得到进一步提升。
检查离群点：对于多元线性回归，其回归线已经成为超平面，无法通过可视化来观测。但我们可以通过残差图，通过预测值与实际值之间的关系，来检测离群点。我们认为偏离2倍标准差的点为离群点：
y_hat_train = lr. predict( X_train_dis)
residual = y_hat_train - y_train. values
r = ( residual - residual. mean( ) ) / residual. std( )
plt. xlabel( "预测值" )
plt. ylabel( " 残差" )
plt. axhline( y= 0 , color= "red" )
sns. scatterplot( x= y_hat_train[ np. abs ( r) <= 2 ] ,
y= residual[ np. abs ( r) <= 2 ] , color= "b" , label= "正常值" )
sns. scatterplot( x= y_hat_train[ np. abs ( r) > 2 ] , y= residual[ np. abs ( r) > 2 ] , color= "orange" , label= "异常值" )
结果：
剔除异常值后，训练模型并查看效果：
X_train_dis_filter = X_train_dis[ np. abs ( r) <= 2 ]
y_train_filter = y_train[ np. abs ( r) <= 2 ]
lr. fit( X_train_dis_filter, y_train_filter)
print ( lr. score( X_train_dis_filter, y_train_filter) )
print ( lr. score( X_test_dis, y_test) )
结果：
0.7354141753913532
0.6302724058812207
可以看到模型效果得到了进一步的提升。
从空气质量总体分布上来说，南方城市优于北方城市，西部城市优于东部城市。 临海城市空气质量整体好于内陆城市 城市是否临海，降雨量，以及维度对空气质量指数影响最大 有95%的把握可以说城市平均空气质量指数在区间(70.63, 80.04)之间 通过历史数据，我们可以对空气质量指数进行预测
