分析印度德里封锁对空气质量的影响#
您将做什么#
计算空气质量指数 (AQI) 并对其进行配对学生 t 检验。
您将学到什么#
您将学习移动平均值的概念
您将学习如何计算空气质量指数 (AQI)
您将学习如何执行配对学生 t 检验并找到
t和p值您将学习如何解释这些值
您需要什么#
在您的环境中安装了 SciPy
对统计术语的基本理解，例如总体、样本、均值、标准差等。
空气污染问题#
空气污染是我们面临的最主要的污染类型之一，它对我们的日常生活有直接影响。COVID-19 大流行导致世界各地实施封锁；这提供了一个难得的机会来研究人类活动（或缺乏活动）对空气污染的影响。在本教程中，我们将研究德里（受空气污染影响最严重的城市之一）在 2020 年 3 月至 6 月封锁前后和期间的空气质量。为此，我们将首先根据收集到的污染物测量结果计算每小时的空气质量指数。接下来，我们将对这些指数进行采样，并对它们进行 配对学生 t 检验。这将从统计学上向我们表明，由于封锁，空气质量得到了改善，从而支持了我们的直觉。
让我们从将必要的库导入到我们的环境开始。
import numpy as np
from numpy.random import default_rng
from scipy import stats
构建数据集#
我们将使用 印度空气质量数据 数据集的精简版本。该数据集包含印度多个城市多个站点的小时和每日空气质量数据和 AQI（空气质量指数）。本教程中提供的精简版本包含 2019 年 5 月 31 日至 2020 年 6 月 30 日德里的每小时污染物测量数据。它包含空气质量指数计算所需的标准污染物的测量数据以及其他一些重要的污染物：颗粒物（PM 2.5 和 PM 10）、二氧化氮 (NO2)、氨 (NH3)、二氧化硫 (SO2)、一氧化碳 (CO)、臭氧 (O3)、氮氧化物 (NOx)、一氧化氮 (NO)、苯、甲苯和二甲苯。
让我们打印出前几行，以便大致了解我们的数据集。
! head air-quality-data.csv
Datetime,PM2.5,PM10,NO2,NH3,SO2,CO,O3,NOx,NO,Benzene,Toluene,Xylene
2019-05-31 00:00:00,103.26,305.46,94.71,31.43,30.16,3.0,18.06,178.31,152.73,13.65,83.47,2.54
2019-05-31 01:00:00,104.47,309.14,74.66,34.08,27.02,1.69,18.65,106.5,79.98,11.35,76.79,2.91
2019-05-31 02:00:00,90.0,314.02,48.11,32.6,18.12,0.83,28.27,48.45,25.27,5.66,32.91,1.59
2019-05-31 03:00:00,78.01,356.14,45.45,30.21,16.78,0.79,27.47,44.22,21.5,3.6,21.41,0.78
2019-05-31 04:00:00,80.19,372.9,45.23,28.68,16.41,0.76,26.92,44.06,22.15,4.5,23.39,0.62
2019-05-31 05:00:00,83.59,389.97,39.49,27.71,17.42,0.76,28.71,39.33,21.04,3.25,23.59,0.56
2019-05-31 06:00:00,79.04,371.64,39.61,26.87,16.91,0.84,29.26,43.11,24.37,3.12,15.27,0.46
2019-05-31 07:00:00,77.32,361.88,42.63,27.26,17.86,0.96,27.07,48.22,28.81,3.32,14.42,0.41
2019-05-31 08:00:00,84.3,377.77,42.49,28.41,20.19,0.98,33.05,48.22,27.76,3.4,14.53,0.4
在本教程中，我们只关心计算 AQI 所需的标准污染物，即 PM 2.5、PM 10、NO2、NH3、SO2、CO 和 O3。因此，我们只使用 np.loadtxt 导入这些特定列。然后我们将 切片 并创建两个集合：pollutants_A 包含 PM 2.5、PM 10、NO2、NH3 和 SO2，pollutants_B 包含 CO 和 O3。这两个集合的处理方式略有不同，我们稍后会看到。
pollutant_data = np.loadtxt("air-quality-data.csv", dtype=float, delimiter=",",
skiprows=1, usecols=range(1, 8))
pollutants_A = pollutant_data[:, 0:5]
pollutants_B = pollutant_data[:, 5:]
print(pollutants_A.shape)
print(pollutants_B.shape)
(9528, 5)
(9528, 2)
我们的数据集可能包含缺失值，用 NaN 表示，因此让我们使用 np.isfinite 快速检查一下。
np.all(np.isfinite(pollutant_data))
np.True_
这样，我们就成功导入了数据并检查了数据的完整性。让我们继续进行 AQI 计算！
计算空气质量指数#
我们将使用印度 中央污染控制委员会 采用的 方法 来计算 AQI。总结步骤如下：
收集标准污染物的 24 小时平均浓度值；CO 和 O3 为 8 小时。
使用以下公式计算这些污染物的子指数：
\[ Ip = \dfrac{\text{IHi – ILo}}{\text{BPHi – BPLo}}\cdot{\text{Cp – BPLo}} + \text{ILo} \]其中，
Ip= 污染物p的子指数
Cp= 污染物p的平均浓度
BPHi= 浓度断点，即大于或等于Cp
BPLo= 浓度断点，即小于或等于Cp
IHi= 与BPHi相对应的 AQI 值
ILo= 与BPLo相对应的 AQI 值任何给定时间的最大子指数即为空气质量指数。
空气质量指数是借助下图所示的断点范围计算的。
让我们创建两个数组来存储 AQI 范围和断点，以便我们以后可以将它们用于计算。
AQI = np.array([0, 51, 101, 201, 301, 401, 501])
breakpoints = {
'PM2.5': np.array([0, 31, 61, 91, 121, 251]),
'PM10': np.array([0, 51, 101, 251, 351, 431]),
'NO2': np.array([0, 41, 81, 181, 281, 401]),
'NH3': np.array([0, 201, 401, 801, 1201, 1801]),
'SO2': np.array([0, 41, 81, 381, 801, 1601]),
'CO': np.array([0, 1.1, 2.1, 10.1, 17.1, 35]),
'O3': np.array([0, 51, 101, 169, 209, 749])
}
移动平均值#
第一步，我们必须计算 移动平均值，pollutants_A 的窗口为 24 小时，pollutants_B 的窗口为 8 小时。我们将使用 np.cumsum 和 切片索引 来编写一个简单的函数 moving_mean 来实现这一点。
为了确保两个集合的长度相同，我们将根据 pollutants_A_24hr_avg 的长度截断 pollutants_B_8hr_avg。这也将确保我们在同一时间段内拥有所有污染物的浓度。
def moving_mean(a, n):
ret = np.cumsum(a, dtype=float, axis=0)
ret[n:] = ret[n:] - ret[:-n]
return ret[n - 1:] / n
pollutants_A_24hr_avg = moving_mean(pollutants_A, 24)
pollutants_B_8hr_avg = moving_mean(pollutants_B, 8)[-(pollutants_A_24hr_avg.shape[0]):]
现在，我们可以使用 np.concatenate 将两个集合连接起来，形成所有平均浓度的单个数据集。请注意，我们必须按列连接数组，因此我们传递 axis=1 参数。
pollutants = np.concatenate((pollutants_A_24hr_avg, pollutants_B_8hr_avg), axis=1)
子指数#
每个污染物的子指数是根据 AQI 与标准断点范围之间的线性关系计算的，公式如上所示。
compute_indices 函数首先在上面创建的数组 AQI 和 breakpoints 的帮助下，获取输入浓度和污染物对应的 AQI 类别和断点浓度的正确上限和下限。然后，它将这些值输入公式以计算子指数。
def compute_indices(pol, con):
bp = breakpoints[pol]
if pol == 'CO':
inc = 0.1
else:
inc = 1
if bp[0] <= con < bp[1]:
Bl = bp[0]
Bh = bp[1] - inc
Ih = AQI[1] - inc
Il = AQI[0]
elif bp[1] <= con < bp[2]:
Bl = bp[1]
Bh = bp[2] - inc
Ih = AQI[2] - inc
Il = AQI[1]
elif bp[2] <= con < bp[3]:
Bl = bp[2]
Bh = bp[3] - inc
Ih = AQI[3] - inc
Il = AQI[2]
elif bp[3] <= con < bp[4]:
Bl = bp[3]
Bh = bp[4] - inc
Ih = AQI[4] - inc
Il = AQI[3]
elif bp[4] <= con < bp[5]:
Bl = bp[4]
Bh = bp[5] - inc
Ih = AQI[5] - inc
Il = AQI[4]
elif bp[5] <= con:
Bl = bp[5]
Bh = bp[5] + bp[4] - (2 * inc)
Ih = AQI[6]
Il = AQI[5]
else:
print("Concentration out of range!")
return ((Ih - Il) / (Bh - Bl)) * (con - Bl) + Il
我们将使用 np.vectorize 来利用向量化的概念。这仅仅意味着我们不必自己遍历污染物数组的每个元素。向量化 是 NumPy 的主要优势之一。
vcompute_indices = np.vectorize(compute_indices)
通过对每种污染物调用我们的矢量化函数vcompute_indices，我们可以得到子索引。为了得到一个具有原始形状的数组，我们使用np.stack。
sub_indices = np.stack((vcompute_indices('PM2.5', pollutants[..., 0]),
vcompute_indices('PM10', pollutants[..., 1]),
vcompute_indices('NO2', pollutants[..., 2]),
vcompute_indices('NH3', pollutants[..., 3]),
vcompute_indices('SO2', pollutants[..., 4]),
vcompute_indices('CO', pollutants[..., 5]),
vcompute_indices('O3', pollutants[..., 6])), axis=1)
空气质量指数#
使用np.max，我们找出每个周期的最大子指数，这就是我们的空气质量指数！
aqi_array = np.max(sub_indices, axis=1)
这样，我们就得到了从2019年6月1日至2020年6月30日的每小时AQI。请注意，即使我们从5月31日的数据开始，在移动平均步骤中我们也截断了这些数据。
配对学生t检验在AQI上#
假设检验是一种描述性统计方法，可以帮助我们根据数据做出决策。根据计算得到的AQI数据，我们想要找出封锁前后平均AQI是否存在统计学上的显著差异。我们将使用左侧检验的配对学生t检验来计算两个检验统计量——t 统计量和p 值。然后，我们将把这些值与相应的临界值进行比较，以做出决策。
抽样#
我们现在将把原始数据集中的datetime列导入到一个datetime64数据类型数组中。我们将使用这个数组来索引AQI数组并获得数据集的子集。
datetime = np.loadtxt("air-quality-data.csv", dtype='M8[h]', delimiter=",",
skiprows=1, usecols=(0, ))[-(pollutants_A_24hr_avg.shape[0]):]
由于德里从2020年3月24日开始全面封锁，因此封锁后的子集为2020年3月24日至2020年6月30日。封锁前的子集是3月24日之前的同等时间段。
after_lock = aqi_array[np.where(datetime >= np.datetime64('2020-03-24T00'))]
before_lock = aqi_array[np.where(datetime <= np.datetime64('2020-03-21T00'))][-(after_lock.shape[0]):]
print(after_lock.shape)
print(before_lock.shape)
(2376,)
(2376,)
为了确保我们的样本近似服从正态分布，我们取大小为n = 30的样本。before_sample和after_sample是在全面封锁前后抽取的一组随机观测值。我们使用random.Generator.choice来生成样本。
rng = default_rng()
before_sample = rng.choice(before_lock, size=30, replace=False)
after_sample = rng.choice(after_lock, size=30, replace=False)
定义假设#
让我们假设封锁前后样本均值之间没有显著差异。这将是零假设。备择假设是均值之间存在显著差异，并且AQI有所改善。数学上表示为：
\(H_{0}: \mu_\text{after-before} = 0\)
\(H_{a}: \mu_\text{after-before} < 0\)
计算检验统计量#
我们将使用t统计量来评估我们的假设，甚至根据它计算p 值。t统计量的公式为：
其中：
\(\mu_\text{after-before}\) = 样本均值差异
\(\sigma^{2}\) = 均值差异的方差
\(n\) = 样本大小
def t_test(x, y):
diff = y - x
var = np.var(diff, ddof=1)
num = np.mean(diff)
denom = np.sqrt(var / len(x))
return np.divide(num, denom)
t_value = t_test(before_sample, after_sample)
对于p值，我们将使用SciPy的stats.distributions.t.cdf()函数。它接受两个参数——t 统计量和自由度(dof)。dof的公式为n - 1。
dof = len(before_sample) - 1
p_value = stats.distributions.t.cdf(t_value, dof)
print("The t value is {} and the p value is {}.".format(t_value, p_value))
The t value is -6.778550885412477 and the p value is 9.648125019613442e-08.
t和p值意味着什么？#
我们现在将计算得到的检验统计量与临界检验统计量进行比较。临界t值可以通过查阅t分布表来计算。
从上表中，对于29个dof，在95%的置信水平下，临界值为1.699。由于我们使用的是左侧检验，因此我们的临界值为-1.699。显然，计算得到的t值小于临界值，因此我们可以安全地拒绝零假设。
临界p值，用\(\alpha\)表示，通常选择为0.05，对应于95%的置信水平。如果计算得到的p值小于\(\alpha\)，则可以安全地拒绝零假设。显然，我们的p值远小于\(\alpha\)，因此我们可以拒绝零假设。
这并不意味着我们可以接受备择假设。它只告诉我们没有足够的证据来拒绝\(H_{a}\)。换句话说，我们未能拒绝备择假设，因此它可能是正确的。
在实践中…#
对于时间序列数据分析，最好使用pandas库。
SciPy stats模块提供了stats.ttest_rel函数，可用于获取
t 统计量和p 值。在现实生活中，数据通常不服从正态分布。对于此类非正态数据，存在诸如Wilcoxon检验之类的检验。