分析印度德里封锁对空气质量的影响#
您将做什么#
计算空气质量指数 (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 日的数据开始，我们也在移动平均步骤中对其进行了截断。
对 AQI 进行配对学生 t 检验#
假设检验是一种描述性统计方法，用于帮助我们根据数据做出决策。从计算得到的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 -10.97328109246125 and the p value is 3.8517128805721366e-12.
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检验这样的检验方法。