感谢Dr.fish一路以来的耐心讲解和细致回答，悄悄的就走到了最后一次，好不舍得大家。
本次课的随堂作业如下：
空气质量分析
直接上代码
# 导入分析包
import pandas as pd
import numpy as np
import scipy.stats
import matplotlib.pyplot as plt
import matplotlib
matplotlib.style.use('ggplot') # 设置图片采用ggplot样式
import seaborn as sns
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
数据导入及清洗
# 导入4年数据
data_2014 = pd.read_csv('homework_data/beijing2014.csv')
data_2015 = pd.read_csv('homework_data/beijing2015.csv')
data_2016 = pd.read_csv('homework_data/beijing2016.csv')
data_2017 = pd.read_csv('homework_data/beijing2017.csv')
# 基本处理
# 加个年份字段，省的每次还要resample
data_2014['year'] = '2014'
data_2015['year'] = '2015'
data_2016['year'] = '2016'
data_2017['year'] = '2017'
# 拼接数据表
data = pd.concat([data_2014,data_2015,data_2016,data_2017],axis=0)
data.columns = ['date','AQI','PM10','PM25','CO','NO2','O3','SO2','year']
data.head()
head()结果
data.info()
info()结果
数据不太整齐，瞧着像是开始只是关注了空气污染指数（AQI），并没有关注空气中的一些……杂质（作为一个文科生，刚上来就词穷……请允许小白自己去墙角冷静下），后面才开始慢慢完善监控体系的。理论上看到这样的数据应该是二话不说直接咔嚓含缺失值的行的，但是考虑到衡量“空气怎么样了”这件事，AQI还是有比较重的话语权的，实在是不舍得删掉上百条数据（啧啧，穷人家的孩子总是能把节俭发挥到任何一个领域，嗯），况且如果只是做单一数据分析的话，缺不缺其实关系也不太大。
所以，经过严肃纠结后，小白做了个严重的决定：后面做高级分析的时候再生成一份严格阉割过的数据。
# 替换日期数据类型，定义index
data['date'] = pd.to_datetime(data['date'],format='%Y%m%d')
data.index = data['date']
数据导入及清洗小结
- 一共4个数据表，每年一张，包括2014-2017年空气质量数据，其中2017年非整年数据；
- 为了便于计算，将4年数据合并为一张df，并增加了year字段，便于探索及分析；
- 数据表共包含9个字段：日期、空气污染指数（AQI）、PM10、PM2.5、一氧化碳（CO）、二氧化氮（NO2）、臭氧（O3）、二氧化硫（SO2）、年份；
- 数据存在比较多缺失值，考虑到数据的完整性，并未对存在缺失值的行进行删除，如后面分析需做缺失值删除，再行生成新df。
数据探索
计划步骤
考虑到空气污染指数（AQI）为标准化后对城市空气的判断标准，故会主要针对AQI进行探索及分析。对于其他几项污染物，计划仅进行相关性分析。
part_1. 观察AQI变化情况
- 观察各年均值
- 计算年度均值置信区间（了解年度数据分布情况）
- 计算两两年置信区间（了解年间数据变化情况）
- 假设检验佐证结论
part_2. 观察污染物相关关系
- 观察两个PM类污染物的相关关系（相互 & 与AQI）
- 观察四个与排放相关污染物的相关关系（与AQI）
part_1. 观察AQI变化情况
1. 计算均值及绘制图示
# 计算AQI年均值（2种方法）
# 方法一：用groupby year计算
data.groupby(['year']).AQI.mean()
# 方法二：用resample计算
data.resample('A').AQI.mean().to_period('A') # 'M'-月, 'A'-年, 'Q'-季度,'W'-周。'BM','BA','BQ'-business
计算年度AQI均值
# 绘制年AQI箱图
data.boxplot(['AQI'], by = 'year', figsize = (10,6))
年AQI箱图
港真，看到统计结果和箱图，小白真的要开始认真的怀疑起自己的认知了。不管是观察箱图中的均值还是箱体，貌似帝都的空气都是在变好的啊。但可是，可但是，作为帝都的原住民，明明小白觉得愈发的无法呼吸了呢。是小白太敏感了，还是被自己的有眼无珠骗了。这得拉数据瞧瞧，看是不是极端天气被算数平均数给匀成“天下大同”了。
# 每天 及 周平均 AQI
data.plot(y =[ 'AQI'],figsize=(20, 6))
plt.title('BY DADY AQI (2014 - 2017)')
data.resample('W').mean().plot(y =[ 'AQI'],figsize=(20, 6))
plt.title('BY WEEK AQI (2014 - 2017)')
每天及每周平均AQI
嗯，这就对了，像蹦极心电图一样的图谱，果然坏天气是“被平均”了。期待后面的分析结果是不是可以佐证小白的这个“预设想法”。
BTW：看到这两张图的第一眼，真是觉得心好累，我们究竟是生活在一个什么样的环境中（@_@）。
特意问了下度娘，AQI解释如下：
空气污染指数，就是根据环境空气质量标准和各项污染物对人体健康、生态、环境的影响，将常规监测的几种空气污染物浓度简化成为单一的概念性指数值形式，它将空气污染程度和空气质量状况分级表示，适合于表示城市的短期空气质量状况和变化趋势。
空气污染指数（AQI）的取值范围与空气质量等级设定如下：
0～50 I级
51～100 II级
101～200 III级
201～300 IV级
＞300 V级
2. AQI 95%置信区间（了解年度数据分布情况）
说明：
事实上数据中存在比较多的空值，如果只是做单字段计算，影响并不大，但是如果进行置信区间的估计和假设检验，还是需要数据是对齐的。为了不影响后面的单字段计算，先将数据集data复制一份，进行缺失值行删除，以便进行区间估计及假设检验。（一下删掉了253条数据，小家子气的小白表示好心疼（>_<））
# 数据清洗-删除缺失值行
df_data = data.dropna()
df_data.info()
删除缺失值行
# 数据按年份分组
df_2014 = df_data[df_data.year == '2014']
df_2015 = df_data[df_data.year == '2015']
df_2016 = df_data[df_data.year == '2016']
df_2017 = df_data[df_data.year == '2017']
# 定义计算函数
def mean_ci(data):
'''给定样本数据，计算均值95%的置信区间'''
sample_size = len(data)
std = np.std(data, ddof=1)
se = std / np.sqrt(sample_size)
point_estimate = np.mean(data)
z_score = scipy.stats.norm.isf(0.025) # 置信度95%
confidence_interval = (point_estimate - z_score * se, point_estimate + z_score * se)
return confidence_interval
# 4年AQI置信度95%置信区间
mean_ci_2014 = mean_ci(df_2014.AQI)
mean_ci_2015 = mean_ci(df_2015.AQI)
mean_ci_2016 = mean_ci(df_2016.AQI)
mean_ci_2017 = mean_ci(df_2017.AQI)
# 4年AQI均值
mean_2014 = df_2014.AQI.mean()
mean_2015 = df_2015.AQI.mean()
mean_2016 = df_2016.AQI.mean()
mean_2017 = df_2017.AQI.mean()
print 'mean_2014: ',mean_2014
print 'mean_ci_2014: ',mean_ci_2014
print 'mean_2015: ',mean_2015
print 'mean_ci_2015: ',mean_ci_2015
print 'mean_2016: ',mean_2016
print 'mean_ci_2016: ',mean_ci_2016
print 'mean_2017: ',mean_2017
print 'mean_ci_2017: ',mean_ci_2017
-----------------------------------------------------------------------------------------------------------------------------------------------
# 输出结果
mean_2014: 109.462381564
mean_ci_2014: (99.719902564963121, 119.20486056288246)
mean_2015: 114.065970584
mean_ci_2015: (105.66456901556916, 122.46737215200551)
mean_2016: 107.196151151
mean_ci_2016: (99.143231702026412, 115.24907060084212)
mean_2017: 97.1970332075
mean_ci_2017: (86.621416349031492, 107.7726500659942)
# 绘制4年均值置信区间示意图
ci_low_2014,ci_high_2014 = mean_ci(df_2014.AQI)
ci_low_2015,ci_high_2015 = mean_ci(df_2015.AQI)
ci_low_2016,ci_high_2016 = mean_ci(df_2016.AQI)
ci_low_2017,ci_high_2017 = mean_ci(df_2017.AQI)
x4, y4 = np.random.random(10000) * (ci_high_2014 - ci_low_2014) + ci_low_2014,[4 for y in range(0, 10000)]
x5, y5 = np.random.random(10000) * (ci_high_2015 - ci_low_2015) + ci_low_2015,[3 for y in range(0, 10000)]
x6, y6 = np.random.random(10000) * (ci_high_2016 - ci_low_2016) + ci_low_2016,[2 for y in range(0, 10000)]
x7, y7 = np.random.random(10000) * (ci_high_2017 - ci_low_2017) + ci_low_2017,[1 for y in range(0, 10000)]
plt.figure(figsize=(15,4))
line_1 = plt.scatter(x = x4, y = y4, alpha = 0.5, marker = '.')
line_2 = plt.scatter(x = x5, y = y5, alpha = 0.5, marker = '.')
line_3 = plt.scatter(x = x6, y = y6, alpha = 0.5, marker = '.')
line_4 = plt.scatter(x = x7, y = y7, alpha = 0.5, marker = '.')
plt.legend((line_1, line_2, line_3, line_4),('mean_ci_2014','mean_ci_2015','mean_ci_2016','mean_ci_2017'))
plt.title('Confidence Interval for the Mean (2014 - 2017)')
plt.xlim(85, 130)
plt.ylim(0, 5)
Confidence Interval for the Mean (2014 - 2017)
通过观察年度均值和置信度95%的置信区间，貌似空气除了在2015年有一次拉低外，2016/2017两年均呈现提升趋势，这个结论真是不可思议，再看看差值之间的置信区间。
BTW：画这个置信区间示意图简直要囧死了。怎么都倒腾不出来美丽的样子，最后还是用最笨最丑的散点图方式画出来的，小白简直想自己挖个坑把自己给埋了。有哪位大神知道更漂亮的线图可以支援置信区间示意图嘛？
3. 两两年之间差值的95%置信区间（了解年间数据变化情况）
# 定义bootstrap方法
def draw_bs_reps(data, func, size=1):
# 初始化空数组
bs_replicates = np.empty(size)
# 重新重复抽样，产生一些列新样本，并计算相应统计量
for i in range(size):
bs_replicates[i] = func(np.random.choice(data, size=len(data)))
return bs_replicates
# 2014 和 2015 年差值置信度95%的置信区间
#一些列重复抽样的均值
bs2014 = draw_bs_reps(df_2014.AQI, np.mean, 10000)
bs2015 = draw_bs_reps(df_2015.AQI, np.mean, 10000)
# 计算均值差，存于bs_diff数组中
bs_diff_2015 = bs2015 - bs2014
diff_mean = df_2015.AQI.mean() - df_2014.AQI.mean()
# 计算95%的置信区间
conf_int = np.percentile(bs_diff_2015, [2.5, 97.5])
# 打印结果
print '2014-2015 diff_mean: ', diff_mean
print '2014-2015 95%: ', conf_int
-----------------------------------------------------------------------------------------------------------------------------------------------
# 输出结果
2014-2015 diff_mean: 4.60358901986
2014-2015 95%: [ -8.26736615 17.55979325]
# 2015 和 2016 年差值置信度95%的置信区间
bs2015 = draw_bs_reps(df_2015.AQI, np.mean, 10000)
bs2016 = draw_bs_reps(df_2016.AQI, np.mean, 10000)
bs_diff_2016 = bs2016 - bs2015
diff_mean = df_2016.AQI.mean() - df_2015.AQI.mean()
conf_int = np.percentile(bs_diff_2016, [2.5, 97.5])
print '2015-2016 diff_mean: ', diff_mean
print '2015-2016 95%: ', conf_int
-----------------------------------------------------------------------------------------------------------------------------------------------
# 输出结果
2015-2016 diff_mean: -6.86981943235
2015-2016 95%: [-18.2783607 4.82231891]
# 2016 和 2017 年差值置信度95%的置信区间
# 2016 和 2017 年差值置信度95%的置信区间
bs2016 = draw_bs_reps(df_2016.AQI, np.mean, 10000)
bs2017 = draw_bs_reps(df_2017.AQI, np.mean, 10000)
bs_diff_2017 = bs2017 - bs2016
diff_mean = df_2017.AQI.mean() - df_2016.AQI.mean()
conf_int = np.percentile(bs_diff_2017, [2.5, 97.5])
print '2016-2017 diff_mean: ', diff_mean
print '2016-2017 95%: ', conf_int
-----------------------------------------------------------------------------------------------------------------------------------------------
# 输出结果
2016-2017 diff_mean: -9.99911794392
2016-2017 95%: [-23.21495792 3.71191372]
# 绘制两两年之间差值置信区间示意图
conf_low_2015,conf_high_2015 = np.percentile(bs_diff_2015, [2.5, 97.5])
conf_low_2016,conf_high_2016 = np.percentile(bs_diff_2016, [2.5, 97.5])
conf_low_2017,conf_high_2017 = np.percentile(bs_diff_2017, [2.5, 97.5])
x1, y1 = np.random.random(10000) * (conf_high_2015 - conf_low_2015) + conf_low_2015,[3 for y in range(0, 10000)]
x2, y2 = np.random.random(10000) * (conf_high_2016 - conf_low_2016) + conf_low_2016,[2 for y in range(0, 10000)]
x3, y3 = np.random.random(10000) * (conf_high_2017 - conf_low_2017) + conf_low_2017,[1 for y in range(0, 10000)]
plt.figure(figsize=(15,4))
line_1 = plt.scatter(x = x1, y = y1, alpha = 0.5, marker = '.')
line_2 = plt.scatter(x = x2, y = y2, alpha = 0.5, marker = '.')
line_3 = plt.scatter(x = x3, y = y3, alpha = 0.5, marker = '.')
plt.legend((line_1, line_2, line_3),('Difference in 2014-2015','Difference in 2015-2016','Difference in 2016-2017'))
plt.title('Confidence Interval for the Difference (2014 - 2017)')
plt.xlim(-25, 20)
plt.ylim(0, 5)
Confidence Interval for the Difference (2014 - 2017)
通过对两两年差均值和置信区间的观察，貌似结论和单独观察每年的均值和置信区间结论差不多。随着负值的慢慢增大，说明天气真的好了？！
最后再放个大招，直接看2014年和2017年的均值差和差的置信区间。
# 2014 和 2017 年差值置信度95%的置信区间
bs2014 = draw_bs_reps(df_2014.AQI, np.mean, 10000)
bs2017 = draw_bs_reps(df_2017.AQI, np.mean, 10000)
bs_diff = bs2017 - bs2014
diff_mean = df_2017.AQI.mean() - df_2014.AQI.mean()
conf_int = np.percentile(bs_diff, [2.5, 97.5])
print '2014-2017 diff_mean: ', diff_mean
print '2014-2017 95%: ', conf_int
-----------------------------------------------------------------------------------------------------------------------------------------------
# 输出结果
2014-2017 diff_mean: -12.2653483564
2014-2017 95%: [-26.87177894 2.25489214]
这么看的话，天气确实是好了，最后再假设检验佐证一下。
4. 假设检验判断 AQI 2014 和 2017 是否有改变（t分布）
# 两个独立样本单边检验
t_statistics, p_value = scipy.stats.ttest_ind(df_2017.AQI, df_2014.AQI)
p_value = p_value / 2
p_value
-----------------------------------------------------------------------------------------------------------------------------------------------
# 输出结果
0.04826124865534643
原假设为“2014年空气质量均值≥2017年空气质量均值”，取α=0.05, 因为 pvalue<α, 所以原假设不成立，得到结论是：2017年比2014年空气质量有所提升。
PS. 其实如果取 α=0.01 的话，就 P-value > α了，就没有充分的理由证明2017年比2014年空气质量有所提升了，所以这天儿吧……其实还是不怎么样！
part_1 结论：
通过从2014年至2017年，对空气污染指数（AQI）从均值、置信度95%的均值置信区间、置信度95%的两两年差值置信区间的观察，此指数呈下降态。但是通过假设检验验证，P-value 为0.048，非常接近α常规阈值0.05（如果取α=0.05，便可以拒绝原假设，认为空气有改善）。但是当α取值为0.01时，便没有充分把握得到“空气质量有提升”这个结论。小白充分考虑“疑点利益归于被告”原则，下的结论为：没有充分的把握证明指数为下降态。
part_2. 观察污染物相关关系
1. 观察PM10和PM2.5的相关关系
观察均值情况
# PM10/PM2.5年度均值变化
data.groupby(['year']).mean()
年度污染物均值变化
# PM10/PM2.5均值箱图
data.boxplot(['PM10','PM25'], by = 'year', figsize = (20,6))# PM10/PM2.5均值箱图
data.boxplot(['PM10','PM25'], by = 'year', figsize = (20,6))
PM10/PM2.5均值箱图
# 每天 及 周平均 PM10/PM2.5
data.plot(y =['PM10','PM25'],figsize=(20, 6))
plt.title('BY DADY PM10/PM2.5')
data.resample('W').mean().plot(y =['PM10','PM25'],figsize=(20, 6))
plt.title('BY WEEK PM10/PM2.5')
每天 及 周平均 PM10/PM2.5
计算相关性
# PM10 和 PM2.5 相关性
# 绘制 PM10 和 PM2.5 的散点图
df_data.plot.scatter(x='PM10', y='PM25',alpha = 0.5,edgecolors= 'w',figsize=(6, 6))
# 添加趋势线
x = df_data['PM10']
y = df_data['PM25']
z = np.polyfit(x, y, 1)
p = np.poly1d(z)
plt.plot(x,p(x),'r',alpha = 0.5)
# 计算相关系数
co = np.corrcoef(df_data.PM10, df_data.PM25)[0,1]
print 'pm10 vs pm2.5_corrcoef: ',co
-----------------------------------------------------------------------------------------------------------------------------------------------
# 输出结果
pm10 vs pm2.5_corrcoef: 0.876348133246
PM10 和 PM2.5 的散点图
确认是否具有显著的相关性
# 定义函数
def draw_bs_pairs(x, y, func, size=1):
"""对配对数据使用bootstrap方法"""
# 设置索引: inds
inds = np.arange(len(x))
# 初始化空数组
bs_replicates = np.empty(size)
# 重复抽样
for i in range(size):
bs_inds = np.random.choice(inds, len(inds))
bs_x, bs_y = x[bs_inds], y[bs_inds]
bs_replicates[i] = func(bs_x, bs_y)
return bs_replicates
def pearson_r(x, y):
"""计算皮尔逊相关系数"""
corr_mat = np.corrcoef(x,y)
return corr_mat[0,1]
# PM10与PM2.5显著性检验
# 获取1000组重新抽样得到的相关系数
bs_replicates = draw_bs_pairs(df_data.PM10, df_data.PM25, pearson_r, 10000)
# 计算95%置信区间
conf_int = np.percentile(bs_replicates, [2.5, 97.5])
# 输出相关系数的95%置信区间
print(conf_int)
-----------------------------------------------------------------------------------------------------------------------------------------------
# 输出结果
[ 0.83981491 0.90752951]
看起来PM10和PM2.5之间是有比较显著的相关性的。那么接下来再看看这两个污染物和AQI之间是不是分别独立相关。
计算PM类污染物与AQI相关关系
# AQI 与 PM10/PM2.5 相关性
# 绘制AQI和PM10/PM2.5的散点图
ax =df_data.plot.scatter(x='AQI', y='PM10',edgecolors= 'w',label='PM10')
df_data.plot.scatter(x='AQI', y='PM25',color = 'r',alpha = 0.5,edgecolors= 'w',label='PM2.5',figsize=(10,10),ax = ax)
plt.xlabel('AQI')
plt.ylabel('PM10/PM2.5')
# 添加趋势线
x = df_data['AQI']
y1 = df_data['PM10']
z = np.polyfit(x, y1, 1)
p = np.poly1d(z)
plt.plot(x,p(x),'b',alpha = 0.5)
y2 = df_data['PM25']
z = np.polyfit(x, y2, 1)
p = np.poly1d(z)
plt.plot(x,p(x),'r',alpha = 0.5)
# 计算相关系数
pm10_co = np.corrcoef(df_data.AQI, df_data.PM10)[0,1]
pm25_co = np.corrcoef(df_data.AQI, df_data.PM25)[0,1]
print 'AQI vs pm10_corrcoef: ',pm10_co
print 'AQI vs pm2.5_corrcoef: ',pm25_co
-----------------------------------------------------------------------------------------------------------------------------------------------
# 输出结果
AQI vs pm10_corrcoef: 0.829003244255
AQI vs pm2.5_corrcoef: 0.873771297298
AQI和PM10/PM2.5的散点图
确认是否具有显著的相关性
# AQI与PM10显著性检验
bs_replicates_AQI_PM10 = draw_bs_pairs(df_data.AQI, df_data.PM10, pearson_r, 10000)
conf_int = np.percentile(bs_replicates_AQI_PM10, [2.5, 97.5])
print(conf_int)
-----------------------------------------------------------------------------------------------------------------------------------------------
# 输出结果
[ 0.79049242 0.86204875]
# AQI与PM2.5显著性检验
bs_replicates_AQI_PM25 = draw_bs_pairs(df_data.AQI, df_data.PM25, pearson_r, 10000)
conf_int = np.percentile(bs_replicates_AQI_PM25, [2.5, 97.5])
print(conf_int)
-----------------------------------------------------------------------------------------------------------------------------------------------
# 输出结果
[ 0.8470402 0.89783092]
PM10和PM2.5的相关关系小结
通过观察PM10和PM2.5的均值，计算相关系数以及确认是否为显著性相关，确认PM10和PM2.5有很强的相关关系，可谓相辅相成。观察PM类污染物与AQI的相关性，发现也是具有强相关的。
2. 观察排放相关的污染物的相关关系
观察均值情况
# 排放相关污染物均值箱图
data.boxplot(['CO'], by = 'year',figsize = (10,6))
data.boxplot(['NO2'], by = 'year',figsize = (10,6))
data.boxplot(['O3'], by = 'year',figsize = (10,6))
data.boxplot(['SO2'], by = 'year',figsize = (10,6))
CO Boxplot
NO2 Boxplot
O3 Boxplot
SO2 Boxplot
# 每天 及 周平均 排放相关污染物
data.plot(y =['CO'],figsize=(20, 6))
plt.title('BY DADY CO')
data.plot(y =['NO2'],figsize=(20, 6))
plt.title('BY DADY NO2')
data.plot(y =['O3'],figsize=(20, 6))
plt.title('BY DADY O3')
data.plot(y =['SO2'],figsize=(20, 6))
plt.title('BY DADY SO2')
BY DADY CO
BY DADY NO2
BY DADY O3
BY DADY SO2
通过观察四个与排放相关的污染物，可以发现除了臭氧（O3）以外，都呈现比较明显的季节性，冬季（11-3月）污染物多，秋季（9-10月）污染物相对较少。
计算相关性
# 排放相关污染物相关性
# 绘制 四项污染物 的散点图
ax = df_data.plot.scatter(x='AQI', y='CO',alpha = 0.5,edgecolors= 'w',label='CO',figsize=(10,10))
df_data.plot.scatter(x='AQI', y='NO2',color = 'r',alpha = 0.5,edgecolors= 'w',label='NO2',ax = ax)
df_data.plot.scatter(x='AQI', y='O3',color = 'g',alpha = 0.5,edgecolors= 'w',label='O3',ax = ax)
df_data.plot.scatter(x='AQI', y='SO2',color = 'y',alpha = 0.5,edgecolors= 'w',label='SO2',ax = ax)
# 添加趋势线
x = df_data['AQI']
y1 = df_data['CO']
z = np.polyfit(x, y1, 1)
p = np.poly1d(z)
plt.plot(x,p(x),'b',alpha = 0.5)
y2 = df_data['NO2']
z = np.polyfit(x, y2, 1)
p = np.poly1d(z)
plt.plot(x,p(x),'r',alpha = 0.5)
y3 = df_data['O3']
z = np.polyfit(x, y3, 1)
p = np.poly1d(z)
plt.plot(x,p(x),'g',alpha = 0.5)
y4 = df_data['SO2']
z = np.polyfit(x, y4, 1)
p = np.poly1d(z)
plt.plot(x,p(x),'y',alpha = 0.5)
# 计算相关系数
co_CO = np.corrcoef(df_data.AQI, df_data.CO)[0,1]
co_NO2 = np.corrcoef(df_data.AQI, df_data.NO2)[0,1]
co_O3 = np.corrcoef(df_data.AQI, df_data.O3)[0,1]
co_SO2 = np.corrcoef(df_data.AQI, df_data.SO2)[0,1]
print 'AQI vs CO_corrcoef: ',co_CO
print 'AQI vs NO2_corrcoef: ',co_NO2
print 'AQI vs O3_corrcoef: ',co_O3
print 'AQI vs SO2_corrcoef: ',co_SO2
-----------------------------------------------------------------------------------------------------------------------------------------------
# 输出结果
AQI vs CO_corrcoef: 0.735547079502
AQI vs NO2_corrcoef: 0.644245282007
AQI vs O3_corrcoef: -0.217972472764
AQI vs SO2_corrcoef: 0.483738644427
四项污染物散点图
确认是否具有显著的相关性
# AQI与CO显著性检验
bs_replicates_AQI_CO = draw_bs_pairs(df_data.AQI, df_data.CO, pearson_r, 10000)
conf_int = np.percentile(bs_replicates_AQI_CO, [2.5, 97.5])
print(conf_int)
-----------------------------------------------------------------------------------------------------------------------------------------------
# 输出结果
[ 0.69157415 0.77446407]
# AQI与NO2显著性检验
bs_replicates_AQI_NO2 = draw_bs_pairs(df_data.AQI, df_data.NO2, pearson_r, 10000)
conf_int = np.percentile(bs_replicates_AQI_NO2, [2.5, 97.5])
print(conf_int)
-----------------------------------------------------------------------------------------------------------------------------------------------
# 输出结果
[ 0.59098479 0.69007658]
# AQI与O3显著性检验
bs_replicates_AQI_O3 = draw_bs_pairs(df_data.AQI, df_data.O3, pearson_r, 10000)
conf_int = np.percentile(bs_replicates_AQI_O3, [2.5, 97.5])
print(conf_int)
-----------------------------------------------------------------------------------------------------------------------------------------------
# 输出结果
[-0.2687816 -0.16408575]
# AQI与SO2显著性检验
bs_replicates_AQI_SO2 = draw_bs_pairs(df_data.AQI, df_data.SO2, pearson_r, 10000)
conf_int = np.percentile(bs_replicates_AQI_SO2, [2.5, 97.5])
print(conf_int)
-----------------------------------------------------------------------------------------------------------------------------------------------
# 输出结果
[ 0.4285412 0.53974655]
看起来这四个和排放相关的污染物和AQI都没有特别显著的相关性，甚至臭氧（O3）还存在着一定程度的负相关。
CO/NO2/O3/SO2的相关关系小结
通过观察四项污染物的均值，发现除了臭氧（O3）外均呈现正相关，而臭氧（O3）与其他三个污染物呈现负相关。计算相关系数以及确认是否为显著性相关， 得出结论与观察均值结论相同。
小白基于臭氧（O3）的特殊形态进行了资料查阅，发现臭氧其实是一个化学反应剂：
O3 + SO2 = SO3 + O2
O3 + CO = CO2 + O2
O3 + NO = NO2 + O2
由此终于可以解释为何臭氧高的时候二氧化硫（SO2）和一氧化碳（CO）都会低了，原来是起了化学反应，被转化为三氧化硫（SO3）和二氧化碳（CO2）了。小白又顺便查了一下三氧化硫（SO3），它居然是酸雨的主要来源之一。
OMG，以后不可以再装酷淋雨走了，再这样下去估计没两天小白就要告别黑长直变身“地中海”了（T.T）
发现了真相的小白，已哭晕在电脑前……
